Spatiotemporal Structure of Cortical Activity: Properties and Behavioral Relevance

Yifat Prut1, 2, Eilon Vaadia1, Hagai Bergman1, Iris Haalman1, Hamutal Slovin1, and Moshe Abeles1

1 Department of Physiology, School of Medicine and the Interdisciplinary Center for Neural Computation, The Hebrew University of Jerusalem, Jerusalem 91120, Israel; and 2 Regional Primate Research Center, University of Washington, Seattle, Washington 98195

    ABSTRACT
Abstract
Introduction
Methods
Results
Discussion
References

Prut, Yifat, Eilon Vaadia, Hagai Bergman, Iris Haalman, Hamutal Slovin, and Moshe Abeles. Spatiotemporal structure of cortical activity: properties and behavioral relevance. J. Neurophysiol. 79: 2857-2874, 1998. The study was designed to reveal occurrences of precise firing sequences (PFSs) in cortical activity and to test their behavioral relevance. Two monkeys were trained to perform a delayed-response paradigm and to open puzzle boxes. Extracellular activity was recorded from neurons in premotor and prefrontal areas with an array of six microelectrodes. An algorithm was developed to detect PFSs, defined as a set of three spikes and two intervals with a precision of ±1 ms repeating significantly more than expected by chance. The expected level of repetition was computed based on the firing rate and the pairwise correlation of the participating units, assuming a Poisson distribution of event counts. Accordingly, the search for PFSs was corrected for rate modulations. PFSs were found in 24/25 recording sessions. Most PFSs (76%) were composed of spikes of more than one unit but usually not more than two units (67%). The PFSs spanned hundreds of milliseconds, and the average interval between two events within the PFSs was 200 ms. No traces of periodic oscillations were found in the PFS intervals. The bins of the matrix that were defined as PFSs were isolated temporally: the spikes that generated PFSs were not associated with high-frequency bursts or rapid coherent rate fluctuations. A given PFS tended to be correlated with the animal's behavior. Furthermore, for 19% of the PFS pairs that shared the same unit composition, each member of the pair was associated with a different type of behavior. The PFSs often appeared in clusters that were associated with particular phases of the behavior. The firing rate of single units did not provide a full explanation for the timing and structure of these clusters. A reduced spike train (RST) was defined for each unit by taking all spikes of that unit that were part of any PFS. In 88% of the cases the degree of modulation of the RST was higher than that of the complete spike train. The results suggest that relevant information is carried by the fine temporal structure of cortical activity. A coding scheme that involves such temporal structures is rich and sufficiently flexible to facilitate a rapid organization of cortical neurons into functional groups. The results can be accounted for by the synfire chain model, which suggests that cortical activity is mediated by synchronous activation of neural groups in a reverberatory mode.

    INTRODUCTION
Abstract
Introduction
Methods
Results
Discussion
References

It generally is accepted that mechanisms of neural coding are based on the spatiotemporal organization of cortical activity. However, the appropriate resolution both in space and in time for such a description is not yet clear. The present work tests the idea that cortical activity is carried in a network that relies on the single neuron resolution and a few millisecond time scale.

The single neuron time-dependent rate function was taken by many as the coding parameter (e.g., Barlow 1972, 1992, 1994; Newsome et al. 1989; Rolls 1991). Others suggested a population coding, based on either the summed activity of neurons (Georgopoulos et al. 1986; Schwartz 1994), or the coherency in firing among cells (Eckhorn et al. 1988; Engel et al. 1991a-c; Gray and Singer 1992; Gray et al. 1989, 1992). Both views ignored the detailed temporal structure of cortical activity, assuming that precision is not compatible with a noisy cortical environment. Despite this notion, several studies further explored the fine temporal structure of neuronal activity at the single spike resolution (Bialek et al. 1991; Eckhorn and Popel 1974; Eckhorn et al. 1976; Steveninck et al. 1997). These also include the study of preferred patterns within the spike train (Dayhoff and Gerstein 1983a,b; Lestienne and Strehler 1987), neural bursts (Bair et al. 1994; Legendy and Salcman 1985), and single-unit oscillatory activity (Kreiter and Singer 1992).

The present study searches for precise patterns of neural activity. The motivation for this search is the synfire chain (SFC) model. This model exploits the presumed sensitivity of cortical neurons to synchronous activation shown both theoretically (Abeles 1982a,b; Nelken 1988) and experimentally (Douglas et al. 1991), the divergent/convergent mode of the cortical network, and finally, the Hebbian synaptic plasticity shown in slices (Alonso et al. 1990; Kirkwood and Bear 1994; Markram et al. 1997) and, more implicitly, in the behaving animal (Ahissar et al. 1992; Ghose et al. 1994).

According to the model, activity in the cortex is naturally organized into groups of cells (pools). The pools are interconnected with a divergent and convergent connections, creating neural chains. Each pool is activated synchronously and activates synchronously the next pool down the chain. The wiring parameters known for the cortex, and its assumed synaptic strength, are both in accordance with the existence of such a network in a way that ensures a secured flow of activity. However, the actual existence of SFCs needs an explicit confirmation.

The model offers a testable prediction: if cortical SFCs exist, simultaneous recording of several cells may provide few cells that are part of different pools composing one chain. In such a case, every time the chain is activated, the participating cells will fire spikes in a sequential manner. The temporal spacing between the spikes will correspond to the relative location of the parent cells along the chain. Such a sequence of spikes is termed a precise firing sequence (PFS).

Several studies, which pursued the SFC model using a dedicated algorithm (Abeles and Gerstein 1988), detected PFSs in cortical activity of the frontal lobe (Abeles et al. 1993; Villa and Fuster 1992) and in the auditory thalamus (Villa and Abeles 1990). Most frequently PFSs found in these studies were composed of spikes of a single unit. Hence it was concluded that the SFC model is not a pure feedforward model but contains feedback connections, creating a reverberatory circuit in which activity, once triggered, can sustain itself for a long period of time. However, these studies were criticized on the basis that neither study established a robust correlation between the precise firing sequences that were found and the behavior of the animal (Shadlen and Newsome 1994). To test the model while paying special attention to its behavioral relevance, recordings were made in areas located in the frontal cortex of awake behaving monkeys. This part of the cortex was selected on the basis of its documented integrative function (Butter and Synder 1972; Damasio 1979; Fuster 1989; Goldman-Rakic 1988; Mishkin 1964; Nauta 1971). It was shown to be involved in planning of movement using "working memory" information, self-initiated movements (Fuster 1991), and divergent thinking skills that require the selection of one solution out of many to a given problem (Milner and Petrides 1984). This area includes cells that are engaged actively in delayed-response tasks (Weinrich and Wise 1982; Weinrich et al. 1984; Wise and Mauritz 1985), delayed-alternation tasks (Goldman and Rosvold 1970), and preparation for movement (Passingham 1988, 1993).

    METHODS
Abstract
Introduction
Methods
Results
Discussion
References

Behavioral paradigm

Two 3-4 kg female rhesus monkeys (Macaca mulatta) were trained to perform a delayed-response paradigm, and one of them additionally was trained to perform a puzzle box opening paradigm.

DELAYED-RESPONSE (DR) PARADIGM. The monkey initiated a DR trial by pressing a central Ready key (Fig. 1A), an event that triggered the onset of a central red light (Ready signal). After3-6 s, a 200-ms visual cue was delivered, in either the left or the right visual field. After a delay period of 1, 2, 4, 8, 16, or 32 s (selected at random), the Ready light was dimmed. This signal (Go signal) instructed the monkey to respond within a given reaction time (RT < 0.6 s) by touching the key from which the cue was emitted. The allowed movement time also was limited to 0.6 s. Correct responses (both temporally and spatially) were reinforced by a drop of juice (Reward). Such a sequence of events, from the Ready light onset, until the reward was given (or earlier in the case of an incorrect performance) was a trial in a Go mode. The monkey also was trained to perform trials in a No-Go mode. In this mode, the sequence of events was the same as in the Go mode, but after the Go signal, the monkey had to hold its hand on the central key for 1.2 s. A set of lights, which were turned on for 3-4 s, instructed the monkey to switch modes (Go to No-Go or No-Go to Go). This signal was given after the monkey performed four correct trials in each mode. Recording sessions at this paradigm typically lasted 3-4 h. For data analysis, we defined six different, nonoverlapping, intervals of correct trials: pre-cue Go or No-Go: the period from the Ready light to the visual cue in the Go (or No-Go) mode; Go delay left or right: the period from the left (right) visual cue to the Go signal in the Go mode; and No-Go delay left or right: the period from the left (right) visual cue to the Go signal in the No-Go mode.


View larger version (26K):
[in this window]
[in a new window]
 
