|Home | About | Journals | Submit | Contact Us | Français|
Cortical areas differ in their patterns of connectivity, cellular composition, and functional architecture. Spike trains, on the other hand, are commonly assumed to follow similarly irregular dynamics across neocortex. We examined spike-time statistics in four parietal areas using a method that accounts for non-stationarities in firing rate. We found that, whereas neurons in visual areas fire irregularly, many cells in association and motor-like parietal regions show increasingly regular spike trains by comparison. Regularity was evident both in the shape of interspike interval distributions and in spike-count variability across trials. Thus Poisson-like randomness is not a universal feature of neocortex. Rather, many parietal cells have reduced trial-to-trial variability in spike counts that could provide for more reliable firing-rate signals. These results suggest that spiking dynamics may play different roles in different cortical areas, and should not be assumed to arise from fundamentally irreducible noise sources.
It is widely held that neocortical neurons in vivo fire action potentials with irregular interspike intervals. Extensive effort has been aimed toward understanding the functional significance of this observation. In one view, irregularity reflects noise: cells fire at certain rates but the exact timing of spikes is random (Shadlen and Newsome, 1998) and is often modeled as a Poisson process with a brief refractory period between action potentials (Gerstein and Mandelbrot, 1964; Kara et al., 2000). Consistent with a Poisson model, variance-to-mean ratios of spike counts in visual cortex often equal or just exceed unity (Britten et al., 1993; Buracas et al., 1998; McAdams and Maunsell, 1999; Tolhurst et al., 1983). It has been posited that Poisson-like spike trains are the fundamental unit of cortical communication (Geisler et al., 2005; Ma et al., 2006; Mazurek and Shadlen, 2002; Shadlen and Newsome, 1998).
The Poisson interpretation of irregularity is inconsistent with some recent observations. For example, in rat auditory cortex many cells fire a single spike in response to pure tone stimuli, and the probability of observing a spike is better modeled as a binomial rather than a Poisson process (DeWeese et al., 2003). Other studies report that neurons in visual cortex can exhibit sub-Poisson variability in anesthetized, paralyzed preparations (Kara et al., 2000), when responding to natural images (Amarasingham et al., 2006), or after accounting for the effects of small eye-movements (Gur and Snodderly, 2006). These data argue that the Poisson model may provide an overestimate of noise levels in the system and that spike-time statistics may vary across the neocortical mantle.
In this light, we decided to systematically examine spike-time statistics in parietal cortex. An important axis upon which spiking neurons vary is the regularity of interspike intervals. At one end of such a spectrum are irregular-spiking cells: neurons that show Poisson-like, or super-Poisson, variability (Figure 1, top left). At the other end are highly regular, clock-like neurons, such as those found in the brain stem or cerebellum, in which there is very little difference between neighboring interspike intervals (Figure 1, bottom right). In this study, we wished to determine where cortical cells lie along this continuum. In particular, we wanted to answer two questions. First, do parietal neurons show characteristic, first-order, interspike interval statistics independent of firing rate or task? Individual cortical cells are known to have distinctive anatomy, connectivity, and intrinsic electrophysiological properties; perhaps single parietal neurons also show idiosyncratic interspike-interval statistics in vivo. Second, if individual cells do show distinctive regularity or irregularity in their interspike intervals, do the spiking statistics of single cells differ across parietal regions?
Whereas regular spiking is common in the peripheral nervous system (Werner and Mountcastle, 1965) and subcortically (Calvin and Stevens, 1968; Pfeiffer and Kiang, 1965), neocortical neurons are thought to fire in a Poisson-like fashion in awake animals. This view, however, arises mainly from recordings in the visual system. Fewer studies have examined spike timing outside of sensory areas (Averbeck and Lee, 2003; Churchland et al., 2006; Compte et al., 2003; Joelving et al., 2007; Lee et al., 1998). Working with awake, behaving monkeys, we recorded from single neurons in four parietal cortical areas (Figure 2A): the middle temporal (MT) and medial superior temporal (MST) visual areas, the lateral intraparietal (LIP) association area, and area 5, a motor-like sub-division of the parietal lobe in which cells activate with contralateral arm movements. We found that cells in MT and MST typically show irregular interspike intervals and their spike-count variance-to-mean ratios are near 1. In contrast, many neurons in LIP, and an even larger fraction of cells in area 5, fire with more regular interspike intervals, with variance-to-mean ratios clearly below 1. We show that cells with increased regularity from Poisson tend to have large, broad extracellular waveforms, suggesting that at least some of these neurons are pyramidal cells (Bartho et al., 2004; McCormick et al., 1985; Nowak et al., 2003) that may function to send a lower-noise spike-rate signal subcortically or to other cortical sites.
We recorded neural activity from two monkeys as they performed the following task. The animals viewed a dot approaching a bar on a computer screen and were trained to make visually guided arm movements in one of two modes (Figure 2B) (Experimental Procedures and Supplemental Experimental Procedures). In the proactive mode, the monkeys pressed a lever to cause the dot to reverse its motion direction before it collided with the bar. In the reactive mode, the dot reversed its motion direction under computer control and the monkeys pressed a lever in response to this abrupt cue. This task allowed us to characterize brain processes underlying the initiation of proactive versus reactive arm movements, described previously (Maimon and Assad, 2006a; Maimon and Assad, 2006b).
In addition, this experiment was advantageous for assessing spiking statistics of cortical neurons. First, we recorded from parietal regions across the visual-to-motor spectrum, using the same methodology in each region. Second, the monkeys performed a large number of identical trials allowing us to robustly estimate statistical distributions. Third, as the task was demanding, spike count variability on correct trials was unlikely to arise from trial-to-trial differences in overall states of arousal/vigilance.
Figure 2C (top) shows activity from two visually responsive cells, aligned to the time of the dot’s motion reversal. The stimulus was oriented such that these direction-selective neurons were stimulated in their preferred direction after the motion reversal. Whereas the two cells fired at similar rates, they exhibited very different patterns of spiking. The MST cell fired irregularly, whereas the LIP cell had more evenly spaced spike times. To the right we show interspike interval (ISI) distributions from a 400-ms time window, starting 100 ms after the turn. The MST neuron had intervals that approximate an exponential distribution, a signature of a Poisson process. The LIP neuron had a more rounded, compact ISI distribution, showing that for a similar mean interval length the LIP cell had fewer long and short intervals about that mean, a phenomenon we call “pulsed firing”. Data from two motor-like cells in area 5 (Figure 2C bottom), aligned to the time of the hand movement, make the same point. Some cortical neurons fire with Poisson-like irregularity, but others fire in a more regular fashion than Poisson.
At high firing rates, where the mean interspike interval approaches the time scale of the refractory period, even a Poisson-like process appears more regular (Gabbiani and Koch, 1998). This is because refractoriness precludes very short intervals and a high firing rate precludes very long intervals, which together give rise to ISI distributions that are necessarily more compact and non-exponential-like. A related issue is that some cortical neurons fire action potentials in bursts of two or more spikes in rapid succession, during which firing rates can be extremely high. Bursting can thereby introduce periodic structure to spike trains (Bair et al., 1994). Neurons we recorded in parietal cortex, however, showed regular spiking at both low rates and high rates, and they did not burst (e.g., Figure 2C, cells in red). We therefore define a “pulsing neuron” as a cell whose interspike intervals are regular independent of firing rate, i.e., not due to simple refractoriness, and whose regularity does not result from a tendency to burst. We use the term “pulsing neuron” by analogy to the pulse of the heart, which changes rate over time, but maintains a fundamental regularity. As will become clear, pulsing neurons have interspike intervals that are more regularly distributed than those derived from a Poisson-like process; however, such cells in cortex do not approach the level of oscillatory precision observed in many sub-cortical neurons, let alone man-made clocks. To emphasize this latter point, we avoid the label “periodic neuron” or “oscillatory neuron”.
As a measure of regularity, we quantitatively characterized ISI distributions. One complication with analyzing ISI plots was the fact that each cell’s spike rate varied considerably over time. This nonstationarity meant that we could not simply plot a single ISI histogram for each cell, as this histogram would have undoubtedly resembled a sum of more basic distributions, one for each rate at which the neuron fired.
To address this issue, we adopted a method for plotting multiple ISI histograms per neuron, in which each ISI histogram consists of interspike intervals associated with a similar mean firing rate (Softky and Koch, 1993) (Experimental Procedures). In this analysis we first calculated a firing-rate histogram across all trials (peri-event time histogram, PETH; 50 ms bins) aligned on a trial-specific event, such as the time of the hand movement or the dot’s motion reversal on the screen (Figure 3A, gray curve). Then, for each trial we uniformly multiplied this overall PETH by a scalar, defined as the number of spikes on the currently analyzed trial divided by the average number of spikes across all trials (Softky and Koch, 1993). For example, if we observed 50 spikes on a given trial, whereas the mean number of spikes was 30, we multiplied each bin of the overall PETH by 1.67 (Figure 3A, red curve), so as to account for the increased response gain on that specific trial. If we observed 17 spikes, then we multiplied each bin by 0.57 (Figure 3A, black curve), so as to account for the lower gain. Next, we associated each interspike interval with the trial-specific firing rate at the central time point of that interval (e.g., dotted red line in Figure 3A associates an interspike interval of 43 ms with a firing rate of 61 spikes/s). Finally, we plotted an ISI histogram for intervals associated with 0–5 Hz firing, a separate histogram for 5–10 Hz firing, and so on, in 5 Hz steps, up to the maximum rate at which the cell fired.
Having generated a set of ISI histograms for each cell, we fit these histograms to simple functions so as to describe their shape with a small number of parameters. Visual inspection of the data, as well as past work (Kuffler et al., 1957; Seal et al., 1983), argued that ISI plots would be well modeled as gamma distributions. Gamma distributions are defined by a scale parameter (θ) and a shape parameter (k), with k providing a useful measure of ISI-regularity (Figure 1 and Figure S1). That is, when k = 1, gamma distributions are mathematically equivalent to exponential distributions, suggestive of a Poisson process. When k > 1, the peaks of the gamma distributions are offset from zero, suggestive of increased regularity from Poisson. When k < 1, gamma distributions have an excess of very short and very long intervals, suggestive of increased irregularity from Poisson. Thus, the shape parameter provides a convenient axis for separating pulsed spike trains, with relatively high k, from more irregular spike trains, with low k.
Before applying this ISI analysis to our complete data set, there were some important issues to consider. First, we only fit ISI histograms with >150 interspike intervals, as histograms with fewer data points were under sampled (e.g., lower-left histogram in Figure 3B). Second, since we wanted to detect neurons whose ISI regularity was not a direct result of refractoriness when firing at a high rate, we focused on the same relatively low spike rates, 10 to 40 Hz, for all neurons. More than 90% of cells provided sufficient interspike intervals in the 10 to 40 Hz range (see below), allowing us to effectively compare the shape of ISI histograms across neurons after normalizing for firing rate, without excluding many cells from analysis. Third, since we were interested in examining neurons whose ISI-regularity was not due to simple bursting or refractoriness, we deliberately ignored very short interspike intervals (≤ 8 ms) when fitting gamma distributions to ISI plots (Experimental Procedures; we return to these short intervals later).
ISI histograms, across the board, were very well described as gamma distributions (Experimental Procedures). In Figure 3B, we show rate-parsed ISI histograms from the four example cells in Figure 2C. For each cell, we found that the shape parameter typically increased slightly with firing rate, but otherwise remained consistent across the fitted histograms (e.g., Figure 3B k1–k6, quantified in more detail below). The area 5 neuron in the bottom row of Figure 3B shows one of the largest variations in shape parameter observed in our data set, ~ ±35% about the mean, yet even this cell’s histograms did not change character from a higher-order Gamma to a more exponential-like appearance across different rates. For a given neuron, we averaged the shape parameter, k, across each of the fitted histograms, generating a single measure of ISI-regularity for each cell. The MST neuron had a mean k of 0.9, consistent with irregular interspike intervals. The LIP neuron had a mean k of 5.1, evidence for pulsed firing. The irregular-firing neuron from area 5 had a mean k of 0.6, well below 1, suggestive of super-Poisson irregularity. The other area 5 neuron had a mean k of 7.6, evidence for pulsed firing.
To verify that our gamma model was not generating a biased picture, we also calculated the coefficient of variation (CV) of each ISI distribution, which provided a parameter-free method of describing the ISI plots. CV is a measure of dispersion, defined as the standard deviation divided by the mean. Exponential distributions have a CV of 1, whereas distributions derived from more regular interspike intervals have CV values below 1. As with the gamma fits, we only calculated CV for well sampled histograms (>150 interspike intervals) and we averaged CV across histograms. In Figure 3B, the two irregular cells had mean CV values of 1.1 and 1.2, whereas the two pulsing cells had CV values of 0.52 and 0.39.
One important consideration is that pulsed spike trains are expected to have lower spike count variance than irregular spike trains; in the extreme, a grandfather clock will always tick ten times in ten seconds, whereas a Geiger counter measuring radioactive decays at a rate of 1 per s will show variability about this value. Thus, to the extent that cortical neurons communicate with spike-rate signals, pulsing neurons should provide a more reliable form of this signal. For each bin of the PETH we calculated the mean and variance of the number of spikes recorded across trials and plotted a single point on a variance-vs.-mean log-log plot (Figure 3B right). The two irregular-firing cells had points that fell along the diagonal, or slightly above it, indicating that spike-count means and variances were roughly equal, as predicted from a Poisson process (Figure 3B first and third rows). In contrast, the two pulsing cells had points that fell below the diagonal, indicating lower variance than predicted from a Poisson process (Figure 3B second and fourth rows). To quantify these data, we calculated the mean residual (vertical distance) of points from the diagonal. The two irregular cells had mean residuals that were slightly positive (0.13, 0.25), whereas the two pulsing cells had mean residuals that were substantially negative (−.75, −1.2). In summary, for each neuron we calculated three metrics to describe its spike-train statistics. Two metrics, k and CV, captured the shape of ISI distributions, whereas the third metric, mean residual, captured spike-count variability across trials. Note that with mean residual we did not parse the intervals in any way, and we examined all firing rates for each cell, not just 10–40 Hz.
Our complete data set consisted of 539 parietal neurons: 102 in MT/MST, 203 in LIP, and 234 in area 5. We analyzed 94% of MT/MST cells (96 of 102), 95% of LIP cells (193 of 203), and 94% of area 5 cells (219 of 234); other neurons did not fire sufficiently between 10–40 Hz to provide at least one ISI histogram with >150 interspike intervals. Overall, spiking in MT/MST was irregular (median k = 1.23), with increasing regularity in LIP (median k = 1.95) and area 5 (median k = 2.70) (Figure 4A left column). All three shape-parameter distributions were significantly different from one another (rank sum test; P < 10−7), demonstrating that first-order spike-time statistics vary across parietal regions. Cells with high k we call “pulsing neurons” and those with low k “irregular neurons”, but cells fell along a continuum with no obvious boundary separating the two extremes. Nevertheless, if one defines an arbitrary threshold at k > 3, one finds that 41% of area 5 neurons exceeded this threshold and had spike trains that appeared definitively regular to the naked eye (Supplementary Information). 21% of LIP cells, but only 4% of MT/MST cells (4 neurons), had k > 3.
Results were comparable with our alternative measure of ISI regularity, mean CV (Figure 4A middle column). As CV is a ratio (s.d./mean), we plotted CV on a logarithmic scale in which values near zero indicated Poisson-like spiking and values below zero indicated pulsed spiking. We observed a median log(CV) of 0.01 for MT/MST, −0.04 for LIP, and −0.13 for area 5. All three distributions were significantly different from one another (rank sum test; P < 10−4). For reference, we also show CV distributions on a linear scale (Figure S2). CV and k were highly correlated (Figure S3) arguing that each provided a similar, reliable measure of ISI regularity.
For each cell we generated a variance-vs.-mean log plot and calculated the mean residual from the equivalence line, averaged across all bins. Neurons with low log(CV) tended to have negative mean residuals (Figure 4A, right column). Correlation coefficients were 0.90 for MT/MST, 0.84 for LIP, and 0.85 for area 5, all highly significant (F-test, P < 10−15). In other words, pulsing cells have lower spike-count variability than irregular cells, demonstrating that pulsing neurons carry a more reliable spike-rate signal than previously thought to exist in the neocortex. This finding was not critically dependent on the bin width used for PETHs or on how we aligned and processed the data (Figure S4).
Our analysis assumes that spike-train irregularity or regularity is a defining feature of a neuron and not something that changes dramatically with firing rate. Stated formally, our working model is that for a given neuron, interspike intervals are drawn from a gamma distribution with a consistent shape parameter (k), but a scale parameter (θ) that changes with spike rate (see Figure S1). If this model is valid, k measured for a given neuron should be the same independent of firing rate, which was largely the case (Figure 4B). Pulsing neurons did not turn irregular at certain spike rates or vice versa.
There was, however, a weak tendency for k to increase with higher rates, and for log(CV) to decrease, likely because of the refractory period (Figure 4B). Thus, if area 5 neurons fired faster than MT/MST neurons this could have contributed to the increased regularity in area 5. To control for firing-rate differences across areas we only analyzed ISI plots in the 10–40 Hz regime. Nevertheless, MT/MST cells might have fired closer to 10 Hz and area 5 cells closer to 40 Hz, for example.
To address this issue, we examined spike rates across areas (Figure 5). We plotted a vertical line from the minimum to the maximum firing rate exhibited by each cell, and positioned this line along the x-axis based on the cell’s log(CV) measure (Figure 5 left). There was substantial overlap in the range of firing rates among individual cells across all cortical areas. We also plotted the maximum, minimum, mean, and standard deviation of spike rates as a function of log(CV) (Figure 5 right). Mean spike rates in MT/MST (mean = 16.8 Hz) and area 5 (mean = 16.6 Hz) did not differ significantly (rank sum test, P = 0.93), and thus cannot explain the observed differences in ISI regularity. Standard deviations of firing rates in MT/MST (mean = 12.3 Hz) were significantly higher than those in LIP (mean = 10.3 Hz) and area 5 (mean = 7.6) (rank sum test, P < 10−4). This larger range of spike rates for MT/MST neurons was likely because we optimally positioned the visual stimulus along their preferred-null directional axis (Experimental Procedures). In contrast, the animals always performed the same lever press, even though many area 5 cells are selective for the direction of arm movements (Kalaska et al., 1983). Only 1 of 12 plots relating firing rate parameters to log(CV) (Figure 5 right) showed a significant correlation (P < 0.05, Bonferroni correction for 12 comparisons); area 5 neurons with low log(CV) tended to have slightly elevated minimum spike rates (r = −0.29; P = 1.8 × 10−5). Overall, these data emphasize the extensive overlap in firing rates between irregular and pulsing cells and across cortical regions. Differences in ISI regularity are thus unlikely to be due to systematic differences in firing rate.
Another concern is that increased ISI regularity could result from bursting (Bair et al., 1994). For example, we observed ISI histograms that had a peak at short interspike intervals, reflective of bursting, followed by a trough and then a second, lower peak (Figure 6A). Since we excluded the first 8 ms in our gamma fits, these cells would have had elevated k parameters simply because of the post-burst trough, which likely reflects a post-burst refractory period (Bair et al., 1994). To examine the role of bursting in pulsed firing, we measured the proportion of spikes in the 0–8 ms bins of each ISI histogram and compared this value with the proportion of spikes predicted based on the gamma fit, which excluded these bins. For example, the histogram in Figure 6A shows data from a cell with a large difference between the actual and predicted number of spikes in the 0–8 ms bins. In contrast to this example, however, we generally observed a good match between predicted (red symbols) and observed (black symbols) values, with perhaps slightly more intervals observed than predicted for rates close to 40 Hz (Figure 6B). The most bursting was in MT/MST, less in LIP, least in area 5. That is, 57% of ISI histograms in MT/MST had more 0–8 ms intervals than predicted by a Poisson process, i.e., an exponential distribution (Figure 6B gray lines), whereas only 27% of histograms in LIP and 17% in area 5 showed this trend. Bursting did not contribute to the pulsed spike trains characterized here. In particular, area 5 cells showed the least bursting and most spike-time regularity.
For two reasons, it is possible we overestimated spike-time irregularity and spike-count variability in MT/MST. First, we typically aligned MT/MST spike trains to the time of the dot’s motion reversal, which had 3–5° of trial-to-trial jitter (Experimental Procedures). Since receptive fields in MT/MST are not homogeneous, this trial-to-trial difference in the visual stimulus could have increased the variance of spike counts across trials. Second, microsaccades (small, rapid eye movements) within the fixation window would have produced irregular movements of the visual stimulus on the retina that could have led to irregular changes in the firing rate of MT/MST cells, particularly since our stimuli had sharp contrast edges. Area 5 cells have very weak, if any, visual responses and would not have been pushed towards irregularity in the same manner.
To address trial-to-trial jitter in the location of the dot’s motion reversal, we analyzed spike trains up to, but not including, the dot’s earliest motion reversal among all trials (Supplemental Experimental Procedures). In this analysis, the visual stimulus on the monitor was identical across trials, yet we observed the same trends of irregularity and variability across the three brain regions (Figure 7A). To address the role of small eye movements, we extracted the time of microsaccades from our data set and found, serendipitously, that one of our monkeys performed microsaccades at twice the rate of the other animal (Figure S5). Despite this difference in saccade rate, the two animals showed very similar spiking statistics in MT/MST (as well as in LIP and area 5) suggesting that microsaccades did not have a strong influence on our results. We also analyzed data from the baseline fixation period (Supplemental Experimental Procedures), during which there was no visual stimulus on the screen and thus microsaccades are less likely to have introduced visually driven irregularity. In addition, we recorded from the same population of cells as the animals performed a memory-saccade task and in the memory-delay epoch of this task there was no visual stimulus in the receptive field (Supplemental Experimental Procedures). In both the fixation period of the main task, as well as the memory period of the memory-saccade task, we observed the expected pattern of irregularity and variability across parietal areas (Figure 7B and 7C). Note that we recorded fewer trials in the memory-saccade task and less cells provided sufficient interspike intervals for analysis. Overall, MT/MST neurons exhibited similarly irregular firing across a broad range of conditions, arguing that spike-time irregularity is a general feature of MT/MST and not an artifact of data alignment or eye movements. Indeed, in all areas, spike-time regularity or irregularity remained a consistent feature of neurons across different tasks and time epochs within each task.
To minimize the influence of the classical refractory period on our results, we focused on relatively long interspike intervals, associated with 10–40 Hz spiking, in our ISI analyses. However, the relative refractory period could still play a role at these rates. One possibility is that pulsing neurons were in fact Poisson-like, but had very long relative refractory periods, lasting tens of milliseconds, which would result in fewer short intervals and thus higher k. If so, one would expect the tails of ISI histograms, after the effects of refractoriness have completely subsided, to resemble exponential distributions, not gamma distributions, but this was not the case (Figure S6). Rather, interspike intervals from pulsing neurons seemed genuinely well captured as draws from gamma distributions with k > 1. One could still consider pulsing cells to have extended refractory periods, in that these cells have reduced probability of very short intervals. However, these neurons also have reduced probability of very long intervals, meaning that an extended refractory period is not the only difference between pulsed and irregular spiking.
Do irregular and pulsing neurons correspond to distinct cellular subtypes? Extracellular spike waveforms can provide some insight into this question. The width of action-potential waveforms in cortex fall into two distinct categories, narrow and broad (Bartho et al., 2004; McCormick et al., 1985; Mitchell et al., 2007; Mountcastle et al., 1969; Rao et al., 1999), and some have argued that narrow spikes correspond mainly to inhibitory interneurons, whereas broad spikes correspond to a wider set of neurons including pyramidal cells (McCormick et al., 1985; Mitchell et al., 2007; Nowak et al., 2003). We measured the time interval between trough and peak of our spike waveforms (Experimental Procedures). 5–12% of neurons in each area had waveforms in which the peak occurred before the trough; these cells were not analyzed further. For the remaining cells, we found a bimodal distribution of spike widths in MT/MST, LIP and area 5 (Figure 8A), much like reported previously in area V4 (Mitchell et al., 2007). We observed a weak, but significant correlation between spike width and log(CV) both in area 5 (r = −0.33, P = 1.6 × 10−6) and LIP (r = −0.29, P = 1.1 × 10−4), but not in MT/MST (r = −0.21, P = 0.06) (Figure 8B).
As a measure of spike amplitude, we also calculated a “signal-to-noise” ratio for each spike waveform (Figure 8C inset) (Experimental Procedures). Mean signal-to-noise ratios were 3.3 for MT/MST, 4.2 for LIP, and 4.8 for area 5 (Figure 8C). All distributions were significantly different from one another (rank sum test, P < 10−3). Cells with large signal-to-noise ratios tended toward pulsed firing. That is, there was a weak, but significant correlation between signal-to-noise and log(CV) in area 5 (r = −0.22, P = 0.0012), LIP (r = −0.21, P = 0.004), and MT/MST (r = −0.45, P = 5.2 × 10−6) (Figure 8D). These correlations remain significant if we exclude one outlier in MT/MST with a signal-to-noise ratio of 12 and one outlier in area 5 with a signal-to-noise ratio of 27 (data point not shown). Big, broad waveforms suggests that many pulsing cells are large pyramidal neurons rather than smaller interneurons.
An influential view holds that spike timing in the neocortex is fundamentally irregular. In contrast to this model, we have identified parietal neurons that fire with increased ISI regularity. These “pulsing neurons” fire in a notably regular fashion across different spike rates (Figure 4), different time periods within a task (Figure 7), and different tasks (Figure 7), suggesting that regularity is defining feature of these cells and not a property that changes with task demands. Pulsed spiking serves at least one important function: increasing the reliability of spike-rate signals (Figure 4A right column). Pulsing neurons are most common in association and motor-like regions of the parietal lobe, less common in visual areas (Figure 4A, left and center column). Although past reports have not emphasized regular spiking in neocortex of awake primates, ISI histograms in rat archicortex (hippocampus) also appear better described by gamma rather than Poisson models, after accounting for the extra parameter in the gamma model (Barbieri et al., 2001). Additionally, a common method for estimating neuronal latencies relies on modeling ISI histograms as gamma distributions (Seal and Commenges, 1985), a method first developed using spike trains from area 5 (Seal et al., 1983).
Since pulsed spiking reduces variance in spike-rate signals, why do so few MT/MST cells fire in this fashion? The answer is not clear; however, if MT/MST cells must respond to visual inputs with millisecond precision, pulsed firing is problematic as there would be interspike periods, sometimes lasting many milliseconds, with extremely low probability of spiking. Motor computations, in comparison, may operate on a longer timescale – as these signals ultimately get filtered by muscles – and thus may be more suitable for pulsing. Consistent with this notion, spike trains in the spinal cord are extremely regular (Calvin and Stevens, 1968). Our results argue that irregular spiking in visual areas does not reflect irreducible noise inherent to all cortical computation (Shadlen and Newsome, 1998), since this noise is reduced in other cortical areas. Rather, widespread irregularity in the visual system is more likely to reflect noise specific to visual circuits or, alternatively, the need for rapid temporal signaling.
Aside from increasing the fidelity of spike-count signals, what functional role is served by pulsed spiking? We examined whether pulsing cells were more likely to participate in the proactive timing of arm movements – the cognitive process assayed by our behavioral task – but found no supporting evidence (Figure S7A). Additionally, there was no detectable relationship between the width of spike waveforms and cognitive firing-rate modulations in our task (Figure S7B). To constrain functional interpretations, it will be important to determine the anatomical inputs and outputs of pulsing neurons. LIP and area 5 send and receive signals from sub-cortical motor circuits, which commonly show regular spiking (Clower et al., 2001; Lynch et al., 1985; Sasaki et al., 1976). Furthermore, it has been shown that motor and cognitive regions of the frontal lobe exhibit closed-loop connectivity with the basal ganglia and cerebellum, regions with extensive rhythmicity (Kelly and Strick, 2003). Perhaps LIP and area 5 are more influenced by this closed-loop architecture (Clower et al., 2005) than MT/MST.
A general concern with extracellular recording techniques is the likely bias toward recording from large cells with high spike rates (Olshausen and Field, 2005). In this study, we made an effort to record from well-isolated units, large or small, that had even a minimal response to the task. Nevertheless, it remains a possibility that spiking statistics were identical across parietal areas, but pulsing cells in MT/MST, for example, were too small for us to consistently isolate. Arguing against this idea, however, pulsing cells in LIP and area 5 tended to have large-amplitude waveforms (Figure 8) and high minimum firing rates (Figure 5). It is unlikely that we missed such neurons in our MT/MST recordings. More probable is a true areal difference in the statistics of firing.
Past studies (Bair & Koch, 1996; Buracas et al. 1998; Masse & Cook 2008) have reported that MT neurons can generate spike trains with high temporal precision across trials, which seems at odds with the notion that spike trains in area MT are fundamentally irregular. In these studies, however, temporal precision in MT responses was a direct result of dynamic visual stimuli, which drove the cells rapidly from a low to high firing rate at multiple time points in a trial. Spike trains were consistent across trials precisely because the stimulus-induced rate changes occurred at consistent moments across trials. However, the underlying spike-generating process in MT neurons was still reported to be well modeled, to first-order, as a fundamentally irregular Poisson-like process (Bair & Koch, 1996; Buracas et al. 1998). For example, variance-to-mean ratios of spike counts in one study were near 1 (Buracas et al. 1998). Indeed, irregular neurons may be more capable than pulsing cells of phase locking to dynamic inputs because irregular cells do not experience a prolonged refractory state between spikes. More generally, one should carefully distinguish stimulus-induced regularity, as found in MT, from a regular underlying mechanism of firing, which appears to be a stimulus-independent feature of cells in LIP and area 5 (Figure 7).
Along similar lines, a notable set of papers have examined the responses of somatosensory cortical neurons to vibrotactile stimuli (Mountcastle et al., 1969; Recanzone et al., 1992; Romo et al., 1998; Salinas et al., 2000). In these studies, it was shown that many neurons in the primary somatosensory cortex (S1) respond with an oscillating firing rate to an oscillating tactile input in the range of 5 to 100 Hz (Mountcastle et al., 1969). In the secondary somatosensory cortex (S2), however, many fewer neurons respond periodically to the oscillating input (Salinas et al., 2000). A progression from more periodic firing in S1 to less periodic firing in S2 might seem at odds with our finding of increased pulsed firing in higher compared to lower parietal regions. Again, however, the periodicity observed in these studies was directly due to periodicity present within the incoming stimulus. The results are to be expected if S1 neurons have intrinsically irregular interspike intervals, with short refractory periods, which allow them to better phase lock to periodic inputs compared to S2 neurons. Indeed, a Poisson model for the underlying spike generating process provided a good fit for information-transfer rates observed in S1 cells (Salinas et al., 2000). An intriguing, albeit speculative, hypothesis is that higher somatosensory regions may show more pulsed firing, and less periodic following of inputs, than lower somatosensory areas.
At this point, one might interpret our results to mean that spike time regularity increases consistently from primary sensory to higher cortical regions. Whereas this trend appears to hold in the parietal areas studied, it is unlikely to be uniformly true. For example, it has been observed that many neurons in dorsolateral prefrontal cortex fire quite irregularly, especially during memory periods of delayed response tasks (Compte et al., 2003). Perhaps a more reasonable interpretation of our data is that neurons involved in oculomotor and especially skeletomotor computations show disproportionate levels of pulsed firing. This interpretation raises the testable prediction of more pulsed firing in primary motor cortex and frontal eye fields compared to MT or MST. Indeed, if information is communicated in Morse-code like patterns of interspike intervals in neocortex, a debated hypothesis, irregular spiking would be associated with higher entropy levels and thus allow cells to transmit more information per unit time. In this view, one might imagine greater irregularity in higher centers of cognitive management, like the dorsolateral prefrontal cortex, which could possibly require a higher capacity for information transfer compared to regions involved in motor control.
Increased ISI regularity could result from membrane properties in single neurons, circuit properties of networks of neurons, or a combination of both. A clearer understanding of the mechanisms for pulsed spiking may require work in more reduced preparations. Specifically, to disambiguate cellular versus network mechanisms, one could test for regular firing in dissociated neurons (e.g., Raman and Bean, 1997), where synaptic input does not contribute to the activity.
A related question is whether pulsing cells show higher order structure in their spike trains, or, whether, as our first-order model suggests, interspike intervals are best captured as independent draws from gamma distributions. Future work will be needed to definitively sort this issue out. Nonetheless, it should be mentioned that global electrical oscillations in the brain provide a potential source for long range order. The relationship between such oscillations and spike timing is an important area of active research (Buzsaki, 2006). We could not directly address how pulsed spiking relates to global brain rhythms because we did not record multiple neurons simultaneously, nor did we measure local-field or scalp potentials. Nevertheless, one way to approach this question in our data set is to examine single-cell autocorrelograms, which plot the probability of observing a spike at a given time delay after another spike. Unlike ISI distributions, autocorrelograms examine all intervals, not just nearest-neighbor intervals. If pulsing cells were locked to a circuit-level oscillation, one might expect distinct oscillatory peaks in autocorrelograms, whereas if interspike intervals were generated cell-autonomously then one would predict autocorrelograms to appear more flat, especially at long time delays (Figure S8). Autocorrelograms of highly regular area 5 spike trains did not exhibit strong oscillatory peaks, suggesting that these interspike intervals may have been divorced from global rhythms. This conclusion is tentative, however, given that very few parietal cells fired with sufficient oscillatory precision so as to expect detectable peaks in autocorrelograms, even if these cells were locked to network oscillations (Figure S8). Additionally, nonstationary firing rates corrupt periodic structure in autocorrelograms and it is not clear how to normalize autocorrelograms for firing rate given that these plots require an extended time window for analysis.
No single analysis provides a complete characterization of neuronal spiking. We therefore also show unprocessed rasters from our entire data set (Supplemental Data). Across cells there was an excellent match between the regularity of spike trains as assessed by eye and as assessed by k or CV of rate-parsed ISI histograms. Examining the entirety of our spike-train data makes clear the pronounced differences in spike-time statistics across parietal neurons and places emphasis on deciphering the origin of this impressive variability.
Two male rhesus monkeys (Macaca mulatta, ~10kg) were surgically implanted with a head post, scleral search coil, and recording chamber. All surgical and experimental procedures followed Harvard Medical School and NIH guidelines. The chambers were dorsally positioned at stereotactic coordinates P3, L10 in the hemisphere contralateral to the hand used for pressing a lever. Extracellular recordings were made from single neurons using tungsten microelectrodes and a guide-tube/grid system (see below). We recorded from four distinct cortical regions: the middle temporal (MT) area, the medial superior temporal (MST) area, the lateral intraparietal area (LIP), and area 5. Recording sites were localized based on physiological responses and structural MRI images, as described previously(Maimon and Assad, 2006a).
We recorded spike times with 1-ms resolution using a dual time-and-amplitude window discriminator (Bak Electronics, Mount Airy, MD). In parallel, the complete voltage trace was streamed to disk and used offline for investigating spike waveforms (see extracellular waveforms below). We recorded neurons one by one, examining action potentials on an oscilloscope as they appeared. Only unambiguous single units, as determined by eye, were analyzed. With 1-ms spike-time resolution, even a perfectly isolated neuron could fire in adjacent time bins; thus we did not employ a quantitative requirement for a refractory period. However, for reference, we show extracellular waveforms from our entire population of cells (Supplemental Data). To minimize bias, we recorded from all neurons, large or small in amplitude, that responded even minimally to our behavioral task.
We recorded a total of 539 cells: 309 in monkey B and 230 in monkey R. From monkey B, we recorded 54 cells in MT/MST, 127 cells in LIP, and 128 cells in area 5. From monkey R, we recorded 48 cells in MT/MST, 76 cells in LIP, and 106 cells in area 5. In early recordings we used tungsten microelectrodes with a 250 μm shaft diameter (FHC, Bowdoin, ME). In later experiments we used 75 μm electrodes (FHC, Bowdoin, ME), as these yielded more stable isolations. With monkey B, 54/54 MT/MST cells, 5/127 LIP cells, and 25/128 area 5 cells were recorded with 250 μm electrodes; the rest were recorded with 75 μm electrodes. With monkey R, all cells were recorded with 75 μm electrodes. Our results in MT/MST were very similar between monkeys (e.g., Figure S5), even though electrodes of different widths were used in each animal.
The electrode signal was amplified, bandpass filtered at 100–6000 Hz (3 db/down) with a notch filter at 60 Hz (Bak Electronics, Mount Airy, MD), and streamed to disk at a 20 kHz sampling rate (CED, Cambridge, England). Offline, we extracted the 80 voltage samples (4 ms) surrounding each spike time and averaged together the first 100 spike waveforms to generate one mean waveform for each cell.
Spike times were recorded with 1 ms of uncertainty by the window discriminator, thus we needed a method of aligning waveforms before averaging. Instead of aligning to the time of peak or trough, which could vary by up to 50 μs (1 sample period), we used an algorithm that took into account all waveform sample points. For each cell, we slid the first waveform past the second waveform and found the time shift that maximized the correlation between the two curves. The pair of waveforms were aligned to this time shift and averaged together. We then slid the third waveform past the mean waveform from the first two spikes and, again, found the time point that maximized the correlation. The third waveform and mean spike were aligned to this time shift and a weighted average was taken, in which the mean waveform was weighted by 2/3 (as it included two spikes) and the third waveform weighted by 1/3. We continued this process for 100 spikes, thus generating the final mean waveform for each cell.
We interpolated the mean waveform with a cubic spline to achieve 10 μs temporal resolution. We then extracted two metrics for each cell, spike width and signal-to-noise. Spike width was defined as the time difference between the maximum and minimum voltages of the mean waveform. Spike signal-to-noise was defined as: spike amplitude/noise amplitude. We calculated spike amplitude by subtracting the minimum voltage of the mean waveform from its maximum voltage. For estimating noise, we calculated the mean voltage in a 0.2 ms baseline period, starting 1 ms before the trough of each spike. We then measured the standard deviation of these measurements across all 100 spikes for a given neuron. Noise amplitude was defined as four times this standard deviation.
The behavioral tasks have been described in detail previously(Maimon and Assad, 2006a; Maimon and Assad, 2006b) and a complete description is provided in Supplemental Experimental Procedures. In brief, monkeys pressed a lever in order to cause a moving dot to reverse its direction of motion, or, in alternate blocks of trials, the monkey pressed the lever in reaction to a computer-generated reversal of the dot’s trajectory. Overall, per neuron we analyzed an average of 199 trials in area 5 (min 78; max 554), 217 in LIP (min 120; max 458), and 227 in MT/MST (min 119; max 495). Only correct trials were analyzed. In generating perievent time histograms (PETHs), spike trains from MT/MST and LIP were aligned to the time of the dot’s motion reversal and spike trains from area 5 were aligned to the time of the arm movement. Our results were not affected by these differences in how we pre-processed the data (Figure S4). Multiple PETHs were generated for each cell, one for each dot turn or arm movement in each task. Overall, 98.7% of PETHs had > 35 trials and thus provided a robust estimate of firing rate. Specifically, 1.3% of PETHs had < 35 trials, 61.5% of PETHs had 35–45 trials, 0.9% of PETHs had 45–55 trials, 21.8% of PETHs had 55–65 trials, and 14.6% of PETHs had > 65 trials.
To evaluate the first-order statistics of cortical spike trains, we examined rate-parsed interspike interval (ISI) histograms (Figure 3). These were generated by applying an algorithm developed by (Softky and Koch, 1993) to each PETH and related set of rasters. For each trial we uniformly multiplied the overall PETH by a scalar, defined as the number of spikes on the currently analyzed trial divided by the average number of spikes on a trial. In other words, we assumed that a cell’s firing rate evolved with the same time course on all trials, but the overall gain changed on any given presentation. Next, we associated each interspike interval with the trial-specific firing rate at the central time point of that interval. Finally, we plotted an ISI histogram for intervals associated with 0–5 Hz firing, a separate histogram for 5–10 Hz firing, and so on, in 5 Hz steps, up to the maximum rate at which the cell fired.
We characterized ISI histograms by measuring their coefficient of variation (CV = s.d./mean) and also by fitting them to gamma distributions. Gamma probability density functions are defined by a shape (k) and scale (θ) parameter, as follows:
We estimated k and θ using a nonlinear least-squares procedure (trust-region algorithm) that minimized the error between the gamma density function and bin values of the ISI histogram (bin width: 1 ms). We excluded the first 8 bins of the ISI histogram in calculating the least squares measure so as to minimize the effects of the absolute refractory period and bursting on the shape of the fitted distributions. As we only examined ISI histograms associated with 10–40 Hz spiking (see Main Text), an 8 ms exclusion zone appeared to provide the best compromise between properly disregarding very short, burst-related or refractory-related, intervals, but also including enough data to properly estimate the overall shape of histograms. In early attempts, when we did not exclude the first 8 ms or when we used a maximum likelihood method, ISI histograms with burst-related peak followed by a trough and then a second peak (Figure 6A) were inappropriately fit with k ~ 1, i.e., as exponentials.
We generated a total of 3234 ISI histograms associated with 10–40 Hz spiking (539 cells × 6 histograms per cell). 2424 of these histograms (75%) had >150 interspike intervals and thus passed our threshold for fitting the histogram to a gamma distribution. We were most interested in k from the gamma fits, as this parameter provided a measure of spike-train regularity (see Main Text). Across all fits, k ranged between 0.25 and 21.3. 95% Confidence intervals for k had a mean width of 0.18 (s.d. 0.12), indicating that our estimates were accurate. Only 1/2424 fits had a 95% confidence interval for k wider than 1. We should note that we view gamma distributions as useful descriptors of ISI histograms; we do not claim any functional significance to the use of this density function as opposed to log-normal or inverse Gaussian distributions, for example, which may provide better fits (Barbieri et al., 2001). Future work will need to determine the process that gives rise to interspike intervals, thus constraining the proper statistical description of ISI histograms.
Importantly, for individual neurons, there was a good match between k and the qualitative appearance of ISI regularity in unprocessed spike trains (Supplemental Data). This correspondence, combined with extensive analysis of simulated spike trains (not shown), makes us confident that the results described are not an artifact of parsing interspike intervals into separate rate-matched histograms. Furthermore, any bias in the rate-parsing algorithm should have applied equally across cortical areas, whereas we observed different results in each area (Figure 4A). Note also that spike-count variance-to-mean plots (Figure 3B, right) were calculated without reference to ISI histograms and thus provide an independent line of evidence for ISI regularity.
We thank J.H.R. Maunsell, C. Padoa-Schioppa, and R. Maimon for helpful comments on the manuscript. K. Irwin, T. Lafratta, M. LaFratta, J. LeBlanc, and D. Averbuch provided technical assistance.
This work was supported by NEI EY12106 and the McKnight Endowment Fund for Neuroscience.
Publisher's Disclaimer: This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.