An important question in neuroscience is whether and how temporal patterns and fluctuations in neuronal spike trains contribute to information processing in the cortex. We have addressed this issue in the memory-related circuits of the prefrontal cortex by analyzing spike trains from a database of 229 neurons recorded in the dorsolateral prefrontal cortex of 4 macaque monkeys during the performance of an oculomotor delayed-response task. For each task epoch, we have estimated their power spectrum together with interspike interval histograms and autocorrelograms. We find that 1) the properties of most (about 60%) neurons approximated the characteristics of a Poisson process. For about 25% of cells, with characteristics typical of interneurons, the power spectrum showed a trough at low frequencies (<20 Hz) and the autocorrelogram a dip near zero time lag. About 15% of neurons had a peak at <20 Hz in the power spectrum, associated with the burstiness of the spike train; 2) a small but significant task dependency of spike-train temporal structure: delay responses to preferred locations were characterized not only by elevated firing, but also by suppressed power at low (<20 Hz) frequencies; and 3) the variability of interspike intervals is typically higher during the mnemonic delay period than during the fixation period, regardless of the remembered cue. The high irregularity of neural persistent activity during the delay period is likely to be a characteristic signature of recurrent prefrontal network dynamics underlying working memory.
The ability to hold and manipulate information in memory for short periods of time has been termed working memory and constitutes a crucial element of higher cognitive functions. Imaging techniques in humans (Courtney et al. 1997; Jonides et al. 1993; Leung et al. 2002; McCarthy et al. 1994) and electrophysiological experiments in monkeys (Funahashi et al. 1989; Fuster and Alexander 1971; Kubota and Niki 1971) have identified the prefrontal cortex as an area critically involved in working memory. A significant proportion of prefrontal neurons show firing rates that are tuned to the position or identity of a sensory stimulus and maintain this tuning persistently across the delay period of a working-memory task, when the stimulus is no longer available to the senses. This tuned persistent activity is regarded as the neural correlate of working memory (Goldman-Rakic 1987, 1995; Wang 2001).
Most previous neurophysiological studies in monkeys have investigated the effect of different stimuli and their maintenance in memory in terms of changes in the mean discharge rate evoked by these events. However, temporal characteristics of the spike train (e.g., oscillations, interneuronal synchronization, or bursting) could, in principle, convey information about task-related events. Oscillatory discharges in particular have received considerable attention over the last decade as the mediator of such functions as attention, figure-ground segregation, memory maintenance, and conscious awareness (Crick 1994; Singer and Gray 1995; Treisman 1996). Oscillations have been reported in the primary visual cortex of the cat and monkey (Eckhorn et al. 1988; Friedman-Hill et al. 2000; Frien et al. 1994; Gray et al. 1989; Livingstone 1996), whereas other studies have contested their prevalence and significance (Ghose and Freeman 1992; Kiper et al. 1996; Young et al. 1992). Fewer studies have examined the existence and possible role of oscillations outside the primary visual cortex, also with conflicting positive (Fries et al. 2001; Kreiter and Singer 1996) and negative (Bair et al. 1994; Cardoso de Oliveira et al. 1997; Tovee and Rolls 1992) results. Notice, however, that positive results are usually observed in population measures like multiunit recordings or local field potentials, whereas single unit activity typically shows little oscillatory behavior. It is conceivable that temporal properties not observed at the neuronal level could emerge from the population collective dynamics (Brunel 2000; Brunel and Wang 2003; Csicsvari et al. 1999).
Recently there are reports suggesting that oscillatory firing may be especially significant in neuronal discharges associated with active memory maintenance in delayed response tasks, compared with baseline spontaneous activity (Pesaran et al. 2002; Tallon-Baudry et al. 2001). Persistent activity during the memory period is believed to be sustained by reverberatory dynamics in a recurrent network (Amit and Brunel 1997; Camperi and Wang 1998; Compte et al. 2000; Durstewitz et al. 2000; Lisman et al. 1998; Wang 1999b, 2001) and it is reasonable to assume that a strongly reverberatory network could promote oscillatory firing under certain conditions, such as when the recurrent excitation is much faster than feedback inhibition (Compte et al. 2000; Tegnér et al. 2002; Wang 1999b). To examine this possibility, we have here applied unbiased analytic methods to investigate the presence of rhythmic activity and its possible correlation with task-related events on a database of neurons recorded in the dorsolateral prefrontal cortex of monkeys as they performed an oculomotor working memory task.
The data used for this study were recorded from 4 male rhesus monkeys (Macaca mulatta), weighing 6-12.5 kg. Detailed descriptions of the experimental techniques as well as the firing rate modulations of recorded neurons in terms of their firing rate modulations were published previously (Chafee and Goldman-Rakic 1998; Constantinidis et al. 2001). Recordings were obtained from the dorsolateral prefrontal cortex that included both the frontal eye fields and area 46 in all 4 monkeys. Two animals (JK and AZ) had additional recordings from areas 7a and 7ip in the posterior parietal cortex. Monkeys were implanted with a scleral eye coil to monitor eye position and a head bolt to stabilize the head during performance of the task (Judge et al. 1980). Surgery and training protocols were in accord with guidelines set by the National Institutes of Health, reviewed and approved by the Yale University Animal Care and Use Committee.
The animals were trained to perform an oculomotor delayed response (ODR) task. They initiated a trial by fixating a central point 0.1-0.2° in size, for 1 s, and continued to maintain fixation as a cue stimulus subtending 0.5-1° flashed for 500 ms at an eccentricity of 13-14°. A delay period lasting 3 s followed the presentation of the cue, at the end of which the fixation point was extinguished and the monkeys were trained to make a saccade to the remembered target location. The cue could appear at one of 8 possible locations, equidistant around the fixation point and randomly interleaved from trial to trial. Ten to 12 correct trials were typically recorded for each location.
Neuronal activity was monitored using varnish-coated tungsten electrodes of 1-4 MΩ impedance, at 1 kHz, as described in Chafee and Goldman-Rakic (1998) and Constantinidis et al. (2001). One or more electrodes were placed in stainless steel guide tubes and advanced into the cortex with the aid of a microdrive. Glass-coated Elgiloy metal microelectrodes (0.6-1.5 MΩ at 1 kHz) were used for a few experiments. Neuronal activity was amplified and band-pass filtered between 400 Hz and 10 kHz. For 2 animals (MAR and COD), the conditioned signal was sampled with a temporal resolution of 30 μs and waveforms were sorted into separate units using a template-matching algorithm (CED, Cambridge, UK). The data acquisition system was configured to classify voltage deflections into separate spikes if the peaks of the deflections were separated by about 1.5 ms, which is shorter than the absolute refractory period of most cortical neurons. Action potentials occurring at shorter intervals could not be correctly resolved. For the other 2 animals (JK and AZ) neurophysiological activity was sampled at 2-ms resolution and action potentials were classified into separate units using a waveform discriminator (Signal Processing Systems, 8701 waveform discriminator, Prospect, South Australia). In practice, action potentials could not be stored as successive samples using this system, so the minimum separation of recorded action potentials was 4 ms.
Neurons were classified as regular spiking (RS) or fast-spiking (FS), putative pyramidal neurons and interneurons, respectively, based on their spike width and baseline firing rate, as described in detail elsewhere (Constantinidis and Goldman-Rakic 2002). Briefly, plotting the baseline firing rate against the width of spike waveform of each unit revealed 2 classes of neurons, one characterized by a longer spike duration and low firing rate and one characterized by shorter spike duration and high firing rate. Approximately 70% of the neurons located in the extreme of the 2 distributions were classified as FS or RS, respectively.
Neurons analyzed in this study were selected from our database of recordings if they satisfied 2 criteria. They displayed activity that was significantly elevated during the delay period relative to the baseline and also exhibited spatial tuning; that is, this activity varied signifi-cantly for the 8 spatial locations tested with the ODR task (ANOVA test, P < 0.05). Some of the analytical techniques for evaluating the spectral properties of the spike train, described below, required a minimum total number of spikes. Given that 10-12 trials were collected per condition in most experiments, this essentially constituted a third requirement of firing rate greater than about 4 spikes/s. A total of 229 prefrontal cortical neurons from the 4 monkeys (169 from COD, 22 from MAR, 21 from JK, and 17 from AZ) that met these criteria were analyzed. For these neurons, we computed autocorrelograms, interspike interval (ISI) histograms, and power spectra for 3 conditions: 1) the last 500 ms of the fixation period; 2) delay period activity for preferred targets, that is, cue locations that elicited a firing rate above the average delay rate across all cues (solid circles in tuning curves of Fig. 4); and 3) nonpreferred cue positions, those eliciting firing rates below the average delay rate (hollow circles in tuning curves of Fig. 4). Additionally, we also applied our analyses to 27 available multiunit recordings (from 2 monkeys) from which 2 or 3 units could be typically isolated.
Power spectrum analysis
In all cases power spectra were computed for each spike train considered, normalized by the expected power of a Poisson spike train at the same firing rate, and then averaged to reduce the variance. Averaging was thus always performed in the frequency domain and on normalized spectra. We computed the power spectrum of a neuronal spike train as described by Jarvis and Mitra (2001). These methods assume the stationarity of the point process analyzed. For this reason, we cut all the spike trains in overlapping 0.5-s windows and treated each window as a different trial of an independent stationary point process. Implicit is thus the assumption that nonstationarities in the data appear only at time scales longer than 0.5 s. A continuous process was built from the spike train by replacing each spike by a delta function and subtracting the average rate, and we subsequently processed these data with Fourier transforms to obtain the power spectra.
The raw estimate of the spectrum, that is, the magnitude-squared Fourier transform of the spike-train data is a biased and inconsistent (i.e., the variance does not decrease with increasing data length) estimate. Windowing the data with a suitable envelope function or taper reduces the bias. This taper function h(t) is nonzero only within the time window of the data and is properly normalized: ∫0T h(t)2dt = 1 (Percival and Walden 1993). To reduce the variance of the estimate (the windowed estimator is distributed as χ2), one can either average over nearby frequencies (lag-window averaging) or use a set of envelope functions that are a complete set of basis functions for time- and frequency-limited functions. One such set of taper functions constitutes the Slepian functions: a set of orthogonal basis functions that maximize the energy in a given time- frequency interval. Because the windows are an orthonormal set, the estimates of the spectra from each window are independent, and hence their average results in a low-variance estimate of the spectrum. We computed Fourier transforms on our spike trains windowed by the first 4 tapers of the Slepian sequences (Percival and Walden 1993). The power spectrum of the spike train estimated in this way was then normalized to the expected power of a Poisson spike train of the same mean rate and spectra from multiple behavioral trials were averaged together. The variance of the normalized spectra was calculated using a jackknife method (Thomson and Chave 1991). The jackknife variance is a powerful method of assessing the variance of samples without any assumptions about the underlying distributions, similar to the bootstrap. The resulting jackknife statistic is t-distributed. For Poisson data, the jackknife variance is a good approximation of the theoretical variance. Specifically, a series of spectrum ensembles were generated by leaving out one trial at a time and averaging over the remaining spectra. The jackknife variance was estimated as the variance across ensembles scaled by the number of degrees of freedom used to compute each member of the ensemble [tapers × (Ntrials - 1)]. The jackknife variance is a robust estimate, even in the presence of nonstationarities.
Coefficient of variation
From the series formed by the ISIs of each spike train we computed 2 measures of the ISI variability: the coefficient of variation (CV) and its local counterpart, CV2. CV is obtained directly from the ISI histogram, by dividing its SD by its mean, . The CV measures how close the spike train is to an ideal Poisson spike train (for which CV = 1), assuming that the data are stationary. Because this assumption was violated by many of our data sets, we used a second measure of ISI variance, which we refer to as CV2. CV2 is computed by comparing each ISI (ISIn) to the following ISI (ISIn+1) to evaluate the degree of variability of ISIs in a local manner (Holt et al. 1996) A Poisson spike train has . For a Poisson spike train with absolute refractory period tr, however, is always smaller than 1. A straightforward calculation departing from the calculations in Holt et al. (1996) yields for a Poisson spike train of rate r and absolute refractory period , which is a nonmonotonic function of rtr. This formula is used to compute tr once and r are known for a given spike train; these values then define the Poisson spike train with absolute refractory period tr that best approximates the data. Furthermore, the probability density of the quantity CV2(n) for adjoining ISIs is given by Thus for a Poisson spike train with absolute refractory period, CV2 depends on the firing rate in a nontrivial fashion.
To assess the degree of burstiness of a given collection of spike trains, we binned all ISIs and defined the following measures: fraction of ISIs <5 ms relative to the fraction of ISIs <5 ms in a Poisson spike train of the same mean firing rate (B1), fraction of ISIs smaller than the mean ISI divided by 10 (B2), or fraction of ISIs smaller than 5 ms (B3). Measures B1 and B2 applied to a Poisson spike train do not change as its rate is varied (in this sense, these measures are rate-independent). However, the measure of burstiness B3 of a Poisson spike train grows with its firing rate. All the results reported here for the autocorrelogram measure A (see following text) were paralleled by B1 and B2. We therefore chose to report only the results for A for the sake of clarity.
We calculated the autocorrelogram of a spike train as the histogram of the time intervals between any 2 spikes, from which the shuffle predictor was subtracted (the shuffle predictor was obtained as the average cross-correlogram of spike trains corresponding to different trials) and the result was normalized to the SD of the shuffle predictor at each time lag (see Aertsen et al. 1989). Thus the autocorrelogram is expressed in units of the shuffle predictor SDs. Our measure of burstiness/refractoriness in a spike train (A) is computed from the autocorrelogram: it is the average height of the central bins that correspond to a time lag <5 ms.
Cell type classification
We used 2 criteria to classify spike trains as “Poisson,” “refractory,” or “bursty”: the departure of the power spectrum in the 5- to 60-Hz band relative to a Poisson spectrum of equal firing rate by more than one SD (bursty: excess power; refractory: suppressed power; Poisson: no structure more than one SD away), and the average height of the central autocorrelogram bins (A). We classified spike trains as bursty if their power spectrum exceeded the Poisson spectrum by more than one SD and A was >1 (central peak more than one SD above the shuffle predictor). Spike trains were classified as Poisson if their power spectrum did not deviate from their Poisson power spectrum by more than one SD and the value of A fell between -1 and 1. Finally, we classified spike trains as refractory if their power spectrum was suppressed relative to the Poisson spectrum by more than one SD and A was smaller than -1 (central trough more than one SD depressed with respect to the shuffle predictor). Notice that some power spectra could satisfy both the “bursty” and the “refractory” criteria and lead to a nondisjoint classification of spike trains. In those cases, we classify spike trains as “bursty” because spectral peaks are often accompanied by power suppression at adjoining frequencies. To classify one neuron according to its spiking characteristics as “Poisson,” “refractory,” or “bursty,” we required that the corresponding classification criterion was satisfied for spike trains recorded in at least one of the 3 task conditions (fixation, delay after preferred targets, and delay after nonpreferred targets). Notice that this classification is nonexclusive: a given cell can belong to various classes. The use of the power spectrum shape to determine the burstiness of a spike train has been discussed by Bair et al. (1994).
Spectral classification techniques: Isomap
To compare spectra in an unbiased way we used the classification scheme Isomap, described by Tenenbaum et al. (2000). As in the case of multidimensional scaling (Mardia et al. 1979), the purpose of Isomap is to provide a visual representation of the pattern of proximities (i.e., similarities or distances) among a set of objects (here, power spectra). However, in contrast to multidimensional scaling (MDS), Isomap does not implicitly assume that distances are taken on a Euclidean manifold but allows for the more general case of geodesics on smooth non-Euclidean manifolds. Thus it is a powerful method to detect the underlying structure of the data without making a priori assumptions about its classifying aspects, and it succeeds in many more cases than classical MDS. In brief, pairwise distances between spectra are determined by using the Itakura-Saito distance where denote the average over frequencies. The Itakura measure is a standard tool for calculating distances between spectra for automatic speech processing. Only segments connecting each spectrum to its 4 closest neighboring spectra are then retained to construct a graph, along which all pairwise distances among spectra are recalculated. The table of distances along this graph is then fed into a classical MDS algorithm (Mardia et al. 1979) to find the best 2-dimensional embedding that preserves the geodesic distance. The output of Isomap consists of the coordinates of each spectrum in this 2-dimensional space. We will focus on the axis that accounts for most of the variance in the population of power spectra (in our graphs, the x-axis) and we will use this coordinate for each spectrum to compare its position in relation to other spectra and thus assess the significance of task dependencies in the shape of the power spectra.
We used various statistical tests to assess the statistical significance of our results. When we used ANOVA tests, we first determined whether the data distributions, or their transformed values, deviated significantly from a Gaussian distribution (Jarque- Berra test, P < 0.01). We also tested the hypothesis of homogeneity of variances of our distributions through Bartlett's test, and we made sure that it could not be rejected at the significance level P < 0.01. Interaction effects of the factors considered in 2-way ANOVA tests were further analyzed by looking at the corrected means (Harwell 1998). When statistical tests showed that the data (or any transformation checked) was not compatible with a Gaussian hypothesis within our signifi-cance level, we used the nonparametric Kruskal- Wallis test (nonparametric version of one-way ANOVA) to assess significance, and we represented the data dispersion as median ± semi-interquartile range. Post hoc analyses were performed by means of a multiple comparison procedure.
All our data analysis was performed in the software package Matlab 6.
A wide range of discharge patterns was observed among the 229 neurons in our data base. Stereotyped patterns of discharge were classified as Poisson, refractory, and bursting as assessed through their autocorrelograms and their power spectra (see Table 1). The most frequent pattern in our database (Poisson discharge) had a virtually flat power spectrum, approaching that which would be expected for an ideal Poisson spike train. An example neuron is shown in Fig. 1A, exhibiting a flat autocorrelogram, an exponentially decaying ISI histogram, and a flat power spectrum. Approximately 89% of neurons (204/229) displayed (see methods for details on neuronal classifi-cation) a flat power spectrum ranging within one SD of the expected power of a Poisson spike train with the same firing rate. Sixty-eight percent of neurons (155/229) had a central peak in the autocorrelogram that did not differ by more than one SD (see methods) from their shuffle predictor (-1 < A < 1). In 63% of neurons (145/229), this stereotyped discharge pattern was present in both measures—the power spectrum and the autocorrelogram. This is the subpopulation of neurons classified as “Poisson” for our subsequent analyses.
Another pattern observed in the sample was characterized by autocorrelograms displaying a central dip and power spectra with significant troughs at low frequencies (Fig. 1B). A significant power spectrum trough in the band 5-60 Hz was identified in 38% of our neurons (87/229). Central dips in the autocorrelogram (A < -1) were clear in 39% of neurons (90/229). Neurons that could be characterized as “refractory” by both measures constituted 26% of all cells in our sample (60/229) and they constitute the class of “refractory” neurons in this study.
Bursting behavior was also present in our database. A cell displaying a representative bursting power spectrum is shown in Fig. 1C. The bursting behavior can be inferred by the sharp central spike of the autocorrelogram, the long tail of the ISI histogram, and by the power spectrum peak in the 20-Hz range. Bursting in this example did not depend on the epoch of the task, despite the marked differences in firing rates during the fixation period and the delay period, or on whether activity was recorded for the preferred or nonpreferred target. Out of 229 neurons 39 (17%) displayed power spectra with the power enhanced significantly over the expected value for a Poisson spike train in the band 5-60 Hz. Autocorrelograms showed a marked central peak (A > 1) in 42% of the neurons (96/229). Overall, 14% of the cells in our sample (32/229, henceforth constituting the class of “bursty” neurons) displayed a “bursting” discharge pattern according to both the power spectrum and the autocorrelogram.
Task-dependent modulation of interspike intervals
We next tested whether the temporal patterns of prefrontal cortical neuron discharge varied systematically across different task conditions. Analysis of the coefficient of variation (either its global CV or local version ; see methods) was higher during the delay than the fixation period but it was not very different between the delay conditions after preferred and nonpreferred cues (Fig. 2). This finding indicated that the coefficient of variation of the ISI was not a function of firing rate alone, as one would expect to observe a similar CV for fixation period and delay period after presentation of a nonpreferred target, which did not differ substantially in rate (13 Hz in fixation and 11 Hz in delay/nonpreferred, average across all cells), and a significantly different CV for delays between preferred and nonpreferred targets, the rates of which were markedly different (18 Hz for delay/preferred and 11 Hz for delay/nonpreferred, average across all cells). This is confirmed in Fig. 2A (bottom right panel), where CV and are shown not to depend on the mean ISI (correlation coefficient -0.03 for CV, and -0.2 for ). The distribution of CV2 for all adjoining ISIs across neurons belonging to each of the classes (Poisson, refractory, and bursty in Fig. 2B, bottom right panels), shows that firing statistics of prefrontal neurons classified as “Poisson” are quite well described by a Poisson spike train with absolute refractory period, neurons in the “refractory” class have some systematic departures from this model but follow it qualitatively, and bursty neurons have a completely different CV2 probability density, not captured by the model (as it should be expected). A second task-dependent effect was associated with the degree of bursting in the spike trains (Fig. 3). Based on the central peak in the autocorrelogram (A), spike trains in the delay period after a nonpreferred target were typically more bursty than either the fixation or the delay period after a preferred target. This also held for neurons classified as Poisson or bursty. However, for cells in the “refractory” class the magnitude of the zero-lag dip increased significantly in the delay period after a preferred cue compared with the other 2 conditions (similar results were obtained from burstiness measures derived from the ISI histogram, not shown). These results indicate that some temporal patterns of spike trains (e.g., burstiness) are independent of firing-rate modulation.
Task-dependent modulations of power spectra
To explore more systematically the modulations in temporal structure of spike trains, we focused on their power spectra. Power spectra can be normalized and error bars can be computed for each point, thus making them a convenient analytical tool for comparisons across task epochs and for determining significant differences. From visual inspection, marked changes in power spectra were not evident in the various epochs of the task. However, weak task modulations were apparent in many cases. Figure 4 shows 3 examples of task-dependent modulations in the power spectrum shape observed in individual neurons of our database. The most frequent task dependency (36/229 neurons, 16%) was a marked deepening of a low-frequency trough in the delay after a preferred cue (Fig. 4A). This was the most salient task-dependent effect observed when we averaged power spectra across the population of neurons (Fig. 5). The effect was consistent across all monkeys and was observed both in the dorsolateral prefrontal cortex, as shown here, and in the posterior parietal cortex (monkeys JK and AZ, data not shown).
In contrast, only a few cells showed power spectra with clear peaks. The neuron in Fig. 4B had a clear peak in the power spectrum that was slightly decreased during the delay period after a preferred target (dependency shared by 6/229 neurons, 3%). Another task-dependency observed (6/229 neurons, 3%) was an enhancement of a peak in the spectrum during the delay, after the preferred target (Fig. 4C).
An averaged spectrogram of 229 neurons was computed for trials in which the cue elicited a delay firing rate above average (preferred cue condition, Fig. 6A) and also for nonpreferred cues (Fig. 6B). We observed very little difference between the spectral properties of neuronal firing before and after cue presentation except for a suppression of power at low frequencies (<20 Hz) after the presentation of preferred cues (Fig. 6A). The spectrogram shows the time course of this suppression of power: it began during cue presentation and persisted throughout the delay period. For trials with nonpreferred cues only a slight power increase at low frequencies was observed (<20 Hz; see Fig. 6B).
Classifying power spectra shapes
To make a more unbiased comparison of the spectra across task epochs, we used the classification scheme Isomap (Tenenbaum et al. 2000). The classification obtained this way is unbiased in the sense that the major differences between the spectra are decided by the algorithm rather than a predetermined spectral feature that we are interested in (for details see methods). The resulting 2-dimensional map is shown in Fig. 7. Note that the principal axis, which accounts for most of the variance (x-axis), separates spectra that show a peak (large positive values of x) from those that show a pronounced trough (large negative values of x). This is shown more clearly in Fig. 8A, where the position of a spectrum on the x-axis of the Isomap was found to correlate remarkably well with the average power in the 5- to 20-Hz band, whereas it correlated only weakly with the average power in the band 20-60 Hz. Even this weak correlation with the 20- to 60-Hz power range can be accounted for by the correlation of the latter with the power in the 5- to 20-Hz band (Fig. 8A, right panel). This result confirms our previous observation that the most significant structures in the power spectra occur at low frequencies (<20 Hz) and it is in that frequency range that the spectra distinguish themselves from one another.
We proceeded to test whether there was any correlation between the principal coordinate of each spectrum on the Isomap and the degree of bursting of the corresponding spike trains. Bair and colleagues (1994) reported that bursty cells typically show a peak in the power spectrum. If peaks in the spectra of our database are primarily caused by bursty firing rather than by nonbursty periodic firing we would expect to see a pronounced correlation between the principal coordinate of the spectra and the bursting of the spike trains. This was indeed so, as shown in Fig. 8B (right panel).
We also investigated how the electrophysiologically identified classes of RS and FS neurons related to the classifications of the spike trains that we have used here: Poisson, refractory, and bursty neurons. This classification showed a significant effect, in that FS neurons are overrepresented among refractory neurons (Fig. 9). This result is similar to the observation from behaving rats that FS putative interneurons tend to be characterized by autocorrelations with a large dip at zero time lag (Csicsvari et al. 1999). Also, an overrepresentation of RS neurons was observed among neurons with very high bursting (for A >5 there are 18 neurons classified as RS and 1 as FS, the same effect is seen in the power spectra; see histogram tails in Fig. 9), in agreement with previous results (Constantinidis and Goldman-Rakic 2002).
We then examined whether the shape of power spectra, as classified by the Isomap algorithm, varied depending on the task conditions (Fig. 10). To assess that, we considered for each cell the Isomap location of its spectra corresponding to fixation, delay after preferred targets, and delay after nonpreferred targets (two examples in Fig. 10, left top panels). For each set of 3 power spectra associated with a given neuron, we calculated their center of mass along the principal axis (crosses in Fig. 10) and recorded the signed distance “r” between the location of each of the spectra and their center of mass, along the Isomap principal axis. A positive value of r indicates that a spectrum corresponding to a particular task condition was shifted toward the “peak”-side of the Isomap with respect to the other task conditions. The 2 examples shown in Fig. 10 indicate the r value associated with the spectra shown in Fig. 4. The 2 spectra displayed opposite sign r values, as the Isomap classification distinguishes peaks and troughs. When the r values for spectra in each of the behavioral conditions were computed, the resulting empirical distributions were tested against the hypothesis of the equality of means (Fig. 10). There was a clear tendency, consistent across monkeys, of power increase at low frequencies (positive r) in the delay period after nonpreferred targets, whereas trough deepening (negative r) was typical in the delay period after preferred targets. Spectra for spike trains recorded during fixation were typically close to the center of mass of the 3 conditions.
Finally, we analyzed a smaller sample of multiunit recordings along the main lines of our single unit analysis. The results are illustrated in Fig. 11. Because of the small sample, some effects observed in single units did not reach significance (task dependency of burstiness, Fig. 11B; and task dependency of spectral shape, Fig. 11C right panel). Other results, however, were indeed replicated: variability of the ISIs (CV and ) also increased significantly in the delay (Fig. 11A), and the main feature of the power spectra was their structure at frequencies below 20 Hz (Fig. 11C, left panel. As for single units in Fig. 8A, the high-frequency index did not correlate more with the principal coordinate than with the low-frequency index, not shown), which correlated very well with the degree of burstiness of the spike trains (Fig. 11C, middle panel). No overall increase of spectral peaks or troughs relative to single-unit recordings was revealed in the multiunit records.
We performed an analysis of the temporal structure of prefrontal neuron spike trains during the fixation and delay periods of a working memory task. Our main findings are threefold. First, the variability of ISIs for a neuron is typically higher during the mnemonic delay period than during the fixation period, regardless of whether the cue is preferred or nonpreferred. Second, the spectral properties of neuronal discharges in most neurons approximated the characteristics of a Poisson process. In a minority of cells that deviated significantly from the expected Poisson pattern, some showed a trough at low frequencies of the power spectrum and a pronounced dip near zero time lag of the autocorrelation function; others displayed a peak (at <20 Hz) in the power spectrum associated with the burstiness of spike discharges. Third, we observed a small but significant task dependency in the spike-train temporal structure. Neuronal firing during the delay period for preferred locations was characterized not only by elevated firing rate, but also by spectral power of low (<20 Hz) frequencies suppressed below the expected power spectrum of Poisson spike train. This power suppression was most prominent in FS neurons (putative interneurons).
General characteristics of spike train temporal properties
We identified 2 types of deviations from the expected power spectrum of a spike train with Poisson characteristics: spectral peaks and troughs. Both of these were almost always observed in frequencies below 20 Hz. Significant spectral peaks were strongly correlated with increased burst firing. Previous studies have shown that firing of action potentials in bursts produces power spectral peaks, even when bursts themselves are not repeated periodically, with a constant interburst interval (Bair et al. 1994). Our current results are in agreement with these findings. Burst discharges are characteristic of a subclass of pyramidal neurons in the cortex (Connors and Gutnick 1990; Gray and McCormick 1996). Recent work suggests that bursts may subserve special types of coding and computations in sensory neurons (Doiron et al. 2003; Kepecs et al. 2002). What may be the significance of burst firing to the working memory circuits in the cortex? One possibility stems from the evidence that excitatory synapses in the prefrontal cortex exhibit short-term facilitation (Hempel et al. 2000). Hence bursts may represent a reliable neural signal through unreliable yet facilitating synapses (Lisman 1997; Wang 1999a), thereby contributing to the synaptic reverberation in the prefrontal circuits.
Our data show that a great proportion of neurons showing a significant dip in the autocorrelation function and power suppression at low frequencies are classified from the shape of the recorded spikes as FS neurons (putative interneurons). This is consistent with a physiological study of behaving rats reporting that FS interneurons tend to show a pronounced dip at zero time lag of their autocorrelations (Csicsivari et al. 1999). Such a dip in the autocorrelation and power suppression at low frequencies are usually related to an effective refractoriness in the neuronal spike train. Therefore interneurons apparently engage more importantly in this task-induced enhancement of the refractoriness of their output. The cellular and synaptic mechanisms, and functional implications, of this effect are currently unknown. Inhibition has been identified by computational models (Brunel and Wang 2001; Camperi and Wang 1998; Compte et al. 2000; Tegnér et al. 2002) as well as neurophysiological studies (Funahashi et al. 1989) as a mechanism responsible for the suppression of firing with respect to baseline for those pyramidal neurons nonselective to the cue during the delay. The contribution of temporal dynamics of interneuronal spike discharges to the generation of such suppressive inhibition is an important issue to be elucidated in the future.
Absence of oscillatory discharges in the autocorrelation and power spectrum
Since the original descriptions of oscillatory firing in the primary visual cortex of the anesthetized cat (Eckhorn et al. 1988; Gray et al. 1989), there have been several studies in the cat and monkey that failed to detect oscillatory firing, or reported that such oscillations were largely unrelated to the properties of visual stimuli (Bair et al. 1994; Ghose and Freeman 1992; Tovee and Rolls 1992; Young et al. 1992). Other investigators detected oscillatory firing but in a smaller proportion of the neurons tested (Bringuier et al. 1997; Cardoso de Oliveira et al. 1997; Nowak et al. 1995; Schwarz and Bolz 1991). Our present results add to this debate because truly oscillatory discharges were very rare in the spike trains of prefrontal cortical neurons. Some obvious factors that could be responsible for discrepancy among the various studies, such as the use of anesthesia, stimulus-induced periodicities, analysis criteria, and species differences, have been addressed by subsequent studies (Eckhorn et al. 1993; Friedman-Hill et al. 2000; Frien et al. 1994; Fries et al. 2001; Gray and Viana Di Prisco 1997).
The present study was designed to test the central hypothesis that elevated firing rate during the delay period following a target in the neuron's memory field was associated with oscillatory discharges, as suggested by Lisman and Idiart (1995) and Tallon-Baudry et al. (2001). Our present results suggest that, at least in individual neurons of the prefrontal cortex, this is clearly not the case. Power enhancement over the expected value for a Poisson spike train occurred only at frequencies below 20 Hz, and they were not related systematically to the behavioral task in any of 4 monkeys. This is in contrast to similar analysis performed on single-neuron recordings and local field potentials in posterior parietal cortex (Pesaran et al. 2002). Pesaran and colleagues detected significant task-related power increase in the gamma band in single parietal neurons. The spectral techniques employed there are analogous to the ones that we apply here. The only differences in the technical aspect of the analysis are their inclusion of non-memory-selective (according to the firing rate) neurons, the database size (40 neurons in their study vs. 229 neurons in the present study), and our unbiased classification through the Isomap method. It seems unlikely that these reasons can entirely account for such contrasting results. This raises the possibility that neurons in the prefrontal cortex have less propensity to display oscillations than those in the parietal cortex. A number of computational studies have shown that oscillations are detrimental to working memory function sustained by local cortical reverberation (Compte et al. 2000; Tegnér et al. 2002; Wang 1999b). For example, suppression of oscillatory activity results in a greater tolerance to intervening stimuli (Compte et al. 2000). There is evidence that working memory-related activity in the prefrontal cortex is more resistant to distracting stimuli than in other cortices, like the inferotemporal cortex and posterior parietal cortex (Constantinidis and Steinmetz 1996; Miller et al. 1996; Sakai et al. 2002). Thus the absence of oscillatory activity in prefrontal single units during a delayed response task could reflect the stability of persistent activity supported by local cortical reverberation, and, at the same time, resistant to distracting stimuli.
Studies with recordings from motor cortex reported maximal rhythmic firing during periods preceding the onset of a cue and reduction of spectral power after the instruction for a particular movement, even though the firing rate modulation of single neurons follows the opposite pattern (Baker et al. 2001; Murthy and Fetz 1996; Sanes and Donoghue 1993). These results led to the suggestion that oscillatory firing in the motor cortex may be associated with a higher level of attention or arousal, in anticipation of a sensory event, rather than with the encoding of stimulus information per se. However, we do not find oscillatory firing in prefrontal neurons recorded during the oculomotor delayed response task, even though this task requires substantial attention allocation throughout each trial. Furthermore, the power suppression detected during delay for preferred cues cannot possibly be associated with a release of attention because that is the period of the task most demanding in terms of attentiveness for the animal.
The possibility remains that oscillations at the population level are not easily detectable at the level of firing patterns of individual neurons (Brunel 2000; Brunel and Wang 2003; Csicsvari et al. 1999). Indeed, the data we analyzed here correspond to single-neuron activities. Moreover, we noted that the results held for a smaller sample of multiunit recordings incorporating 2 or 3 units each (Fig. 11). On the other hand, the identification of task-related oscillatory activity has found its strongest evidence in measures of global activity of a large ensemble of neurons, such as local field potential (Pesaran et al. 2001) or intracranial EEG recordings in the occipital and parietal lobes (Raghavachari et al. 2001; Tallon-Baudry et al. 2001) during delayed response tasks. Oscillations may be more readily apparent in multiunit recordings, given that such recordings increase the probability that one of the units will display oscillations, and even more so in local field potentials. Unit isolation therefore may be partly responsible for the large differences in percentage of neurons displaying oscillatory firing in different studies. Therefore it would be interesting to undertake local field potential recordings from the prefrontal cortex during activation of working memory. The functional significance of rhythmic firing, detectable only in the activity of large numbers of neurons, as exemplified by local field potentials, but not in the spiking output of individual neurons, remains unclear at the present time.
Power spectral changes related to memory function
Neurons in our database showed little change in the spectral characteristics of their spike trains between the fixation and the delay periods, irrespective of whether the cue was preferred or nonpreferred. Assuming that the spectral features of a spike train are primarily affected by the degree of attentiveness required for a task, the fact that there were no major differences in modulation of power spectra between the fixation and delay periods in the present database could be explained on the basis that the fixation requirement in an ODR task also requires substantial sustained attention.
Nevertheless, small but significant spectral effects of the task could be detected. We used the Isomap method to quantify the power spectrum's structure. We found that the first principal coordinate of the Isomap correlates well with the average power in the 5- to 20-Hz band, and distinguishes the presence of a trough or peak in that frequency range. Using this measure, and consistently for all 4 monkeys, we found that the spectral signature of active memory maintenance in single prefrontal neurons is a suppression of power at low frequencies (<20 Hz) if the cue presented is preferred, and a slight enhancement of power if the cue is nonpreferred (Fig. 10). A similar phenomenon has been observed in posterior parietal cortex during this same task (Fig. 5 in Pesaran et al. 2002), although there a power increase in the gamma band was also very prominent. It is thus possible that refractoriness is a shared signature for single neurons during working memory across cortical areas, whereas gamma band oscillations are a result of the local network connectivity and much more robust in parietal cortex than in prefrontal cortex. The comparison of our results with those of Pesaran et al. (2002) leads to the hypothesis that, in spike trains recorded during working memory, task-dependent power suppression below 20 Hz and gamma-band power increases are attributed to distinct mechanisms. As discussed in the following text, we also provide suggestive evidence that interneurons might be involved in the generation of the power suppression effect below 20 Hz.
Temporally irregular nature of persistent neural activity
We found that the ISI's coefficient of variation (CV) of a prefrontal neuron is consistently high (see also Shinomoto et al. 1999). Previous studies have reported that the ISI CV of visual cortical neurons in response to a fixed sensory stimulus is close to 1, the value for a Poisson process (Sofkty and Koch 1993). When the stimulus is naturalistic and time-varying (videos), the CV of V1 and inferotemporal neurons is as high as 1.8-1.9 (Baddeley et al. 1997). In the latter case, the large CV may be directly related to the slow temporal changes of the neural firing rates. Here we report that during the delay period of an oculomotor task, when there is no external stimulus and the monkey's eye position is maintained, the CV of mnemonic activity of prefrontal neurons is above 1, on average, and larger than during the fixation period, regardless of whether the cue is preferred or nonpreferred. Single-neuron modeling shows that, generally, the CV is larger when the neuron fires at a lower rate (Liu and Wang 2000; Sofkty and Koch 1993; Troyer and Miller 1997) or is more bursty (Wilbur and Rinzel 1983). Therefore when the cue is nonpreferred, a larger CV of the delay activity may be explained by a lower firing rate (compared with the fixation period) and perhaps by a slightly higher propensity of burstiness across the neural population (as observed in Fig. 3). On the other hand, when the cue is preferred, the neuronal firing rate is higher (up to more than 60 Hz), and burstiness is smaller, yet the CV is still larger. Whether this high CV of prefrontal neurons is a result of slow temporal changes of delay period activity (see, for example, Chafee and Goldman-Rakic 1998; Romo et al. 1999) or reflects the balanced dynamics between synaptic excitation and inhibition (Shadlen and Newsome 1993) in the working memory circuit, remains to be explored in the future. Also, the incidence of task-dependent burstiness could be significant because it has been shown that fluctuations of neural activity and synaptic inputs at low or high frequencies may have differential impacts on the reverberatory network dynamics (Nowak et al. 1997).
The observed irregularity of persistent activity has potentially important implications. Because high variability of neural discharges is likely to originate from stochastic synaptic bombardments (Shadlen and Newsome 1993; Shu et al. 2003; van Vreeswijk and Sompolinsky 1996), our result is consistent with a network mechanism for generating mnemonic delay activity in the prefrontal cortex. By contrast, it is unclear how the high ISI variability could be explained if persistent activity is a single-neuron phenomenon (Egorov et al. 2002). Moreover, the scenario of balanced synaptic excitation and inhibition (Shadlen and Newsome 1993; van Vreeswijk and Sompolinsky 1996), proposed to account for a high CV close to 1, does not seem to be compatible with the premise of multistability between a resting state and active memory states (Renart 2000; Renart et al. 2003). Therefore the fact that delay activity of prefrontal neurons in behaving monkeys indeed shows a high variability, as reported here, highlights a major inadequacy of the existing working memory models. Resolution of this issue will help to elucidate the microcircuit organization and dynamics of working memory cortical networks.
This work was supported by the National Science Foundation, the National Institute of Mental Health, the Alfred P. Sloan Foundation, and the Swartz Foundation (to X.-J. Wang); MH-38546 Grant to P. S. Goldman-Rakic; a McDonnell-Pew Program in Cognitive Neuroscience Award to C. Constantinidis; and support from the Wennergren Foundation, the Swedish Medical Research Council, and the Fernstroms Foundation to J. Tegnér.
We thank M. N. Franowicz who participated in some of the recordings, Y.-H. Liu for contribution to the initial phase of data analysis, and B. J. Pesaran for helpful comments on the manuscript.
The costs of publication of this article were defrayed in part by the payment of page charges. The article must therefore be hereby marked ”advertisement” in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.
- Copyright © 2003 by the American Physiological Society