FIG. 1. Behavioral paradigms. A: schematic illustration of the sequence of events and signals that composed the delayed response paradigm. Trial was initiated when the monkey pressed the central key, an event that triggered the onset of a red light (Ready). After a pre-cue period of 3-6 s, a light (Visual-Cue) was turned on for 200 ms from either a left or a right target key. After a delay period of 1-32 s, a Go signal was given to the monkey. According to the behavioral mode, the monkey either had to touch the target key (in the Go mode) or hold its hand on the central key (in the No-Go mode). Correct performance was reinforced by a reward. After the performance of 4 correct trials in a given mode, a switching mode signal arrived, instructing the monkey to switch the behavioral mode either to Go or to No-Go, depending on the current mode. B: examples of 2 of the 9 boxes that were used in the puzzle box opening paradigm. Each box had a unique opening configuration (the arrows in the figure indicate the specific opening mechanism of the illustrated boxes). External surfaces of all boxes were identical in texture. C: schematic illustration of the event sequences composing the puzzle box opening paradigm. Box was presented to the monkey (Serve), which tried to open it (Touch). After the box was opened (Open), the monkey ate the reward inside (Eat), and then the box was removed (Remove). No time restrictions were applied in these trials.

Note that during the six intervals, the monkey's behavior is identical (resting its palm on the central key), whereas the behavioral context of these intervals is different.

PUZZLE BOX (PB) PARADIGM. Nine boxes, identical in shape but having different opening mechanisms (Fig. 1B) were used. The external faces of the boxes were etched in a grid manner, to avoid tactile recognition of the opening strategy. A box was selected randomly and presented to the monkey, which had to open it to obtain the reward inside (a piece of apple). The sequence of events along a trial is shown in Fig. 1C. Two intervals of the trials were used for analysis: pre-trial: the period before box presentation and trial: the period starting at box presentation until box opening.

The monkey sat in a soundproof, dimmed room, and its behavior was monitored and recorded through an infrared video camera. The events during a trial were keyed manually into a computer. Recording sessions at this paradigm lasted ~2 h.

Surgery

After the monkeys were trained fully, surgery was performed in aseptic conditions. Anesthesia was induced by ketamin injection (8 mg/kg im), followed by doses of sodium pentobarbital (accumulative amount of 30 mg/kg iv). The skull above the frontal lobe was exposed, and a hole was trephined. Stainless steel recording chambers were centered over the superior limb of the Arcuate sulcus of either one or both hemispheres. The chambers were cemented in place with dental acrylic. A metal rod was attached to the occipital bone to enable fixation of the monkey's head during the recording sessions. Two Ag-AgCl cup-electrodes were implanted in the orbital bone to measure horizontal eye movements of the monkey. After recovery and retraining sessions, recording were started. All treatments and maintenance of the monkey adhered to the regulations of the Hebrew University.

Data acquisition setup

Recordings were made using an array of six glass-coated tungsten microelectrodes (impedance of 0.5-1.5 MOmega , measured at 1 kHz). The electrodes were arranged in a circular array, with one electrode located in the center of this circle. The distance between neighboring electrodes was 500 µm with a maximal distance of 1,000 µm. The electrodes were advanced either together or individually by remote-controlled microdrives. Spike-sorting devices enabled recording of <= 16 single units simultaneously. Spiking times (sampled at 1-ms resolution), Electrooculogram and behavioral events were all digitized and stored in files for further analysis.

Data analysis

Previously, an algorithm was developed for a global search for PFSs of any complexity (i.e., number of participating spikes) (Abeles and Gerstein 1988). Here we use a more restricted searching algorithm based on the following definition of PFSs: a PFS is made by three spikes and two intervals. It is defined by its unit composition (i.e., the units that fired the spikes composing the PFS) and the time delays between its spikes. Figure 2 schematically illustrates one PFS; a spike is fired by unit S1, then after t1 ms, a spike is fired by unit S2, and after t2 ms (t2 > t1) after the first spike, a spike is fired by unit S3. This precise firing sequence is designated according to the following convention
(<IT>S</IT><SUB>1</SUB>, <IT>S</IT><SUB>2</SUB>, <IT>S</IT><SUB>3</SUB>; <IT>t</IT><SUB>1</SUB>, <IT>t</IT><SUB>2</SUB>)
Note that additional spikes (of any neuron, including the participating units) may appear within the intervals (t1, t2), as schematically illustrated in Fig. 2. According to this definition, S1, S2, and S3 could represent spikes of different units or the same unit.


View larger version (20K):
[in this window]
[in a new window]
 
FIG. 2. Illustration of a precise firing sequence (PFS). Spike trains of 3 different units are shown. PFS is defined by its unit composition (here S1, S2, and S3) and the time interval between these units (here t1 and t2). Spikes in the figure are presented as rectangles. Among the 3 spikes that constitute the PFS (black-square), other, nonrelated, spikes may appear ().

The first stage in the search for PFSs is constructing a threefold correlation matrix for each triplet of units (1 "trigger" unit and 2 "reference" units). In this matrix, each bin stands for a pair of delays between spikes of the three units. The matrix is computed with a bin size B, and a total length of T, and hence contains N2,[N = (T/B)] bins. This concept of correlation matrix was developed previously by Palm, Aertsen, and Gerstein (Aertsen et al. 1989; Palm et al. 1988) for the joint-peristimulus time histogram (JPSTH) method. However, here instead of using a stimulus as a trigger event, we used a spike emitted by the "trigger" unit as a trigger event. The next step is to locate those bins (i.e., PFSs) in the matrix at which the counts were significantly above chance level. For this, we had to compute a matrix of the expected counts. An explicit formulation of the algorithm to compute the expected matrix is given in the APPENDIX. In short, this computation took into account the pairwise correlation and the firing rate of the participating units. Having both the correlation matrix and the expected count matrix, the probability of obtaining the observed counts in each bin was computed, assuming that the statistic of the counts in each bin obeys Poisson distribution. The case wheret1 t2 (i.e., the bins along the main diagonal) was ignored for reasons described later on (see APPENDIX ). For each such correlation matrix, the probability per bin was used to determine which of its bins deviates significantly from the expected counts.

After recording neural activity from two monkeys, suitable data for further analysis were selected. In total, 40 recording sessions were obtained from the first monkey (B) and 92 from the second monkey (C). For analysis, sessions were selected to fulfill several criteria: 1) stationarity of the data across trials. Stimulus-locked modulations of rate were not taken as signs for nonstationarity. 2) Large number of units in a given recording session. This point emerges from the basic model, which describes a network organization, and thus units' interactions are of extreme interest. Sessions with at least four isolated units were taken for analysis. 3) Well-isolated units were preferred (57 of the 127 units), although units that represent a mixture of activity of several single units (usually 2-3 units) were taken for analysis as well. And 4) cells with relatively high firing rates were preferred to enhance the robustness of the statistical test carried out while searching for PFSs. In most of the cases, such cells demonstrated behaviorally locked modulation of firing rates, which made them more likely candidates for active participation in local network processes.

Altogether 25 sessions fulfilled these criteria (9 DR sessions and 5 PB sessions from monkey B and 11 DR sessions from monkey C). For each triplet of units, a threefold correlation matrix was computed separately for each of the behavioral interval (6 for DR, and 2 for PB, as described above). The time window used for the search (T) was 450 ms, and the bin size (B) was 3 ms. Such a bin size allows a jitter of ±1 ms around the time intervals of the PFSs. This jitter was selected on base of results obtained in a previous study carried in our lab. In this study, it was found that a further increase in the jitter (above ±3 ms up to ±10 ms) resulted in a small increase of new PFS-creating unit triplets, whereas some PFSs ceased to be significant at these bin size values. This result was further validated by others (Lestienne 1996). Last, we were encouraged to keep this low level of jitter by the fact that already in this case many PFSs were found.

Histology

At the completion of the experiments, the monkeys were euthanized with an overdose of pentobarbital (100 mg/kg) and perfused transcardially with normal saline, followed by 4% formaldehyde. The brains then were blocked, frozen, sectioned in the parasagittal plane, and stained with cresyl violet. Recording sites were reconstructed based on the linear gliosis associated with the microelectrode tracks and identified anatomic markers.

    RESULTS
Abstract
Introduction
Methods
Results
Discussion
References

Existence of PFSs in cortical activity

PFSs were found in all but one of the 25 recording sessions analyzed. In 15 sessions (total of 77 cells), >30 different PFSs per session were found. These sessions were used for further analysis. The number of different PFSs per session varied considerably, ranging from several tens to >1,000 different PFSs. The mean number of PFSs per session was 623.5 ± 572.9 SD [a median of 357 and median of absolute deviation from the median (MAD) of 331].

Unit composition of PFSs

The PFSs that were found in the data were divided, based on their unit composition, into three different types (Fig. 3): 3-U PFS composed of spikes emitted by three different units (Fig. 3A); 2-U PFS composed of spikes emitted by two different units (Fig. 3B); and 1-U PFS composed of three spikes emitted by one unit (Fig. 3C). Note that in the figure, the appearance of the PFS spikes as sharp vertical lines is an inevitable outcome of having plotted only those cases of a match between the intraspike intervals and the precise timing criteria of the requested template (i.e., the PFS). The proportion of each of the three different PFS types was computed for each recording session. Then the average proportions across different recording sessions, as well as the weighted-average (obtained by weighing the single day proportions by the total number of PFSs found at that day) were computed (Fig. 3). Most of the PFSs (~76%) were 2-U or 3-U, namely, at least two units participated in their composition. This suggests that a network phenomenon underlie the production of the PFSs. However 67% of the PFSs were 1-U or 2-U, where the same unit took part at least twice in the same PFS. This is consistent with previous results (Abeles et al. 1993).


View larger version (52K):
[in this window]
[in a new window]
 
FIG. 3. Examples of PFS types. A: a 3-unit (3-U) PFS <5,2,6;31,52> composed of spikes fired by 3 different units; B: a2-U PFS <1,4,4;55, 148> composed of spikes fired by two different units; C: a1-U PFS (<13,13,13; 208,313>), composed of spikes emitted by a single unit. Each of the PFSs in these examples was found in the PFS searching process as a significant event. Then a research for their occurrences was carried out through all the data. Each piece of data, containing the appropriate PFS, was aligned so that the 1st spike of the PFS would be at time 0. Vertical straight line at time 0 is composed of the 1st spikes of the PFS. Observed jitter of the 2nd and 3rd spike columns reflects the allowed ±1 jitter in the PFS searching process. Numbers under each of the type names indicates relative proportion of PFS types. For each of the types, both the mean proportion per session and the weighted mean (numbers inside parenthesis) are given. Mean proportion was computed by first computing the percentages of the 3 types in each session and then averaging across sessions. Weighted mean was computed by first pooling together all the PFSs (of all sessions) and then computing the percentages.

When examining the different PFSs that were found (ignoring the number of repetitions of each single PFS), no significant correlation was found between the number of PFSs in which a given unit took part and the unit's firing rate. Nor did we find a significant correlation between the tendency of two units to take part in the same PFS and the distance, the strength of coherence in firing, or the rate of these units. The latter correlations were measured as the slope of the linear regression computed between the probability of two units to be part of a PFS and the tested parameter (e.g., interunit coherence strength, measured by the peak area of the cross-correlation computed for this unit pair). However, one should keep in mind that the algorithm for detecting PFSs corrects for correlation in firing among units. Therefore the results may be biased against PFSs that were generated by units that exhibited coherent firing. As different studies reported that nearby units are more likely to be correlated in their firing (Eggermont 1992; Vaadia et al. 1991), the algorithm then may lead to a similar bias against PFSs that are composed of spikes from neighboring cells.

Temporal structure of PFSs

We studied the temporal characteristics of the timing of spikes within PFSs pursuing the following issues: 1) what are the maximal and minimal delays between spikes within a PFS? 2) Are there any signs of oscillation in the temporal intervals composing a PFS? 3) Are there any preferred intervals found in PFSs between any given pair of units? And 4) What is the temporal precision of PFSs?

These questions were studied in the following way. For a given pair of units {Si, Sj} (where the order of appearance is important), all the PFSs that included these two units, namely (Si, Sj, S), (S, Si, Sj), and (Si, S, Sj) were collected. Then the time interval between Si and Sj was taken from each of these PFSs. Following our definition of PFSs, t1 was taken in case of (Si, Sj, S), t2 - t1 in case of (S, Si, Sj), and t2 in case of (Si, S, Sj), as the appropriate time intervals between the spike of Sj and the spike of Si. The unit pairs {Si, Sj} either contained spikes fired by the same unit (SU pairs where i = j) or by two different units. Note that SU pairs, are not synonymous with 1-U PFSs, as they are found both in 1-U and 2-U PFSs. All by all we had 326 pairs of units (of 492 possible pairs of recorded units) that participated in at least one PFS. Of them, 207 pairs defined >= 10 PFS intervals (e.g., a pair of units which took part in 10 different PFSs was included in this group). This definition considers only the number of different PFS intervals the pair is part of, regardless the number of repetition of each single PFS. These 207 pairs were used for further analysis.

STATISTICS OF THE INTERUNIT, INTRA-PFS TIME INTERVALS. For each group of group of intervals (of the 207 groups), three intervals were computed: the minimal value, the mean value, and the maximal value.

Minimal interval distribution. The resulting distribution had a long (right) tail distribution with a median of 10 ms. This means that for many unit pairs, at least one case (i.e., one PFS) of tight synchrony was found.

Mean interval. This distribution was symmetric and concentrated around an average value of 200 ms.

Maximal interval. The distribution of the maximal values was again skewed, but this time having a left tail. The median value was 436 ms, where 25.6% of the pairs had the maximal interval (448 ms) allowed by the parameters used in the search for PFSs.

RELATION BETWEEN PFSS AND PERIODIC OSCILLATIONS. To trace any signs of neuronal rhythmic activity, we analyzed the distribution of intervals found for SU pairs. Usually these histograms were too sparse for computing their power spectra; therefore they were checked manually and found to contain no traces of oscillations in terms of peaks located around intervals of a specific frequency and its harmonics. Furthermore, the global distribution of the intervals (t1, t2), computed for all the PFSs of all the recording days showed no tendency of any interval to be dominant. It should be mentioned that in all the data (including the 25 currently analyzed sessions), 3.8% of cells were found to demonstrate clear oscillatory activity in the range of 1-2 Hz. This frequency is not expected to be reflected in the structure of PFSs with a maximal time window of 450 ms. We thus conclude that the PFSs in our data are not related to single unit oscillatory activity.

SHAPE OF THE INTERVAL DISTRIBUTION. A distribution of intervals was computed for each of the 207 groups of intervals (where each group corresponds to a pair of units {Si, Sj}). Two examples of these histograms are shown in Fig. 4. The histogram in Fig. 4A illustrates a uniform distribution with no tendency for any specific interval to dominate. The second histogram (Fig. 4B) illustrates a different distribution (computed for a different pair of units), with a clear tendency for some intervals to be dominant. For example, the interval of 76 ms appeared 20 times in PFSs that included the two corresponding units. The existence of dominant intervals was quantified by measuring the normalized coefficient of variation (CV) of each histogram. This measure was obtained by first smoothing the raw histogram (H) by convolving with a Gaussian curve F.1 Then the normalized CV for each histogram was defined as
CV = <FR><NU>σ<SUP>2</SUP>(<IT>H − H</IT>*<IT>F</IT>)</NU><DE><IT>E</IT>(<IT>H</IT>)</DE></FR>
The subtraction emphasizes sharp peaks in the histogram (similar to applying a low-pass filter) and therefore, high CV values are expected for histogram with sharp peaks, whereas flat (or slowly varying) histograms are expected to have low CV values. Indeed in Fig. 4, the flat histogram (Fig. 4A) has a low CV value (0.84), whereas the "peaky" histogram (Fig. 4B) has a high CV value (10.48). The distribution of the CV values that were computed for the 207 different histograms is given in Fig. 4C. The distribution does not decay smoothly with increased CV values. Rather there are excess counts for CV values >4. Accordingly we defined two types of histograms: flat histograms, which contained no dominant intervals (CV < 2), and peaky histograms, which had one or few dominant intervals (CV > 4).


View larger version (12K):
[in this window]
[in a new window]
 
FIG. 4. Characteristics of PFS intervals. A and B: examples of 2 distributions of PFS intervals. We first selected all the PFSs that contained a given pair of units. Then the intervals that were found between the corresponding units in the selected PFSs were pooled, and their distribution was computed. In A (the intervals between units 13 and 11), the resulted histogram is "flat" and its coefficient of variation (CV) value is low (0.84). In B (intervals between units 5 and 12), the histogram is "peaky" and its CV value is high (10.45). In C, the overall distribution of CV values is plotted. It can be seen that the distribution does not decay monotonically at large CV values but rather has an excess of counts for CV values >4. Two dashed vertical lines represent the thresholds selected for defining flat and peaky histograms.

Of the 207 histograms, 89 were flat (CV < 2) and 61 were peaky (CV > 4).

Of the 207 pairs, 40 were SU pairs and the rest (167 pairs) were made of spikes emitted by two different units. In the latter group, 33 pairs were "unidirectional." These were pairs at which the same unit order always was preserved (Si right-arrow Sj), having no intervals of the mirror order (Sj right-arrow Si).

TEMPORAL PRECISION OF PFSS. A way to approach this issue is to evaluate the extent to which PFSs appear as temporally isolated events. This is illustrated in Fig. 5, which presents two extreme alternatives. The first is a case where PFSs ride on top of elevated "hills" in the threefold correlation matrix. These hills correspond to slowly varying rate modulations or higher order correlations. In the alternative situation, PFSs appear as isolated bins in the matrix, emerging from a flat background.


View larger version (25K):
[in this window]
[in a new window]
 
FIG. 5. Illustration presenting 2 alternatives for the temporal preciseness of PFSs. The figure illustrates a cross section through a row (or a column or a diagonal) in a counts-expected correlation matrix. Stepwise plot represents the bins of the matrix. bullet , PFSs. A: case where PFSs are the tips of "hills" of counts excess. These hills have moderate slopes. In this case, the exact location of the threshold determines which bins will be considered as PFSs. Accordingly, the PFSs are part of a much slower phenomenon. B: different case, where the PFSs are isolated peaks. In this case, the PFS bins outstand from their surroundings, being much higher than the background level. Therefore, the exact threshold value has no crucial impact on the PFSs definition, and the PFSs represent a phenomenon of a fast time scale. C: top view of the counts-expected matrix, illustrating the concept of "belts" defined around PFS bins. Each belt, filled with a different gray shading, is 1 bin wide and surrounds the PFS.

The relationships between PFSs and their surroundings were studied in the count-expected matrices. There, for each PFS bin, we defined a set of six square shaped belts (of 1 binwidth) around it (Fig. 5C). Each belt was identified by its distance from the central bin (i.e., the PFS). In case of alternative A in Fig. 5, it is expected to find a general decrease of values as we go further away from the central bin, reflecting the presumed hills on top of which the PFSs are riding. Such a decrease will be manifested in a significant negative slope of the regression line computed between the value of the bins around a PFS, and their distance from it. If alternative B is the correct one, then the regression line should be flat. The analysis, carried over all the recording sessions, yielded zero slopes in all but one case [in this case, negative outliers in the matrix yielded a negative slope; removal of these outliers, based on their leverage coefficient (Sokal and Rohlf 1981), resulted again in a 0 slope]. As this study extended <= 21 ms from the central bin, it indicates that PFSs do not ride on locally elevated hills at that range.

Second, in addition to the PFS bins themselves, which were (as expected) significantly different from the background level, the bins composing the first belt around the PFSs also tended to be elevated compared with the background level (although they were still much lower than the PFS itself). A similar result was reported in a previous study (Abeles et al. 1993). Such a behavior is expected from a phenomenon of a temporal precision in the same order of magnitude as the bin size. Thus the "arbitrary" binning of the data (as done for generating the interval matrix) leads to extra counts at bins adjacent to the PFS bins.

As a control test, the same approach was applied on selected sets of random bins (pseudo-PFSs), arbitrarily defined in these matrices. The results then were compared with those found for real PFSs. None of the above phenomena was found for the pseudo-PFS bins.

Relationships between PFSs and the spike trains

So far we focused on the properties of unit composition and temporal structure of PFSs. Next we studied the relationships between PFSs (which include a fraction of the spikes emitted by any given unit) and the entire spike trains. This approach is useful in extracting information about where along the spike train PFSs appear. To study these relationships, the PFSs, which were detected as significantly repeating events in any of the six behavioral segments, were used as templates. Approximate (±1 ms) matches to these templates were searched all over the data. Then for each unit, all the spikes emitted by this unit, which were part of any PFSs, were marked as PFS spikes. Three spike trains were defined for each unit: all the spikes of the unit, PFS spikes [which were part of (any) PFSs], and Non-PFS spikes (which were not part of any PFS).

Accordingly, for any pair of units (a trigger unit and a reference unit), three cross-correlation functions were computed using one of the above three groups as trigger events, and all the spikes of a reference unit. Figure 6 shows an example of three such functions, computed for one neuronal pair. When taking as trigger events only PFS spikes (dark gray line), there is a clear, central, double-sided peak, reflecting a tendency of the reference cell to increase its firing rate around the time of PFS spikes of the trigger unit. However, when taking all the spikes of the trigger unit as trigger events (black line), the correlation function is almost flat, reflecting independence in the firing between the complete spike trains of the two units. When using as triggers non-PFS spikes (light-gray line), the cross-correlation function consists of a central, double sided, trough, reflecting an anticorrelation in firing between these spike trains. This is a qualitative change in the interunit relationships. This analysis was conducted for 320 pairs made of two different units. Of these, 118 (36.87%) demonstrated such a qualitative change. The change in all cases was that hills, found when taking PFS spikes, disappeared or reversed when taking non-PFS spikes as triggers. This indicates that PFSs appear in periods of coherently elevated firing rates. A result in the same spirit was reported already by Vaadia et al. (1995). This result was found despite the fact that firing rate and coherency were not found as a strong influence on the actual statistical definition of PFSs. The difference is that here not only the PFSs identity is considered but also their number of repetitions. In general, we found that when searching the data for PFS templates (the templates of those PFSs that previously were found as significantly repeating events), the resulted matches often were associated with periods of elevated firing rates. However, there was no one-to-one relationship between the two phenomena: we found many cases where excess of PFSs appeared without rate modulations and periods of rate modulations with no PFSs.


View larger version (34K):
[in this window]
[in a new window]
 
FIG. 6. PFSs associated cross-correlation functions. For each pair of units, 3 correlation functions were computed, based on the selection of trigger spikes: all the spikes of the trigger unit (black histogram); only those spikes that were part of PFSs (dark gray histogram); and only those spikes that did not take part in any PFS (light gray histogram). All the spikes of the reference unit were considered for all 3 histograms. It can be seen that although for group 1 (middle), the correlation histogram is almost flat, it has a pronounced central peak for group 2 (top), and a central trough for group 3 (bottom).


View larger version (33K):
[in this window]
[in a new window]
 
FIG. 7. Behavioral profiles of PFSs. A and B: examples of PFSs with correlation to behavior. Each graph ( bullet ) shows the number of times that a given PFS occurred in each of the 6 behavioral intervals. Intervals were pre-cue period in Go mode (PG) and No-Go mode (PN), delay period after left- or right-directed visual cue in the GO mode (GDL or GDR) and delay period after left- or right-directed visual cue in the No-Go mode (NDL or NDR). Expected counts (· - · and open circle ) were computed using the assumption of uniform rate of the PFS along the recording session. Different expected values in different intervals result from difference in total recording time among these intervals. In A, the PFS <7,9,7; 88,109> tends to accumulate in the pre-cue GO period, whereas in B the PFS <8,12,12; 31,193> tends to accumulate in the delay period, Go mode after visual cue coming from the right. In both cases, the counts were found to significantly deviate from the expected (P < 0.01). Presentation of the behavioral profile and its expectation might be somewhat misleading. A continuous line plot was used, as it helps in comparing the actual counts to the expected values. C and D: 2 examples that demonstrate the differential behavioral profiles of PFSs. Each example contains 2 behavioral profiles of PFSs that share the same unit composition but have a different temporal structure. Each pair of profiles were shown (using a chi 2 test) to be significantly different.

PFSs and behavior

The correlation between PFSs and the animal's behavior is of a major interest. Existence of such correlation is crucial to support the hypothesis that PFSs are related to internal representation and cortical information processing. The study of the correlation between PFSs and the animal's behavior was carried out using two different approaches: the first approach views the PFSs as unitary events in time with a singular composition and temporal structure which carry specific information. The second approach focuses on the spikes that compose the PFSs, leaving aside the detailed temporal structure they create. According to this approach, the PFSs mark spikes from the whole spike train that are more likely to be involved in behaviorally related network activity. This dual approach mainly reflects our uncertainty about the cortical mechanism which underlies PFSs. Thus it is not clear whether each PFS should be treated individually or all of them can be pooled together.

PFSS AS UNITARY EVENTS. The correlation between single PFS and the animal's behavior was quantified by measuring the tendency of each PFS to appear in different behavioral segments. As most of the PFSs did not appear in all the trials performed in the recording session at which the PFSs were defined, we pooled together all the trials. Then for each PFS, we listed the number of its occurrence in each of the six intervals of the DR paradigm (described in METHODS). This list was termed the "behavioral profile" of the PFS and was used to estimate the correlation between single PFSs and behavior (similarly to the use of PSTH). The expected behavioral profile was computed based on the null assumption of constant rate of PFS appearance across intervals. Thus the expected counts for the ith PFS, in the interval m was termed Ei(m), and defined as
<IT>E<SUB>i</SUB></IT>(<IT>m</IT>) = <IT>N<SUB>i</SUB>T<SUB>m</SUB></IT>/<IT>T</IT>;
<IT>m</IT>∈ {PG, PN, GDL, GDR, NDL, NDR}
where Ni is the total number of repetitions of the ith PFS, Tm is the overall recording time during the interval m, and T is the total recording time of all intervals (i.e., T = T1 + ··· T6).


View larger version (29K):
[in this window]
[in a new window]
 
FIG. 8. Clusters of PFSs. A: appearance of PFSs in 3 different single trials that were recorded during the same recording session. Each chart presents the delay period of a single trial. All 3 trials included a 32-s delay period in the No-Go mode. X axis is time from the visual cue. Y axis is a random value set between zero and full scale, which was used to "spread" the plot at this axis. Each dot represents a single PFS. In this case, we considered all the PFSs that were found in the trial. Occurrence time associated with each PFS is the appearance time of its 1st spike. Sequential number of the trial, and the total number of PFS it contained (N) are given for each plot. In all 3 trials, there is a clear tendency of PFSs to cluster and not to appear uniformly along the trial. When summing all the charts plotted for a given recording day (considering only those trials that contained 32 s of a delay period and sorting the trials according the behavioral mode), we got 2 PST histograms (B). Difference between these histograms and the peristimulus time histogram (PSTH), commonly used for evaluating unit response, is that instead of counting the spikes, we counted occurrence of PFSs (represented by the time of the 1st spike of the PFS). Left: was computed for GO mode trials; right: computed for No-Go mode trials. There were many more PFSs in the No-Go mode, and these had a clear tendency to cluster.

The deviation of the behavioral profile from the "expected profile" was measured using a chi 2 test. Figure 7, A and B, provides examples of behavioral profiles (), and their expected profiles (- - -) computed for two single PFSs. In both cases, the behavioral profiles significantly deviated from the expected. Such PFSs were considered to have a correlation with the animal's behavior. On average, 54.76% (weighted2 average 59.62%) of the PFSs (n = 7842) had a behavioral profile significantly deviating from the expected profile (P < 0.01). Note that while computing the expected behavioral profile, the changes in firing rate of the units were ignored. To test the possibility that correlation to behavior is a byproduct of firing rate modulations, we examined the similarity of behavioral profiles computed for different PFSs sharing the same unit composition (and hence differ only in their temporal structure). For example, Fig. 7C shows the behavioral profile of two PFSs. Both PFSs share the same unit composition, <1,11,1>, and differ in their temporal structure (dt1 = 172 ms, dt2 = 214 ms for the first PFS and dt1 = 196 ms, dt2 = 223 ms for the second PFS). The behavioral profiles of these two PFSs were very different: whereas the first PFS tended to occur more often in the No-Go mode after a visual cue from the right, the second PFS tended to appear in the No-Go mode after a visual cue from the left. The difference in the profiles was tested using a chi 2 test and was found to be statistically significant (P < 0.01). This phenomenon was quantified in the following way: at each session, all the different PFSs that had the same unit composition (where the units appeared in the same order) were pooled together. Then for each such group of PFSs, a chi 2 test was carried out to measure the similarity of behavioral profiles among its members. The test was conducted for all the pairwise permutations of PFSs of any given group [i.e., for a group that includes N different PFSs of the same unit composition, a total of N(N - 1)/2 tests were carried out]. The average number of tests per session was 29,376.7. Of this, an average value of 18.96% (weighted average of 24.61%) pairs demonstrated a significant differential behavior, namely contained two PFSs that were different in their behavioral preference. This indicates that not only the unit composition determines the behavioral specificity of the PFSs but also the temporal structure of selected spikes fired by a given set of neurons has a role in determining this specificity.


View larger version (82K):
[in this window]
[in a new window]
 
FIG. 9. Relationships between PFSs and spikes in the box opening paradigm. A: a single trial from the puzzle box opening paradigm. Top 7 traces, represent the activity of 7 different units. Also plotted are the analogue signals encoding the horizontal eye movements of the monkey [electrooculogram (EOG)], and the periods of when the monkey was actively touching the box, trying to open it (Touch). This trace goes upward, whenever the monkey touched the box. This figure contains the complete spike train (CST) of each unit, whereas B contains the reduced spike train (RST) of the same units. Relation of spike-activity to box touching is clearer when examining the RST (B) than when examining the CST (A). C-F: raster displays and PSTH of 2 single units (recorded in a different session than the unit in A and B): unit 8 (C and D) and unit 6 (E and F). Both units were recorded in the puzzle box opening paradigm. Activity is aligned on the time of touching the box. Figure compares between the CST (C and E) and the RST (D and F) of the 2 units. For both units, the degree of modulation (DM) around the trigger event is higher for the RST than for the CST. Note, however, that for unit 6 there is dissociation between the firing rate modulation (the CST in E) and the PFS spikes modulation (the RST in F).

Examination of the specific occurrence time of PFSs along single trials revealed that PFSs tended to appear in clusters. This phenomenon is demonstrated in Fig. 8. Figure 8A presents the onset times of PFSs along the delay period (between the visual cue and the go signal) of three different single trials. All the PFSs that were found during this period were considered. It can be seen that the PFSs were not evenly distributed but appeared in clusters. The clusters were not locked to the stimulus, as in each trial they may appear in different times but their frequency of appearance was similar across trials. The clusters showed correlation with behavior as can be seen in Fig. 8B. In the figure, while clusters can be observed clearly in the No-Go mode, they are much less frequent in the Go mode. The clusters could not be explained by the PSTH of the participating units. It was found that although the PSTH of all the PFSs in the No-Go mode shows very marked peaks (Fig. 8B), the combined-units PSTH, obtained by superimposing all the spike trains (i.e., summing the raw spike trains and then computing the PSTH of the resulted combined spike train), was almost flat. Raising the combined-units PSTH to the third power (as the search for PFSs requests existence of 3 spikes within a given time window and thus the probability of this event is proportional to the multiplication of the 3 corresponding rate functions) produced peaky histograms. However, the peaks in this PSTH did not match the peaks of the PFSs PSTH. The PFSs clustering phenomenon was found in all the recording sessions, of the two behavioral paradigms (i.e., DR and PB).

REDUCED SPIKE TRAIN. To study the possible special significance of those spikes that participated in PFSs, we measured the extent of the stimulus-locked rate modulations of any given unit. This measure was evaluated for the complete spikes train (CST) of the unit, and compared with the value obtained for the corresponding reduced spike train (RST), which contained only those spikes that took part in PFSs. The motivation for this is best explained by the following examples.


View larger version (9K):
[in this window]
[in a new window]
 
FIG. 10. Logarithmic plot of the distribution of the DM improvement index, DMii. Index was computed for 121 cases where responses were observed for both the reduced and the complete spike trains. This index was defined as the ratio between the DM of the RST and the DM of the CST. In 14 cases, the ratio was slightly lower than 1 (with mean of 0.91), whereas for the rest of the cases, this improvement index was >1.

Figure 9, A and B, shows a single trial taken from the PB paradigm. The top traces show the activity of seven single units recorded simultaneously during the performance of this trial. The Touch trace encodes the events of touching and releasing the box while attempting to open it. Whenever the touch line goes up, it means the monkey was touching the box. Very little modulation of the firing rate can be seen in relation to the touch/release act of the monkey. However, when plotting only those spikes that participated in any PFS during this trial, the correlation between this subset of spikes (RST) and the behavior is much more obvious. A similar phenomenon can be seen also when pooling trials together (Fig. 9, C and F). The trigger event for these raster plots is the time when the monkey was touching the box. A PSTH is plotted at the top of each of the raster displays. The CST and RST of two neurons (units 6 and 8) are shown. For unit 8 (Fig. 9, C and D), the ratio between the level of activity after the trigger and its level before the trigger is higher in the RST histogram (Fig. 9D) than in the CST histogram (Fig. 9C). This is an example of enhancement of response-locked rate modulations. For unit 6 (Fig. 9, E and F), the CST shows a moderate reduction in firing rate after the trigger. However, the RST demonstrates an increased level of firing around the trigger. This is an example of dissociation between the characteristics of the unit rate modulation, and the behavior of a certain subset of spikes emitted by this unit: the spikes that are part of PFSs. It can be seen further that the fraction of spikes, which composes the RST, of the total spike train changes in different periods along the trial. This indicates that the modulation of the RST is larger than the firing rate modulation. To test whether the modulation of the RST is not only larger but also more specific in its relation to the behavior, the stimulus-locked degree of modulation (DM) of the RST was compared with the DM of the CST in the following way.


View larger version (75K):
[in this window]
[in a new window]
 
FIG. 11. Single trial activity in the pre-cue Go period of a single cell. Each row in this figure shows activity from the Ready signal to the visual cue. Duration of this period was selected randomly from the range of 3-6 s, and therefore the rows are unequal in length. Top: activity recorded during the 1st trials that were performed after the switching mode signal was given. Bottom: activity during the remaining (non-first) trials. When plotting all the spikes of the neuron (A), a tendency to an elevated activity during the 1st trials can be observed. Difference in activity between 1st trials and the rest of the trials is more pronounced when plotting only those spikes of the unit that were part of any PFS found at that day (B). In C, we plotted those spikes of the unit that were not part of any PFSs (namely the spikes that remain after deleting the spikes of B from A). In this case, no clear difference can be observed between 1st and other trials.

The PSTH around the four different types of visual cues (GO/No-Go, left/right) was plotted for both the CST and the RST. In cases where response existed in both histograms, the onset and offset of the response were defined using the PSTH of the CST. The background level (either at the pre-cue period, or, if the rate modulated in the pre-cue period, after the cue, when the response to the cue returned to background level) was determined as well. The maximal DM for each signal was computed by dividing the maximal signal response value by its mean background level. Then theDM improvement index (DMii) was defined as the ratio ofthe DM computed for the RST and the DM computed forthe CST
<IT>DM</IT><SUB>ii</SUB><IT> = DM</IT><SUB>RST</SUB>/<IT>DM</IT><SUB>CST</SUB>
An improvement index >1 indicates that the DM of the RST is better compared with the DM for the CST. When this value is <1, the opposite is true. The DMii was computed for 121 cases, which fitted the above criteria (Fig. 10). Of these, in 14 cases the index was <1 (although only slightly with an average 0.91), pointing to cases where the RST did worse. For the remaining 107 cases, the index was >1 (the PFSs had a better DM). The average DM ratio was 2.22 (median 1.4) and was found to be significantly different from 1 (t-test, P <<  0.001). Thus for the majority of cases, using only those spikes that are part of PFSs improved the DM. This supports the assumption that the modulation of the RST is behaviorally more specific than the modulation of the CST. The RST also was used to study the correlation between PFSs and behavior on a single trial basis (in contrast to the previous test where data were averaged across trials to yield the PSTH of the units). In this way, it was possible to estimate the variability among trials which are presumably identical in their behavioral context. The motivation for this analysis is presented in Fig. 11. The figure shows a difference between the pre-cue period of trials that occurred after switching from No-Go to Go mode, and the rest of the trials in the Go mode; The firing rate of the unit is different in these two sets of trials. This difference is more pronounced for the RST (Fig. 11B) than for the CST (Fig. 11A). The remaining spikes, left after selecting out all the spikes that were part of PFSs (namely subtracting the RST from the CST), show no clear rate difference between the two types of pre-cue period (Fig. 11C). A way to quantify the improvement in the ability to distinguish between the two sets of trials is using analyses of variance (ANOVAs). The test compared two populations: one containing the firing rates computed for single trials in the pre-cue period of the first trials in the GO mode and the second containing the firing rates computed for single trials in the pre-cue period of other trials in the GO mode. The results of the test done using the RST were compared with the results obtained when using the CST. Two examples of this analysis are shown in Fig. 12. Figure 12, A and B, present the distribution of two rate populations: one was computed during the pre-cue periods of first trials () and the second during the pre-cue period of the remaining trials (- - -). When computing these distributions for the CST (Fig. 12A), the two rate distributions practically overlapped. However, when computing them for the RST (Fig. 12B), there is a clear separation between the two distributions. We designed eight one-way ANOVA tests,3 each testing the ability to distinguish between two related behavioral modes. The behavioral modes and the tests are listed in Table 1. These tests were done for each recording session, and the results varied considerably between sessions, and behavioral intervals. For all (n = 54) except one unit, there was at least one case (average of 3.6 cases) of the eight conducted tests at which the RST yielded an improvement, and for one unit, the RST provided a better separation in all eight tests. This means that the RST improvement of separation was behaviorally specific. This can be seen in Fig. 12, C and D, plotted for the same unit as in Fig. 12, A and B, but now instead of comparing ability to distinguish between first and other trials, the test was the ability to distinguish between pre-cue Go trials () and pre-cue No-Go trials (- - -). Here, in contrast to the previous case (Fig. 12, A and B), the usage of the RST firing rate provides a poorer discrimination. This result further challenges the idea that PFSs are the mere byproduct of firing rate modulation. If indeed the RST represents only some polynomial function of the CST, then the RST is expected to improve any separation already achieved by the rate function (i.e., the CST). However, the behavioral specificity of the improvement counters this argument.


View larger version (24K):
[in this window]
[in a new window]
 
FIG. 12. A comparison between the RST and the CST based on single trial activity. Ability to separate 2 populations of pre-cue intervals was tested. Average firing rate in the pre-cue period of individual trials was computed and taken as a random variable. In A and B, we plotted the 2 distributions of the rates in the 1st trials in each block of trials (First, ) and the other (Others, - - -) trials in the block. A: distribution of rates for the CST; B: distribution of rates for the RST. Difference between the 2 distributions is more pronounced for the RST than for the CST (at which no statistically significant difference could be found). A similar analysis, for the same unit, was done in C and D, but this time the rates were grouped according to the mode: Go () vs. No-Go (- - -). Note that here both for the CST (C) and the RST (D), the difference between the two distributions is small although it is clearer for the CST than for the RST.

 
View this table:
[in this window] [in a new window]
 
TABLE 1. Designing of one-way ANOVA tests between pair of trial groups

    DISCUSSION
Abstract
Introduction
Methods
Results
Discussion
References

This study presents a search for precise spatiotemporal events in cortical spike trains to test the hypothesis that cortical activity is governed by a successive synchronous activation of neuronal groups as suggested by the SFC model. The results indicate that PFSs are found frequently in cortical activity. Further, quantitative analysis of their properties and their relation to the animal behavior suggests that PFSs reflect processes of neural computations. These results support and extend previous reports of precise temporal firing patterns in various regions of the brain (Eckhorn et al. 1988; Gray and Singer 1989; Gray et al. 1989, 1992), and further emphasize the importance of temporal structure of neural activity, as an additional coding dimension beyond the neuronal rate modulations that were taken by many as the only relevant information source. The abundance of PFSs was reflected not only in their high frequency but also in the considerable fraction of spikes (of a given neuron) that were part of PFSs. It was reported previously (Abeles et al. 1993), and confirmed in our results (data not shown), that the fraction of PFS spikes is variable along the trial and could reach <= 100% of the spike train in periods of response.

PFSs and behavior

The results show that PFSs are correlated with the animal's behavior. This correlation was in many cases either stronger than that of the rate modulation (Figs. 9, C and D, and 12) or could not be predicted from it (Fig. 9, E and F). This was true both when treating a PFS as a singular event (Fig. 7) as well as when pooling together the PFSs and focusing on a reduced version of the spike trains that includes only those spikes that were part of PFSs (Fig. 11). It was shown previously that coherency in firing between units (Konig et al. 1993; Roelfsema et al. 1994) and the temporal dynamic of this coherency (Vaadia et al. 1995) can be correlated with perception and behavior even when the firing rate of these units is not correlated. Our results support the above conclusion that interaction between units may convey additional information beyond the one carried by the single neuron firing rate. In addition, our analysis suggests that during a given epoch, some spikes in the train are more informative than other spikes. Similar indications were provided by analysis of bursts in spike trains (Bair et al. 1994; Legendy and Salcman 1985). Here we have further shown that different spikes fired by a single neuron could have different behavioral significance. This could be inferred from the fact that different PFSs with the same unit composition did not always have the same behavioral profile (Fig. 7, C and D). This means that not only the unit identity of the PFSs determines their behavioral profiles but also their fine temporal structure. This offers an additional mechanism for a neural code, which is richer than mere rate code. Accordingly, the spike train of a unit can be subdivided into components, each relating to a different behavioral context. The subdivision is achieved by the temporal structure created by the individual spikes. The information supplied by the temporal structure, in addition to the unit identity, enriches the coding capacity of the neurons. In such a scheme, single cells do not code one unique function but have a dynamic coding capacity.

Are PFSs the outcome of single cell properties?

The results described here do not support the hypothesis that PFSs are the result of single cell properties. We analyzed this possibility by focusing on two kinds of mechanisms: dendritic processes and regular firing. According to the first alternative, PFSs are generated by slow dendritic processes (Barlow 1996). This is unlikely as it is hard to account for a precision of ±1 ms after a pause of 200 ms (the average interunit intra-PFS intervals), while the membrane time constants are in the range of 8-20 ms (Douglas et al. 1988; Shadlen and Newsome 1994), and the effect of events is expected to vanish after two to three time constants. Last it is hard to imagine that a neuron can intrinsically "remember" a given spike for a delay of hundreds of milliseconds and then fire a second, precisely timed, spike while between these two spikes, the very same cell fires other, nonrelated spikes. The fact that most PFSs found in this study were generated by spikes of more than one unit casts further doubts on this explanation.

The second alternative we considered is superpositioning of periods of regular firing across different neurons. This alternative relies on the documented evidence showing that most pyramidal neurons in cortical slices fire regularly (McCormick et al. 1985) and produce intrinsic oscillations (Amitay 1994; Lopes de Silva 1991) in vitro. Accordingly, PFSs may be interpreted as an outcome of temporal overlap of such periods across different cells. This alternative is in conflict with the activity transfer suggested by the SFC model as the source of PFS. Although a direct test is needed to evaluate the relative merits of the two models, the results of this work posed some questions regarding the option of superpositioning.

1) One would expect that such a mechanism will mostly generate 1-U PFSs (i.e., made by spikes of a single unit), whereas this was not the case (Fig. 3).

2) There was no excess of intervals in the range 15-100 ms, which is the corresponding period to the observed in vitro oscillations, which were in the frequency range of 10-60 Hz (Amitay 1994; Llinas et al. 1991; Nunez et al. 1992).

3) About a third of the pairs of units creating PFSs were "unidirectional," i.e., their spikes appeared in PFSs in one order only. Such behavior is not predicted by the superposition model.

4) Even in the case of two spikes in a PFS, emitted by a single unit (i.e., SU pair), other spikes of the same unit were found in between the PFS's spikes (e.g., Fig. 3C).

5) In 1-U PFSs (where all 3 spikes were fired by 1 unit), the regular firing pattern would have been evident, but we found no structure in the firing interval within 1-U PFSs.

We therefore conclude that it is not very likely that the PFSs represent intracellular or single cell mechanisms, and suggest instead that network interactions generate these events.

Are PFSs a byproduct of (coherent) rate modulations?

Special efforts were made here to prevent defining events that reflect rate and coherence modulation as PFSs. Both the searching algorithm, which corrected for the pairwise correlation among the units, and the actual search for PFSs, which was conducted for each segment of the behavior separately, contributed to this effort. Indeed there were cases of cells that fired vigorously and/or exhibited high degree of correlation in firing with other cells but did not create PFSs. When searching again for occurrence of PFSs in the data, the matches that were found often were associated with elevated rates. However, these were not one-to-one relationships (Figs. 7, C and D, 8, and 9, A, B, E, and F). In addition, when looking on a finer temporal scale, the resulting PFSs were shown to be temporally isolated events in the count-expected matrix, namely they were not the outcome of rapid rate modulations (e.g., bursts).

The time locking within PFSs was more precise than their time locking to behavioral events. This means that PFSs are not a byproduct of strongly locked responses across different cells. Furthermore, such a trivial relationship between firing rate and PFSs would yield similar behavioral profiles for PFSs with the same unit composition. However this is not the case. We also found that the RST may provide information which was not found using the CST (Fig. 12, A and B). The additional information was mode specific (Fig. 12, C and D). We therefore conclude (as was previously suggested (Abeles et al. 1993) that PFSs are not a byproduct of rate modulations although these are related phenomena.

Reverberatory SFCs as a source for PFSs

The finding of PFSs in the data lead to the study of the agreement between these events and the SFC model. The original SFC (Abeles 1982b) model was a pure feedforward model. However, the fact that a single neuron usually tended to take part in many different PFSs while most of the PFSs were made of spikes of more than a single unit (Fig. 3), suggested that the mechanism underlying PFSs is a reverberatory SFC network. In such a network, a single cell can participate in different links of the same or different chains (Abeles et al. 1993; Bienenstock 1995; Herman et al. 1995). Activity propagates through groups of cells based on synchronous transmission. These groups can be bound together to create a more complex network characterized by a rich variety of trajectories, which then may lead to a large number of different PFSs, as indeed was found in this study. The reverberatory nature of the circuitry is further demonstrated in our results by the existence of PFS bursts. Accordingly, a set of bursts reflects a sequential activation of modules (sets of links in the chain), whereas each burst represents sustained activity within a module.

The temporal characteristics of the PFSs may shed some light on the extent of the underlying network. Note, however, that in reverberatory circuits the time span of PFSs may be misleading, as it is not clear how many loops of activity traverse through the network while generating a single PFS. The mean interval between two units (200 ms) is a much larger value than the monosynaptic transmission time, which is ~2-10 ms (Bullier and Henry 1979; Ferster and Lindstrom 1983; Martin and Whitteridge 1984). This means that the observed time intervals of PFSs encompass many successive events of interlink transmission of activity volley. It is the limited and dilute sampling as well as the bias against short intervals that was imposed by the searching algorithm that are the sources of the relatively long average interval. These temporal parameters of PFSs raise the question of the ability of SFC to maintain a few-millisecond level of precision while activity traverses through tens of links in view of the expected accumulation of temporal error. The theoretical aspects of this issue were studied (Diesman et al. 1996; Herman et al. 1995), and the ability of neurons to fire precisely timed spikes when the appropriate excitation is provided was demonstrated in vitro (Mainen and Sejnowski 1995). According to these studies, it is safe to assume that the jitter within the spike volley that arrives at each link is short and stable along the SFC (see Gewaltig et al. 1996 for a detailed theoretical discussion). The fact that the intervals between a unit pair were either uniformly distributed (Fig. 4A) or tended to accumulate within a limited number of different intervals (Fig. 4B) can be interpreted as an indication for two extreme modes of reverberating activity: simple mode, characterized by a limited number of activity routes and therefore only a few possible intervals between units that are part of PFSs (peaky interval histogram); and complex mode, characterized by many different routes of activity in the network and thus many different intervals between a pair of units (flat interval histogram). The existence of a complex mode of activity may reflect interaction between different chains, serving for tasks such as binding (Bienenstock 1995).

In any case, for the SFC model PFSs are not a coding entity in a way that requires the brain to measure intervals. PFSs serve as a method to learn about the characteristics of the local network. For the brain itself, it is the synchronous transmission mode in the local network that is of importance. This means that there is no need for any specific readout mechanism tuned to decode the PFSs' intervals into external cues or internal states. The chains themselves are sensitive to sequences of firing configurations and thus serve as their own readout mechanism. The results of this study are therefore in agreement with the SFC model, although they cannot disprove other alternative models. The data, however, call for an emphasis on the temporal structure of the spike trains as an information-carrying dimension of brain activity.

The question of the necessity and usage of a complex delicate temporal coding in the cortex is still open. Many claimed that the rate function of the single neuron serves as a reliable source of information that can account fully for animal performances (Barlow 1994; Newsome et al. 1989; Shadlen and Newsome 1994). On the other hand, others pointed out the necessity and richness of temporal coding (Bienenstock 1995; Hopfield 1996; Von der Malsburg 1985). It should be pointed out that the SFCs provide both elevation of firing rate of single neurons, and specific sequences of firing constellation (as described by MacGregor 1991). Therefore, temporal coding does not necessarily replace rate coding or any other type of coding but may be considered as an additional carrier of information. In such a scheme, some brain systems may readout just the elevation of firing rates, whereas others may exploit the specific information provided by the precise spatiotemporal firing sequences generated by the SFCs.

    ACKNOWLEDGEMENTS

  We thank V. Sharkansky and M. Nakar for technical help.

  This work was supported in part by grants from the Israel-USA Binational Science Foundation, the Basic Research Fund administrated by the Israel Academy of Sciences and Humanities, and the Human Frontiers Science Program.

    APPENDIX

This section presents the algorithm used for detecting PFSs in the data.

Stage 1: computing the counts matrix

The detection of PFSs is based on the threefold correlation matrix of each unit triplet. This triplet correlation among the units was represented by a count matrix Czxy:{c(i,j)}. This matrix contained in its {i,j} cell the number of times in which the sequence (z, x, y, j) occurred, where z, x, and y are the names of the trigger unit and the two reference units, respectively, and i and j are the quantized delays. The delays are given by i = (tx - tz)/B, j =(ty - tz)/B, where tz, tx, and ty are the occurrence times of the spikes of the corresponding units and B is the bin size. The count matrix has a dimension of N × N and therefore counts sequences with a maximal delay T = NB - 1.

Stage 2: computing the expected matrix

The estimation of the expected counts for each bin in the interval matrix was based on the null hypothesis that the only processes governing the accumulation of counts in the matrix are the firing rates and the pairwise correlation in firing between the units. The computation of the expected bin count was done in two steps.

The first step took into account only the contribution of the pairwise correlations between the trigger unit and each of the reference units, namely z right-arrow x and z right-arrow y. These correlation-predictors were designated as
<IT>C<SUB>zx</SUB></IT>(<IT>j</IT>), <IT>C<SUB>zy</SUB></IT>(<IT>i</IT>);
<IT>j</IT>, <IT>i</IT>∈ {1, 2, . . . , <IT>N</IT>}
Each bin (i or j) contains the counts of the total number of times the reference cell fired at that delay after the trigger unit fired a spike. Assuming there were no interactions beyond the pairwise correlation of z to x and z to y, the expected counts in a bin should have been
<IT>e</IT><SUB>(<IT>i,j</IT>)</SUB><IT> = <FR><NU>c<SUB>zy</SUB><UP>(</UP>i<UP>)</UP></NU><DE>n<SUB>t</SUB><UP></UP></DE></FR> × <FR><NU>c<SUB>zx</SUB><UP>(</UP>j<UP>)</UP></NU><DE>n<SUB>t</SUB><UP></UP></DE></FR> × n</IT><SUB><IT>t</IT></SUB>
<IT>e</IT><SUB>(<IT>i,j</IT>)</SUB><IT> = <FR><NU>c<SUB>zy</SUB><UP>(</UP>i<UP>)</UP>c<SUB>zx</SUB><UP>(</UP>j<UP>)</UP></NU><DE>n<SUB>t</SUB><UP></UP></DE></FR></IT>
where Czx and Czy are given in counts and nt is the number of triggers found for z. After obtaining the expected matrix E:{e(i,j)} which is the XY-corrected matrix, the contribution of the correlation between the two reference units (i.e., x left-right-arrow  y) remained and manifested as diagonal bands of counts excess. This was corrected in a second step. The second step of computation was carried out to obtain the global expected matrix E':{e'(i,j)}, in which all the pairwise correlations are corrected for. Here we shall describe this process for one triangular half of the matrix. Generalization to the other half of the matrix is trivial. The bins in the main diagonal were ignored in this analysis, as the counts there could have been the outcome of positive correlation of y right-arrow x as well as x right-arrow y. An ambiguity, which interferes with the statistical estimation of the corresponding expected values).

In the half of the matrix where Tx < Ty, the correlation effect resulted in modulations of the firing rate of y as a function of the delay from a spike fired by x (i.e., x right-arrow y). For the ith diagonal, which contained the following set of ni bins
{<IT>c</IT><SUB>(<IT>j</IT>,<IT>j+i</IT>)</SUB>‖<IT>i</IT>, <IT>j</IT>∈ [0, . . . . , <IT>N</IT>), <IT>i + j < N</IT>}
it was assumed that the average counts should be equal to the average expected
<IT><A><AC>E</AC><AC>¯</AC></A></IT><SUP></SUP>′<SUB><IT>i</IT></SUB>= <IT><A><AC>C</AC><AC>¯</AC></A></IT><SUB><IT>i</IT></SUB> (A1)
Any deviation from this equality was interpreted as a contribution of the correlation Cxy. The contribution was assumed to be constant along each diagonal but presumably different for different diagonals. The average count along the ith diagonal was
<IT><A><AC>C</AC><AC>¯</AC></A><SUB>i</SUB></IT>= <FR><NU><LIM><OP>∑</OP><LL><SUB><IT>j=</IT>0</SUB></LL><UL><SUP><IT>N−i</IT></SUP></UL></LIM><IT>c</IT><SUB>(<IT>j</IT>,<IT>i+j</IT>)</SUB><IT></IT></NU><DE>n<SUB><IT>i</IT></SUB></DE></FR> (A2)
The contribution of the Cxy correlation at each bin along the ith diagonal (j, i + j) was dependent on the number of triggering spikes of unit x, which is in fact the counts at Czx(j), and the average correlation strength per spike of Cxy after a delay i, designated as cri. Note that the units of measure for both the matrices (count and expected) and the correlation vectors are given in counts. Only vector cr has no dimensions. If cri was known, then the value of the expected counts in the {i, j} bin could have been computed based on the assumption that the correlation effects are additive so that the expected value is composed of the XY-contribution [e(j, i+j)], computed in the first stage, and the diagonal correlation effect
<IT>e</IT><SUP></SUP>′<SUB>(<IT>j</IT>,<IT>i+j</IT>)</SUB><IT>= e</IT><SUB>(<IT>j</IT>,<IT>i+j</IT>)</SUB><IT>+ x</IT><SUB><IT>j</IT></SUB><IT>cr</IT><SUB><IT>i</IT></SUB> (A3)
<IT><A><AC>E</AC><AC>¯</AC></A></IT><SUP></SUP>′<SUB><IT>i</IT></SUB>= <FR><NU>1</NU><DE><IT>n<SUB>i</SUB></IT></DE></FR><LIM><OP>∑</OP><LL><SUB><IT>j=</IT>0</SUB></LL><UL><SUP><IT>N−i</IT></SUP></UL></LIM><IT>e</IT><SUP></SUP>′<SUB>(<IT>j</IT>,<IT>i+j</IT>)</SUB>= <FR><NU>1</NU><DE><IT>n<SUB>i</SUB></IT></DE></FR><LIM><OP>∑</OP><LL><SUB><IT>j=</IT>0</SUB></LL><UL><SUP><IT>N−i</IT></SUP></UL></LIM><IT>e</IT><SUB>(<IT>j</IT>,<IT>i+j</IT>)</SUB><IT>+ x</IT><SUB><IT>j</IT></SUB><IT>cr</IT><SUB><IT>i</IT></SUB> (A4)
= <FR><NU>1</NU><DE><IT>n<SUB>i</SUB></IT></DE></FR><LIM><OP>∑</OP><LL><SUB><IT>j=</IT>0</SUB></LL><UL><SUP><IT>N−i</IT></SUP></UL></LIM><IT>e</IT><SUB>(<IT>j</IT>,<IT>i+j</IT>)</SUB>+ <FR><NU>1</NU><DE><IT>n<SUB>i</SUB></IT></DE></FR><LIM><OP>∑</OP><LL><SUB><IT>j=</IT>0</SUB></LL><UL><SUP><IT>N−i</IT></SUP></UL></LIM><IT>x</IT><SUB><IT>j</IT></SUB><IT>cr</IT><SUB><IT>i</IT></SUB><IT>= <A><AC>C</AC><AC>¯</AC></A></IT><SUB><IT>i</IT></SUB> (A5)
The last equality is based on Eq. A1. As cri is still unknown, rearranging these equations yields its value
<IT>cr<SUB>i</SUB></IT>= <LIM><OP>∑</OP><LL><SUB><IT>j=</IT>0</SUB></LL><UL><SUP><IT>N−i</IT></SUP></UL></LIM><IT>c</IT><SUB>(<IT>j</IT>,<IT>i+j</IT>)</SUB><IT>− <FR><NU>e<SUB><UP>(</UP>j<UP>,</UP>i+j<UP>)</UP></SUB><UP></UP></NU><DE><LIM><OP>∑</OP><LL><SUB>j=<UP>0</UP></SUB></LL><UL><SUP>N−i</SUP></UL></LIM>x<SUB>j</SUB><UP></UP></DE></FR></IT> (A6)
As a conservative measure, only positive values of the cri were considered (i.e., those that increase the expected value as computed in the first stage). Negative values were set to zero.

Once the estimation of the correlation factor along the diagonal was obtained (cri), Eq. A3 was used to calculate e'(i,j).

Stage 3: computing the probability matrix

We assumed that the counting process within each bin {i, j} is Poisson. This assumption is based on the fact that at every time step (1 ms in our data), we looked for combined events of a spike fired by unit z, a spike of unit x, fired at a delay of i × B and also a spike of unit y fired at a delay of j × B. This combined event has a very low probability, but we look for it repeatedly great many times (~106 trials for 1,000 s of recording). However, the assumption that the counting process is Poisson does not imply that the spike trains themselves are Poisson. Under the assumption of Poisson counting, the probability of finding c(i,j) counts or more, given an expected value of e'(i,j) is given by
<IT>p</IT><SUB>(<IT>i</IT>,<IT>j</IT>)</SUB><IT>= p</IT>[<IT>c</IT>(<IT>i</IT>, <IT>j</IT>) or more‖<IT>e</IT>′(<IT>i</IT>, <IT>j</IT>)] = <LIM><OP>∑</OP><LL><SUB><IT>k=c</IT>(<IT>i</IT>,<IT>j</IT>)</SUB></LL><UL>∞</UL></LIM><FR><NU><IT>e<SUP>−e<UP>′(</UP>i<UP>,</UP>j<UP>)</UP></SUP>e</IT>′(<IT>i</IT>, <IT>j</IT>)<SUP><IT>k</IT></SUP><IT></IT></NU><DE>k</DE></FR> (A7)

Stage 4: defining PFS in the matrix

Next we had to determine in which of the bins of the matrix the counts were significantly higher than the expected level. Three factors were taken into account: the probability level of the single bin in the matrix, the size of the matrix, and the distribution of significant bins throughout the matrix. The searching algorithm was designed to detect groups of PFSs that are aligned along columns, rows, or diagonals (in both directions) of the matrix. The search for such PFSs layout is based on the observation that often 3 units generated several significant PFSs. These PFSs tended to share common intervals (t1, t2, or t2-t1). Alignment of this kind is expected in a case where PFSs represent fragments of a more complex temporal structure. In such a case, different repetitions of the same complex pattern will yield different triplets, sharing several common intervals.

At each probability level, pi, we first counted for each row the value nrpi, which is the number of bins in the ith row with a probability lower than pi. The expected number of such bins per single row was e = Npi (as there are N bins in a row). The probability per row (pr) for having this number of bins was computed assuming Poisson distribution.

Next, the probability of the rows was tested, starting from a cutoff level of p' < 0.01. This time, for each probability level p'j, the number of rows that are significant at a level less than p'j was counted, where expected chance number was e' = Np'j. If the following inequality was true
<IT>p = P</IT>(<IT>n<SUP>p′j</SUP></IT>‖<IT>e</IT>′) < 0.01 (A8)
all the rows that were significant at a p'j level were revisited, and those bins that were significant at the pi level were taken as PFSs. If Eq. A8 did not hold (i.e., the number of rows that had chance probability below p'j was not significant), we moved to a higher significant level (i.e., p'j+1 < p'j where the ratio is p'j/p'j+1 = radical 2) and checked Eq. A8 again for the new expected value. This process was repeated until we reached the lowest value of p'j (arbitrarily selected to be 8.88e-18). Then we moved to a higher significant level of the individual bins and repeated the whole process forpi+1 = pi/radical 2, recounting the number of bins, their expected numbers, and the probability. After testing the rows, the same process was applied on the columns and the diagonals. For symmetric matrices (where the two reference cells are the same unit), we tested only the rows and half of the diagonals (i.e., from the main diagonal upward).

For a bin to be taken as a PFS, it should have met several criteria.

First, the bin should have contained a minimum of four counts. This threshold was set to avoid selecting PFSs in a sparse matrix where the probability of any event is extremely low. Such PFSs might emerge only due to the discrete nature of our data; second, in cases where the same unit contributed more than one spike to a PFS, (1-U or 2-U PFSs), we set a minimal refractory distance of 9 ms between any two spikes of this unit. This condition was added to reject cases where two spikes of a single PFS are in fact a member of a burst in the single unit activity.

    FOOTNOTES

1   F(t) = Ae{-[t2/(2sigma 2)]}, t is in  (-30,30), sigma  = 21, where A is normalization factor that assures that F(t) has a unit area.

2   Weighted average was obtained by factorizing the contribution of the data of each session to the grand average based on the number of PFSs that were found in that session.

3   Initially, a two-way ANOVA test was used, but it usually yielded significant interactions between the rows and the columns. Therefore a one-way ANOVA test was used instead between each of the four possible pairs of subpopulations.

  Address for reprint requests: Y. Prut, Regional Primate Research Center, Box 357330, University of Washington, Seattle, WA 98195.

  Received 13 March 1997; accepted in final form 10 February 1998.

    REFERENCES
Abstract
Introduction
Methods
Results
Discussion
References

0022-3077/98 $5.00 Copyright ©1998 The American Physiological Society