| HOME | HELP | FEEDBACK | SUBSCRIPTIONS | ARCHIVE | SEARCH | TABLE OF CONTENTS |
|
|
|
|||||||
|
|||||||||
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
INNOVATIVE METHODOLOGY
1The Clinical School, Addenbrooke's Hospital, Cambridge, United Kingdom; 2Department of Neuroscience, University of Pennsylvania, Philadelphia, Pennsylvania; and 3University of Newcastle upon Tyne, Sir James Spence Institute, Royal Victoria Infirmary, Newcastle upon Tyne, United Kingdom
Submitted 30 September 2005; accepted in final form 9 March 2006
 |
ABSTRACT |
|---|
 TOP ABSTRACT INTRODUCTIONMETHODSRESULTSDISCUSSIONGRANTSACKNOWLEDGMENTSREFERENCES |
|---|
|
INTRODUCTION |
|---|
|
TOPABSTRACTINTRODUCTIONMETHODSRESULTSDISCUSSIONGRANTSACKNOWLEDGMENTSREFERENCES |
|---|
; Softky and Koch 1993). In particular, it is not known whether irregularity can change in time in a similar way to firing rate. Such modulation of irregularity along a behavioral sequence could indicate the presence of temporal coding. However, measurement of irregularity from experimental data is a challenging problem. Further progress requires the development and application of new metrics and is the purpose of this paper.
Previous measurement of irregularity has largely been confined to evaluation of the coefficient of variation (CV) of the intervals, which is the ratio of the SD of the intervals to their mean. Physiological data on CV are available for the somatosensory system (Werner and Mountcastle 1963), muscle spindles (Burke et al. 1979; Nordh et al. 1983; Stein and Matthews 1965), auditory system (Stabler et al. 1996; Young et al. 1988), visual system (Holt et al. 1996), and motor system (Prut and Perlmutter 2003; Stern et al. 1997). Possible sources of spike timing irregularity have been sought at the membrane and synaptic levels both in experimental (Mainen and Sejnowski 1995) and modeling studies (Salinas and Sejnowski 2001; Stein 1965 for a review). Modeling studies have examined the effects of membrane excitability on CV and output interspike interval structure (Gutkin and Ermentrout 1998), and of statistical assumptions about multiple synaptic inputs with regard to excitation and inhibition balance (Christodoulou and Bugmann 2000, 2001; Feng and Brown 1998) and correlation (Feng and Brown 1999 2000; Salinas and Sejnowski 2000, 2002; Svirkis and Rinzel 2000).
However, the CV is not an ideal metric. We wish to study spike train irregularity when the data consist of some number of trials aligned to stimulus or behavioral events with the goal of a time-resolved measure conceptually similar to the peristimulus time histogram (PSTH) for firing rate. The CV could be used to characterize the distribution of intervals taken within a particular time bin and across trials. However, it will then measure not only the spike firing irregularity but also the extent to which firing rate varies across trials. Such trial-to-trial rate variation is a common feature of experimental data (Baker et al. 2001). In this paper, we examine several alternative metrics of irregularity that are independent of firing rate, only slightly affected by firing rate slope (both across and along trials) and have good time resolution. We characterize one such metric in considerable detail before applying it to data recorded from neurons in monkey motor cortex during a behavioral task.
|
METHODS |
|---|
|
TOPABSTRACTINTRODUCTIONMETHODSRESULTSDISCUSSIONGRANTSACKNOWLEDGMENTSREFERENCES |
|---|
We define metrics of irregularity in the context of spike train data taken under conditions of repeated trials involving stimulus or behavioral events. The data may then be represented as a raster where each trial is a line of spike times, and all trials are aligned on a particular selected stimulus or behavioral event. The time coordinate of the raster can later be binned. We seek metrics of irregularity the performance of which is comparable to that of a PSTH so that we can obtain a time-resolved average measure of irregularity relative to the alignment event. Obviously adequate estimation of any metric will involve averaging over some number of interspike intervals. Such averages can be taken along trials (horizontally in the raster); this loses time resolution relative to the alignment marker. Alternatively the average can be taken within a time bin but across trials (vertically in the raster); this preserves the desired time resolution but may create estimation difficulties if the trials are statistically different (for example in underlying firing rate). The desired irregularity metric should therefore be independent of firing rate and also insensitive to firing rate change, i.e., d(rate)/dt.
A metric that meets the above requirements and uses two adjacent intervals (Ii and Ii+1) is
 | (1) |
 | (2) |
We have examined a number of other metrics that are similarly independent of the average firing rate but less localized in time, so that modulation of rate either within or across trials can make the metric values problematic. An example of such a metric depends on a single interval (Ii) and local mean taken over n adjacent intervals on either side
 | (3) |
 | (4) |
Depending on details of individual datasets, there may be problems in calculating the local mean interval. Maximum time resolution would be attained by using three adjacent intervals (n = 1 in Eqs. 3 and 4; the irregularity metric then relates to the middle interval), and averaging across trials within time bins. Although not explicitly shown in this paper, in many stationary cases, the second irregularity metric produces results very similar to the first metric. However, we prefer the metric of Eqs. 1 and 2, since it is better localized and more rate insensitive. A further possible metric is that introduced by Shinomoto et al. (2003) to calculate a measure they refer to as LV
 | (5) |
Investigation of the measures we developed required generation of spike trains with known properties. For this purpose, we simulated point processes with gamma-distributed intervals. The gamma density has previously been used to describe interval distributions because it has two parameters that allow both rate and regularity (the "order" of a gamma distribution) to be independently varied. At an order value of one, the gamma interval distribution is exponential (as for a Poisson process, the most irregular case). At high order values >20, the distribution is narrow and almost symmetric about its mean interval (i.e., 1/rate); this corresponds to the intervals produced by a pacemaker neuron (the regular case). In this paper, the rate and/or order of the gamma distribution was modulated, using simulation methods previously described (Baker and Gerstein 2000).
To determine how the irregularity measure was affected by correlation between successive interspike intervals, we produced simulated spike trains as follows. N gamma distributed intervals were generated, and the median of the first M of these was determined. Denote ki as the indices of those intervals that are smaller than the median, and ji as the indices of those that are larger than or equal to the median. The intervals were then reordered according to
 |
 | (6) |
After characterizing the measure properties with simulated data, they were applied to data from experimental recordings. All calculations were carried out in the MATLAB language (Mathworks) on standard PC Windows computers.
Behavioral task
Experimental data in this paper were recorded from two female macaque monkeys (monkey T, 6.1 kg; monkey E, 5.4 kg), purpose-bred for research and trained to perform a precision-grip task for food reward. The behavioral task was a bilateral version of that used by Baker et al. (2001). The animal sat facing two precision-grip manipulanda. Access to these manipulanda was initially obstructed by clear plastic barriers. A trial was initiated by placing each hand on home pad switches, located just in front of these barriers. After a delay, a cue indicated whether left, right, or bimanual performance would be required on that trial; these three possibilities were chosen at random with equal probability. The cue (0.5- to 1-s duration) consisted of a 6-Hz vibration of the barrier and flashing light-emitting diodes (LEDs) on the relevant side or sides. After the cue was an instructed delay period (0.5- to 1-s long), during which the animal had to maintain both hands on the home pads. Both barriers then dropped allowing access to the manipulanda (the go cue). The animal was required to reach out with the cued hand or hands and squeeze the levers between finger and thumb in a precision grip. Both levers had to be moved more than a criterion distance and held with this displacement for 1 s before they could be released to obtain a food reward. The levers were attached to optical encoders that measured lever displacement and torque motors that provided a resisting force. Under computer control, the motors simulated the action of a spring. A 0.15-N force was required to move each lever from its end stop; further movement produced increased resisting force with a spring constant of 0.03 N/mm of lever displacement. Auditory feedback on task performance was provided using tones that indicated the different task phases. Early release of a home pad switch, or movement of an arm that had not been cued, resulted in failure of that trial and a return to the initial condition. In this paper, we report only data from unimanual trials. We refer to the trials as contralateral or ipsilateral, indicating whether the hand contralateral or ipsilateral to the recorded neurons was moving. All analysis was aligned to the end of the hold phase of the task. Trials were discarded from analysis if the animal moved out of target prior to successful completion of the hold phase, if the delay between the go cue and entry into the target zone was >0.95 s, or if the delay between the end of the hold phase and lever release was >0.65 s.
Surgical preparation
Once behavioral training was complete, the animals were implanted with a stainless steel headpiece to allow atraumatic head fixation (Lemon 1984). A recording chamber was placed over a craniotomy centered at stereotaxic coordinates A13 L18, providing access to the hand representation of primary motor cortex. Finally, fine stimulating electrodes (Microprobe, type LF01G) were inserted into the pyramidal tract (PT) at the level of the medulla, to permit antidromic identification of pyramidal tract neurons (PTNs) (see Baker et al. 1999; Lemon 1984). After recordings were complete in this chamber, a second chamber was implanted at the same coordinates on the other side, with corresponding PT electrodes, and further data were gathered from that hemisphere. All surgical operations were performed under deep general anesthesia (2–2.5% isoflurane in 50:50 O2:N2O) using aseptic technique and were followed by a full course of antibiotic (co-amoxiclav 140/35, 1.75 mg/kg clavulanic acid, 7 mg/kg amoxycillin, Synulox, Pfizer) and analgesic (buprenorphine; Vetergesic, 10 µg/kg, Reckitt and Colman Products) treatment. At the end of the experiment, the animals were killed by overdose of anesthetic and perfused through the heart with fixative. Subsequent histological processing of the tissue confirmed that recording and stimulating electrodes had been correctly located. All procedures were carried out under appropriate licenses from the UK Home Office.
Single-unit recordings
Single-unit activity was recorded from the motor cortex using a 16-channel "Eckhorn" microdrive (Eckhorn and Thomas 1993; Thomas Recording, Giessen, Germany). After amplification (2–10 K gain) and filtering (band-pass 300 Hz to 10 kHz), signals were continuously sampled at 25 kHz and stored for off-line analysis. Cells that responded antidromically to PT stimulation were confirmed as PTNs using a collision test (Baker et al. 1999). Spike waveforms were discriminated into the occurrence times of single-unit action potentials using custom-written cluster cutting software. Only units with consistent spike waveforms and no interspike intervals <1 ms were used for subsequent analysis.
|
RESULTS |
|---|
|
TOPABSTRACTINTRODUCTIONMETHODSRESULTSDISCUSSIONGRANTSACKNOWLEDGMENTSREFERENCES |
|---|
order. The difference between the utility of IR and CV becomes apparent when we consider the estimation problem. Estimation of mean and SD of intervals, as required for CV, requires that the sample intervals be drawn from a process of constant rate. However, in real data, firing rate varies both along and across trials, making it generally impossible to estimate CV. By contrast, estimation of IR requires only that firing rate be constant during the two adjacent intervals used in the m metric; it is essentially a local measure in time that is not affected by firing rate variations within or across trials. Problems with IR estimation do arise if there are very fast rate variations along a trial so that rate constancy over two adjacent intervals is violated; we will examine that situation below and in Fig. 2.
|
|
To satisfy our goal of a time resolved irregularity measure, we will average the irregularity metric values for intervals within a particular time bin over many aligned trials. Any trial-to-trial variation of underlying rate will not affect this local estimate of the IR. The confidence limits of the IR estimate [mean(m) in the bin] will obviously depend on N, the number of interval pairs across trials in the particular time bin. A complicating factor, however, is that in computing the IR we use overlapping interval pairs (1,2; 2,3; 3,4;...), thus introducing some dependence into successive samples. Although detailed numbers and effects of this dependence will vary from bin to bin, on average there are twice as many interval pairs in a bin with the overlap calculation as there would be with a calculation that used independent samples by jumping (1,2; 3,4; 5,6;...). We evaluated the upper and lower confidence limits for the overlap case using simulation. Figure 1E shows these results for 95 and 99% confidence limits as a function of N, the number of interval pairs in the bin. Although the figure shows results using data simulated with gamma order 4, the curves were similar for other orders. The figure is expressed as multipliers of the quantity Q = SD(m)/N for the bin. If we were dealing with independent and normally distributed values, this quantity would be the SE for N independent events, and the multiplier would be 1.98 for the 95% confidence limit. In the present situation, N is twice what it would be for independent interval pairs, and the SD of m is smaller. Q is therefore smaller than the corresponding SE for independent interval pairs. Thus although the multiplier values in Fig. 1E are larger than the values appropriate for the independent situation, the overall confidence limits will be smaller than if nonoverlapping interval pairs had been used. The curve for the upper 95% confidence limit multiplier is well fitted by [1.21/log(N) + 2.29]. The multiplier for the corresponding lower confidence limit has a constant value of 2.25 for n > 5. For gamma orders 4 and 8, the actual upper and lower 95% confidence limits as a function of N per bin are shown in Fig. 1F. The asymmetry of the PDF mean distribution in Fig. 1D underlies the difference between upper and lower confidence limits. Figure 1F provides an estimate of the sensitivity of the measure: it shows the smallest change in IR that will be detectable given a certain size of dataset.
An alternative method of computing IR would be to average over a fixed number of values rather than using however many values fall in a fixed time bin. This would mean that the temporal resolution changed depending on the firing rate: during periods of high firing rate, a narrow time window would be used, whereas at low rates, the temporal smearing would be greater. However, unlike the current method, the confidence limits on IR would be of relatively constant size throughout. We prefer the current approach using a fixed bin width, rather than using a fixed number of intervals, because it is straightforward to show a variable-sized confidence limit on a standard line plot. It is much harder to show clearly a variable time window associated with each IR estimate.
Although the irregularity metric is independent of firing rate (Fig. 1A), because it measures a ratio of two adjacent intervals, it is possible that its value could depend on the derivative of rate d(rate)/dt. This is addressed in Fig. 2, where we illustrate the effects on the metric of different profiles of rate change—in all cases, data were simulated using a Gamma process of constant order 4 so that irregularity of firing was actually fixed throughout. Figure 2A shows a step change in rate from 20 to 200 Hz. The IR measure (Fig. 2B) is artifactually increased in the single bin which encompasses the rate step (thick line). Note that the confidence limits are narrower after the step, as the higher firing rate yields more spikes per bin.
Figure 2C illustrates a ramp rise in firing rate (from 20 to 400 Hz in 0.5 s). The time resolved IR (Fig. 2D, thick line) is essentially flat except for a positive error in the single bin near the ramp onset. The critical variable here is the underlying rate change over the time of the two adjacent intervals used for evaluating the m metric. Thus although the rate derivative is constant throughout the ramp, the nonstationarity only perturbs the metric significantly at low firing rates. The thin lines in Fig. 2, B and D, shows the result of applying a correction procedure that will be described in the following text. For the step rate change, this effectively compensates for the error; for the ramp change, it results in an underestimate of the irregularity.
Figure 2, A—D, pertains to IR values in a bin near a rate disturbance. To obtain greater insight into how rate changes effect the IR measure, we examined the metric m obtained in the immediate vicinity of a step rate change. We evaluated m for each of the two interval pairs that involve the distorted interval, one pair with the last normal interval before the step, one with the first normal interval after the step (see Fig. 2E). Averaging many such step events produces Fig. 2F, which shows the corresponding error in the IR measure versus the ratio of the rate before to the rate after the step. The error has been calculated as the difference between the IR measured from the distorted intervals (Fig. 2E), and the IR which is measured from pairs of intervals generated with stationary rate and the same gamma order (as in Fig. 1A). The top curve shows results for the interval pair on the high rate side of the step, and the bottom for the pair on the low rate side of the step; their mean is also shown. The ordinate is the error. The error increases as the step size increases and is larger on the high rate side of the step. We have verified that the curves in Fig. 2F are independent of the gamma order used in the simulation and of the absolute rate values used for the step; only the rate ratio is relevant.
As an aside, note that even at a step rate ratio of one (corresponding to no rate change), there is a small positive error in m. This is an artifact of the method used to calculate this figure. The middle interval (I2 in Fig. 2E) has been selected as the one that straddles the assumed rate step. The process of selecting an interval based on a fixed time biases the distribution toward longer intervals. The intervals on either side (I1 and I3 in Fig. 2E) will show no such bias, and the result is a slightly elevated value of the metric m, even if there is no actual rate step. This effect is not seen when we average all metric values falling within a certain bin, which is how we will normally use the irregularity measure.
A further potential source of error in the IR measure is if successive intervals are correlated because the metric m depends on the joint probability distribution of Ii and Ii+1 (Eq. 2). To assess the effect of this, we simulated data with different levels of serial correlation, using the interval reordering method described in METHODS. Figure 2G shows a joint-interval scatter plot for such a simulation using gamma order 8, in which the first 50% of the intervals had been reordered. For these data, the correlation between successive intervals so produced was r = –0.30. Figure 2H illustrates the effect of different levels of interval correlation on the IR measure for spike trains generated with varying underlying regularity. Increasing the magnitude of the correlation elevates the estimate of IR; the numerical size of this effect is greatest for the most irregular spike train (gamma order 2). Note, however, that even with interval correlations spanning the relatively high range from r = 0 to –0.3, the curves relating to gamma orders 2, 4, 8, and 16 have nonoverlapping values of IR.
Figure 2 indicates that the IR measure can be substantially affected by rate steps; rate ramps by contrast are less problematic except when they start from low initial rates. Several features of experimental data will act to lessen the IR errors caused by such rate nonstationarities. First, because metric values are averaged in a bin, there will be a "dilution" effect: metrics corrupted by rate changes will be only a fraction of those averaged within the bin. This dilution will be greater at high firing rates as more uncorrupted intervals will fall within a single bin. Second, real data are also likely to exhibit latency jitter, such that the timing of the rate step is not precise relative to the alignment marker. This will further dilute the IR error. However, it is clear that some caution must be exercised in interpreting IR values in the presence of rapid rate changes.
Given the steep modulations in rate that are commonly observed in experimental data, it would be useful to develop a correction for the distorted IR values near firing rate steps. Such a correction is demonstrated in Fig. 3. Figure 3A shows the IR error within a 100-ms-wide bin for simulated data at various gamma orders and rate steps of various sizes. The curves initially rise steeply but then diminish once the rate rise is greater than fivefold. The reason for the decline is dilution: at the high rates there will be more "normal" intervals (i.e., drawn from periods with stationary rate) and hence more "normal" metrics m in the bin.
|
The one-sided distortion of the PDF of ms shown in Fig. 3B suggests a simple way to correct the IR error for bins containing a rate step. Whenever the measure S exceeds some positive or negative threshold, the IR will be computed using only the intervals and ms metric values that fell into the unscathed half of the ms PDF. Thus for the 5x step condition shown in Fig. 3B we would compute the IR using only interval pairs corresponding to the right half of the PDF. In each case we take the absolute value of this estimate. The procedure reduces the number of metric values to be averaged, usually by over half, and therefore the confidence limits at such a bin will be wider than at adjacent bins where the correction is not necessary.
The correction procedure strongly reduces the IR error as shown for various gamma orders and step sizes in Fig. 3E. Note that the maximum span of the error is about 1/10 of that for the uncorrected IR values in Fig. 3A. The choice of threshold value on S to trigger the correction procedure is arbitrary, but Fig. 3F demonstrates that the correction procedure is relatively insensitive to the chosen value.
Application of the IR measure to experimental data recorded from monkey motor cortex during performance of a precision grip task is shown in Figs. 4–6. In each case, trials have been aligned on the end of the hold phase of the task, indicted by zero on the time scale (see METHODS). A bin width of 0.1 s was used, chosen to be sufficiently small to represent the profile of IR modulation but also large enough to yield narrow confidence limits on the IR estimates. The columns show averaged finger and thumb position trace (Fig. 4A), the spike raster (B), PST histogram (C), time-resolved S measure (D), time-resolved IR measure (with and without correction, thick and thin lines respectively, E), and comparison of the interval histograms from shaded regions of the IR graph (F). The shaded region in Fig. 4D marks the threshold for IR correction. When S falls within this range, the IR measure is the average of all available m metrics. If S is lower than threshold, only positive m values contribute to the corrected IR; if S is above threshold, only negative m values contribute. In either case, the number of m values available for averaging is reduced; the confidence limits for the corrected bins are therefore wider than the adjacent uncorrected bins. For this cell, there were profound step changes in firing rate, which produced large artifactual changes in the uncorrected IR value around 1.1 s before the end hold marker. These were successfully removed by the correction procedure. The time course of firing rate (Fig. 4C) was quite different from the time course of the IR (Fig. 4E). Further confirmation that firing regularity did indeed change in a task-dependent way was provided by examining interval histograms from different regions of the task (shaded regions in Fig. 4E). Figure 4F plots these histograms; the abscissa is log(I) – log(I), which effectively normalizes the intervals for firing rate. The interval distributions clearly differ in shape: intervals from region (a) in Fig. 4E, where IR is low, are more narrowly distributed than those from the region with high IR (b).
|
|
|
A population analysis of the available dataset is shown in Fig. 7. Figure 7A plots the mean IR modulation during contralateral hand performance calculated over all 108 cells recorded. This shows that firing became on average more regular during the hold phase of the task. However, the single-cell examples already illustrated in Figs 4–6 demonstrate that such a simple population mean hides a more complex underlying pattern. To analyze these patterns of IR modulation across different cells in more detail, we first determined for each cell whether IR modulated significantly during contralateral task performance. To achieve this, for a single cell's IR modulation, all possible pairs of IR values were compared, and the number counted that were significantly different (as judged from their calculated confidence limits). Because there were 30 IR values (a 3-s analysis window, using 0.1-s bins), there were 435 pairwise comparisons. Using simulation data in which regularity did not change, we determined that this count of significantly different pairs would not exceed 50 more than 5% of the time by chance. Sixty-one cells had a pair count >50 and were therefore assumed to have a significantly task-modulated IR value (P < 0.05). Figure 7B shows the average IR for the remaining 47 cells that were not significantly modulated. This average still shows a small modulation, consistent with that from the significantly modulated population (Fig. 7C). It is likely therefore that some of the "unmodulated" cells did in fact change their firing regularity with the task, but that it was too small relative to the noise level to be detectable.
|
) and unidentified cells (
) are shown. Although no clear clustering was apparent, it is instructive to view how rate and the IR measure modulate for different subsets of the population indicated by the rather arbitrary boundaries in 7G. This allows a concise overview of the diversity of patterns in modulation across the population, without exhaustively presenting results from all cells. The panels of Fig. 7, H and I, show this; each plot corresponds to different regions (
) in Fig. 7G. The firing rate patterns from Fig. 7I, 1–4 are fairly similar with elevated firing during the movements before and after the hold portion of the task. The associated mean IR patterns in Fig. 7H, 1–4, show considerable variety: irregularity may be lower during hold than during pre- and post-movements (H, 1–3) or more rarely, irregularity may be higher during hold (H4). The rises in IR values around time –1 s also show considerable latency differences that are not in step with the associated latency differences in the firing rate panels. Thus in H1 and I1, the initial IR peak is earlier than the rest of the population, whereas the initial firing rate peak is later. The reverse situation pertains for H2 and I2. Clearly the precise details of these plots would be altered by changing the location of the boundaries placed on Fig. 7G. However, the key conclusion is that irregularity changes do not simply mirror the rate modulation. A similar analysis was carried out for the IR modulation during ipsilateral task performance. Only 15/108 cells showed a significantly modulated IR measure. The mean IR (Fig. 7J) and mean firing rate (Fig. 7K) across the entire population of recorded cells showed little modulation.
It was of interest to determine whether the reduced irregularity observed during the hold phase of the precision grip was specific to this task or occurred during other periods of steady contraction. We accordingly also calculated the IR measure for spikes falling between the time when the home pads were first pressed and the time when the go cue was given for trials performed with the contralateral hand. The mean IR for this phase, over all 108 cells, was 0.57. This compared with a mean IR of 0.63 during the hold phase of the precision grip (taken as the 0.8-s period prior to the end-hold marker). In only 1/108 cells were the two estimates of IR significantly different as judged from their 95% confidence limits, a number well within that expected by chance. The irregularity therefore appeared similar during these two different periods of steady contraction.
|
DISCUSSION |
|---|
|
TOPABSTRACTINTRODUCTIONMETHODSRESULTSDISCUSSIONGRANTSACKNOWLEDGMENTSREFERENCES |
|---|
The concept of the irregularity of neural discharge is intuitively straightforward; its experimental measurement is more problematic. Many previous studies have used the CV of the interspike interval to measure irregularity. This gives good results if neural firing rate is constant. However, if rate changes over the measurement period, the CV will include a contribution from rate variation as well as spiking irregularity. The problem is especially acute in data from awake behaving animals when firing rate shows not only modulation during task performance but also considerable variability from trial to trial (e.g., Baker et al. 2001).
The novel irregularity metric m and the IR measure that we introduce here have several benefits in analyzing experimental data. The m metric has good temporal resolution and is independent of firing rate. Values of m may therefore be validly averaged over periods with very different firing rates. We have used this property in the analysis of our experimental data: m values were binned relative to a behavioral marker, and averaged across trials to yield a time resolved IR measure. Trial-by-trial variation in discharge rate could not affect the result.
We have found only two situations where the IR measure can be corrupted. First, over time periods where there is a rapid change in the firing rate, the IR value may be artifactually increased. Slow changes in rate from a moderate baseline do not appreciably alter the measure (Fig. 2) because the variability in interspike interval is usually large compared with the rate trend. The problem comes with very rapid, near stepwise rate changes particularly from low firing rates (Fig. 2). To overcome this problem, we have introduced an additional procedure based on the unsigned ms metric, which identifies and largely compensates for such errors.
Second, if there is a correlation between successive interspike intervals, this can elevate the IR measure (Fig. 2G). Although experimental data often show a nonzero joint interval correlation, this is usually positive and caused by the profound task-dependent rate modulations. Short intervals occur during periods of high firing rate and hence are usually surrounded by short intervals; likewise, long intervals come from task phases with low firing rate and are predominantly preceded or succeeded by long intervals. Such correlation due to rate modulation will not effect IR, which is robust even if rate changes greatly. No current method exists to measure experimental interval correlation while correcting for modulating rate, and so it remains uncertain if the sensitivity of IR to correlated intervals is a practical concern. However, it could be argued that negative interval correlation (as in Fig. 2, G and H) is just a further form of irregularity. In this case, the fact that it elevates IR could be seen as an advantage of the method.
The time-resolved nature of the IR measure makes it well suited to analyses that seek to generate "surrogate" data; such surrogates are used to define the expected results of a particular novel analysis method. Previous work (Baker and Lemon 2000) has used surrogates with a firing rate that modulated as in the experimental data, and gamma distributed intervals of fixed regularity. These surrogates were then used to determine the number of spatiotemporal repeating patterns expected in a given dataset by chance. Our present results demonstrate conclusively that irregularity is not fixed (in a similar task to that used by Baker and Lemon 2000) but modulates during task performance. The instantaneous irregularity measure could be used to refine the gamma-based surrogate generation process so that surrogates matched not only the single-trial firing rates of the experimental data but also the trial-averaged modulation in irregularity. We note that other methods of surrogate generation based on spike time dithering (Date et al. 1998, 2000) automatically preserve both rate and irregularity profiles (G. L. Gerstein, unpublished observations). How better surrogates would affect the predicted number of spatiotemporal repeating patterns (if at all) remains to be seen (Hatsopoulos et al. 2001).
An alternative potential method for assessing irregularity is to fit the observed spike train to a parametric model (e.g., a gamma process as in Baker and Lemon 2000). Once optimal model parameters have been determined (Kass et al. 2005; Truccolo et al. 2005), their values directly measure regularity. During the exploratory phase of this project, we initially tried several such model-based approaches. None were as practically effective as the metric-based analysis reported in this paper, and all were computationally more complex. In addition, it is theoretically more attractive to carry out analysis in a model-free framework if possible, because even relatively good models such as the gamma process are unlikely to represent fully the subtleties of an experimental spike train.
Causes of task-dependent changes
The high variability in interspike intervals that is commonly seen in central neurons results from random fluctuations in synaptic inputs (Shadlen and Newsome 1998; Softky and Koch 1992; Stein 1965). During the precision grip task performed with the contralateral hand, many cells showed more regular firing during the hold phase compared with the movement phase. This probably implies a reduction in the transmembrane synaptic voltage noise during holding compared with movement. Two distinctly different factors could contribute to produce this.
First, the amount of synaptic input occurring per unit time could be reduced. The average firing rate of many of these neurons was reduced during the hold phase (Figs. 4–7), so that a reduction in excitatory drive to the motor cortical cells during this period appears reasonable. However, the synaptic noise will depend on the total number of synapses active, both excitatory and inhibitory; it is not clear how the level of inhibitory input might change during the hold phase. Second, the voltage noise will depend on the input conductance of the neuron, which determines the size of voltage excursion generated by a single synaptic activation. It is conceivable that the complex interplay of voltage-gated channels responsible for the afterhyperpolarization (AHP) in these cells could produce a change in input conductance during the time of lowest firing rate.
Two other factors could contribute to the observed changes in firing irregularity. Wetmore and Baker (2004) and Chen and Fetz (2005) recently showed that many motor cortical neurons have a depolarizing peak in their membrane potential 10–60 ms after a spike. At both low and high rates, spiking from such a cell will be relatively unconstrained, yielding irregular firing. However, if the cell fires at intermediate rates, the mean interval will coincide with the postspike potential peak. In this situation, spikes will occur preferentially over the narrow range of intervals around the peak, resulting in more regular discharge. The average firing rate often fell to 
25 Hz during the hold phase (Figs. 4–7). This is in the correct range for intervals to align with the membrane potential peaks, and the low irregularity seen during the hold phase may be partly caused by this biophysical mechanism.
However, not all of the modulations in irregularity seen can be explained as a biophysical byproduct of the profound rate modulations that occur during task performance. For example, a subset of cells showed a significant increase in irregularity during the hold phase of the task (Figs. 5 and 7); the rate profile for these cells, as for many others, was low during the hold compared with the movements before and after.
An alternative explanation for irregularity changes is that the cortical network has entered a different state. During steady contractions, the motor cortex can generate synchronous bursts of oscillations, at frequencies around 10, and 15–35 Hz (Baker et al. 1997; Kilner et al. 2000; Murthy and Fetz 1992; Stancak and Pfurtscheller 1996). During the hold phase of the precision grip task, a considerable fraction of the input to motor cortical neurons is synchronized to these network oscillations (Baker et al. 2003). As mentioned in the preceding text, the average neural firing rate in the hold phase is 25 Hz. Cells that spike around the peak of one network oscillation cycle will therefore be coming close to threshold just as the next oscillation peak arrives. The range of intervals will then be contracted to lie close to the oscillation period, and the discharge will become more regular. Interestingly, a similar modulation in network oscillation amplitude is seen during performance of a precision grip with the ipsilateral as the contralateral hand (Kilner et al. 2003). Our data, however, show far less modulation of irregularity in this situation (Figs. 5 and 7). This implies that network oscillations are not the only explanation for the regularity modulation which we report.
Impact of irregularity changes on downstream neurons
The changes in irregularity that we have observed are interesting in their own right as they reflect the changing state of the cortical network. However, they may also have important effects on downstream neurons receiving synaptic input from these cells. Irregularity is unlikely to exert an effect when downstream neurons simply sum over many inputs. The union of a large number of spike trains with gamma-distributed inter-spike intervals tends toward a Poisson process (i.e., highly irregular), regardless of the regularity of the summed inputs (Tiesinger 2004). However, the synapses made by the cortical cells may display history-dependent modulation of the postsynaptic potential amplitude. In this case, irregularity will become important, as it will influence the relative proportions of long and short intervals seen at a given firing rate.
One important class of neurons receiving input from the cells reported here is other cortical cells, both locally and in distant areas. Cortical neurons show an exquisite repertoire of history-dependent modulation of synaptic efficacy (reviewed in Thomson 2003), including facilitation, augmentation, depression, and selective enhancement after inter-spike intervals falling in a narrow window (the "notch filter"). The relative importance of these history effects depends on both the type of pre- and postsynaptic cell. Additional complications to the effectiveness of a particular cortico-cortical connection can come from the "state" of the postsynaptic cell as influenced by its other inputs (Petersen et al. 2003; Sachdev et al. 2004). The modulation in irregularity reported here is therefore very likely to alter the firing of other cortical neurons.
In addition, the motor cortical neurons that modulated their irregularity included identified PTNs, which make monosynaptic connections to spinal motoneurons (Porter and Lemon 1993). Motoneurons show a facilitation of synaptic inputs for short interspike intervals (Porter 1970), although the reported time constant of this effect is 10 ms. Because the PTNs that we report generally fire at rates <100 Hz, there are usually few very short intervals for either regular or irregular neurons. An exception would be when bursting is present (as in the unidentified cell illustrated in Fig. 5); changes in irregularity caused by alterations in the burst probability will almost certainly have an effect on the magnitude of excitatory postsynaptic potentials (EPSPs) produced in target motoneurons.
One report (Lemon and Mantel 1989) suggested that the corticomotoneuronal synapse might show selective facilitation for interspike intervals from 40 to 70 ms. If correct, the modulation in irregularity that we have observed could greatly influence motoneurons. The mean firing rate during the hold phase of our task is 25 Hz (Fig. 7I), and changes in irregularity will alter the fraction of interspike intervals in the enhanced 40- to 70-ms range. However, Lemon and Mantel (1989) pointed out that their indirect results could also be explained by cortical synchronization rather than a property of the corticomotoneuronal synapse. Direct evidence for a "notch filter" at the corticomotoneuronal synapse is lacking. Jackson et al. (2006) recently reported spinal monosynaptic field potentials after paired corticospinal tract stimulation; these did not show any significant enhancement for intervals 50 ms. It therefore remains uncertain whether changes in irregularity could have functional significance for motoneuron activation.
Importance of irregularity changes for cortical computation
There are currently two conflicting theories of how neurons encode and processes information. Conventionally, it is assumed that a neuron's firing rate encodes information (Lettvin et al. 1959). To read this code, a downstream neuron must effectively count spikes. Irregular firing introduces noise, requiring averaging over longer counting windows and multiple neurons to decode the information accurately (Shadlen and Newsome 1998). On this view, the periods of low firing irregularity are when the neuron is best able to transmit information with high signal-to-noise ratio.
By contrast, more recent work has suggested that precise spike timing may carry information (Abeles 1991; Gerstein et al. 1989). The observed variability in the interval distribution then represents not superfluous noise but temporal encoding of the signal. Highly regular neural discharge then corresponds to a limited encoding space, and more information can be transmitted if a wide range of interspike intervals are permitted (irregular discharge). The relatively irregular discharge seen in cortical neurons has been used as evidence to support both the temporal and rate coding hypotheses by different authors (Shadlen and Newsome 1994; Softky and Koch 1992).
The predicted consequences of the irregularity changes observed here therefore depend on the model of neural representation that is adopted. Although our present data cannot directly resolve which model is correct, the task relationship of the irregularity changes is of interest. Our task incorporates complex and highly skilled motor acts, involving targeted reaching, fractionated hand shaping, and lever squeezing. After lever release, the hand is reshaped and reaches out for the reward. During these most-demanding phases of the task, the motor cortical discharge is on average most irregular. It becomes regular only during the hold phase, which is by contrast straightforward: the animal must merely maintain a steady contraction for a criterion time that leads to reward. Our results therefore associate periods of irregular discharge with complex motor processing. This is contrary to the prediction of the rate coding hypothesis but agrees with what would be expected if the motor cortex were using temporal coding during the solution of this demanding motor problem.
|
GRANTS |
|---|
|
TOPABSTRACTINTRODUCTIONMETHODSRESULTSDISCUSSIONGRANTSACKNOWLEDGMENTSREFERENCES |
|---|
|
ACKNOWLEDGMENTS |
|---|
|
TOPABSTRACTINTRODUCTIONMETHODSRESULTSDISCUSSIONGRANTSACKNOWLEDGMENTSREFERENCES |
|---|
|
FOOTNOTES |
|---|
Address for reprint requests and other correspondence: S. N. Baker, University of Newcastle upon Tyne, Sir James Spence Institute, Royal Victoria Infirmary, Newcastle upon Tyne NE1 4LP, UK (E-mail: stuart.baker{at}ncl.ac.uk)
|
REFERENCES |
|---|
|
TOPABSTRACTINTRODUCTIONMETHODSRESULTSDISCUSSIONGRANTSACKNOWLEDGMENTSREFERENCES |
|---|
Baker SN and Gerstein GL. Improvements to the sensitivity of gravitational clustering for multiple neuron recordings. Neural Comput 12: 2597–2620, 2000.[CrossRef][Web of Science][Medline]
Baker SN and Lemon RN. Precise spatiotemporal repeating patterns in monkey primary and supplementary motor areas occur at chance levels. J Neurophysiol 84: 1770–1780, 2000.
Baker SN, Olivier E, and Lemon RN. Task dependent coherent oscillations recorded in monkey motor cortex and hand muscle EMG. J Physiol 501: 225–241, 1997.
Baker SN, Philbin N, Spinks R, Pinches EM, Wolpert DM, MacManus DG, Pauluis Q, and Lemon RN. Multiple single unit recording in the cortex of monkeys using independently moveable microelectrodes. J Neurosci Methods 94: 5–17, 1999.[CrossRef][Web of Science][Medline]
Baker SN, Pinches EM, and Lemon SN. Synchronisation in monkey motor cortex during a precision grip task. II. Effect of oscillatory activity on corticospinal output. J Neurophysiol 89: 1941–53, 2003.
Baker SN, Spinks R, Jackson A, and Lemon RN. Synchronization in monkey motor cortex during a precision grip task. I. Task-dependent modulation in single-unit synchrony. J Neurophysiol 85: 869–885, 2001.
Burke D, Skuse NF, and Stuart DG. The regularity of muscle spindle discharge in man. J Physiol 291: 277–290, 1979.
Chen D and Fetz EE. Characteristic membrane potential trajectories in primate sensorimotor cortex neurons recorded in vivo. J Neurophysiol 94: 2713–2725, 2005.
Christodoulou C and Bugmann G. Near Poisson-type firing produced by concurrent excitation and inhibition. Biosystems 58: 41–48, 2000.[CrossRef][Web of Science][Medline]
Christodoulou C and Bugmann G. Coefficient of variation (CV): vs mean interspike interval (ISI) curves: what do they tell us about the brain? Neurocomput 38- 40: 1141–1149, 2001.
Date A, Bienenstock E, and Geman S. On the temporal resolution of neural activity. Technical Report, Division of Applied Mathematics, Brown University. www.dam.brown.edu/people/elie/publications.html, 1998.
Date A, Bienenstock E, and Geman S. A statistical technique for the detection of fine temporal structure in multi-neuronal spike trains. Soc Neurosci Abstr 25: 1411, 1999.
Date A, Bienenstock E, and Geman S. A statistical tool for testing hypothesis about the temporal resolution of neural activity. Soc Neurosci Abstr 26: 828.6, 2000.
Eckhorn R and Thomas U. A new method for the insertion of multiple microprobes into neural and muscular tissue, including fiber electrodes, fine wires, needles, and microsensors. J Neurosci Methods 49: 175–179, 1993.[CrossRef][Web of Science][Medline]
Feng J and Brown D. Impact of temporal variation and the balance between excitation and inhibition on the output of the perfect integrate-and-fire model. Biol Cybern 78: 369–376, 1998.[CrossRef][Web of Science][Medline]
Feng JF and Brown D. Coefficient of variation of interspike intervals greater than 0.5. How and when? Biol Cybern 80: 291–297, 1999.[CrossRef][Web of Science][Medline]
Feng J and Brown D. Impact of correlated inputs on the output of the integrate- and-fire model. Neural Comput 12: 671–692, 2000.[CrossRef][Web of Science][Medline]
Gerstein GL, Bedenbaugh P, and Aertsen AMHJ. Neuronal assemblies. IEEE Trans Biomed Eng 36: 4–14, 1989.[CrossRef][Web of Science][Medline]
Gutkin BS and Ermentrout GB. Dynamics of membrane excitability determine interspike interval variability: a link between spike generation mechanisms and cortical spike train statistics. Neural Comput 10: 1047–1065, 1998.[CrossRef][Web of Science][Medline]
Hatsopoulos NG, Amarasingham A, Bienenstock E, Geman S, and Donoghue J. Assessing precise temporal patterns of spikes among cortical neurons. Soc Neurosci Abstr 27: 63.1, 2001.
Holt GR, Softky WR, Koch C, and Douglas RJ. Comparison of discharge variability in vitro and in vivo in cat visual cortex neurons. J Neurophysiol 75: 1806–1814, 1996.
Jackson A, Baker SN, and Fetz EE. Tests for presynaptic modulation of corticospinal terminals from peripheral afferents and pyramidal tract in the macaque. J Physiol 573: 107–120, 2006.
Kass RE, Ventura V, and Brown EN. Statistical issues in the analysis of neuronal data. J Neurophysiol 94: 8–25, 2005.
Kilner JM, Baker SN, Salenius S, Hari R, and Lemon RN. Human cortical muscle coherence is directly related to specific motor parameters. J Neurosci 20: 8838–8845, 2000.
Kilner JM, Salenius S, Baker SN, Jackson A, Hari R, and Lemon RN. Task-dependent modulations of cortical oscillatory activity in human subjects during a bimanual precision grip task. Neuroimage 18: 67–73, 2003.[CrossRef][Web of Science][Medline]
Lemon RN. Methods for Neuronal recording in Conscious Animals. IBRO Handbook Series: Methods in Neurosciences. London: Wiley, 1984.
Lemon RN and Mantel GWH. The influence of changes in the discharge frequency of corticospinal neurons on hand muscles in the monkey. J Physiol 413: 351–378, 1989.
Lettvin JY, Maturana HR, McCulloch WS, and Pitts WH. What the frog's eye tells the frog's brain. Proc Inst Radio Eng 47: 1940–1951, 1959.
Mainen ZF and Sejnowski TJ. Reliability of spike timing in neocortical neurons. Science 268: 1503–1506, 1995.
Murthy VN and Fetz EE. Coherent 25- to 35-Hz oscillations in the sensorimotor cortex of awake behaving monkeys. Proc Natl Acad Sci USA 89: 5670–5674, 1992.
Nordh E, Hulliger M, and Vallbo AB. The variability of inter-spike intervals of human spindle afferents in relaxed muscles. Brain Res 271: 89–99, 1983.[CrossRef][Web of Science][Medline]
Petersen CC, Hahn TT, Mehta M, Grinvald A, and Sakmann B. Interaction of sensory responses with spontaneous depolarization in layer 2/3 barrel cortex. Proc Natl Acad Sci USA 100: 13638–13643, 2003.
Porter R. Early facilitation at corticomotoneuronal synapses. J Physiol 207: 733–745, 1970.
Porter R and Lemon RN. Corticospinal Function and Voluntary Movement. Physiological Society Monograph Series. Oxford, UK: Oxford Univ. Press, 1993.
Prut Y and Perlmutter SI. Firing properties of spinal interneurons during voluntary movement. I. State-dependent regularity of firing. J Neurosci 23: 9600–9610, 2003.
Sachdev RN, Ebner FF, and Wilson CJ. Effect of subthreshold up and down states on the whisker-evoked response in somatosensory cortex. J Neurophysiol 92: 3511–3521, 2004.
Salinas E and Sejnowski TJ. Impact of correlated synaptic input on output firing rate and variability in simple neuronal models. J Neurosci 20: 6193–6209, 2000.
Salinas E and Sejnowski TJ. Correlated neuronal activity and the flow of neural information. Nat Rev Neurosci 2: 539–550, 2001.[CrossRef][Web of Science][Medline]
Salinas E and Sejnowski TJ. Integrate-and-fire neurons driven by correlated stochastic input. Neural Comput 14: 2111–2155, 2002.[CrossRef][Web of Science][Medline]
Shadlen MN and Newsome WT. Noise, neural codes and cortical organization. Curr Opin Neurobiol 4: 569–579, 1994.[CrossRef][Medline]
Shadlen MN and Newsome WT. The variable discharge of cortical neurons: implications for connectivity, computation, and information coding. J Neurosci 18: 3870–3896, 1998.
Shinomoto S, Shima K, and Tanji J. Differences in spiking patterns among cortical neurons. Neural Comput 15: 2823–2842, 2003.[CrossRef][Web of Science][Medline]
Softky WR and Koch C. Cortical cells should fire regularly, but do not. Neural Comput 4: 643–646, 1992.[CrossRef][Web of Science]
Softky WR and Koch C. The highly irregular firing of cortical cells is inconsistent with temporal integration of random EPSPs. J Neurosci 13: 334–350, 1993.[Abstract]
Stabler SE, Palmer AR, and Winter IM. Temporal and mean rate discharge patterns of single units in the dorsal cochlear nucleus of the anesthetized guinea pig. J Neurophysiol 76: 1667–1688, 1996.
Stancak A and Pfurtscheller G. Mu-rhythm changes in brisk and slow self-paced finger movements. Neuroreport 7: 1161–1164, 1996.[Web of Science][Medline]
Stein RB. A theoretical analysis of neuronal variability. Biophys J 5: 173–194, 1965.[Web of Science][Medline]
Stein RB and Matthews PB. Differences in variability of discharge frequency between primary and secondary muscle spindle afferent endings in the cat. Nature 208: 1217–1218, 1965.[CrossRef][Medline]
Stern EA, Kincaid AE, and Wilson CJ. Spontaneous subthreshold membrane potential fluctuations and action potential variability of rat corticostriatal and striatal neurons in vivo. J Neurophysiol 77: 1697–1715, 1997.
Svirskis G and Rinzel J. Influence of temporal correlation of synaptic input on the rate and variability of firing in neurons. Biophys J 79: 629–637, 2000.[Web of Science][Medline]
Tiesinga PHE. Chaos-induced modulation of reliability boosts output firing rate in downstream cortical areas. Phys Rev E Stat Nonlin Soft Matter Phys 69: 031912, 2004.[Medline]
Thomson AM. Presynaptic frequency- and pattern-dependent filtering. J Comput Neurosci 15: 159–202, 2003.[CrossRef][Web of Science][Medline]
Truccolo W, Eden UT, Fellows MR, Donoghue JP, and Brown EN. A point process framework for relating neural spiking activity to spiking history, neural ensemble, and extrinsic covariate effects. J Neurophysiol 93: 1074–1089, 2005.
Werner G and Mountcastle VB. The Variability of central neural activity in a sensory system, and its implications for the central reflection of sensory events. J Neurophysiol 26: 958–977, 1963.
Wetmore DZ and Baker SN. Post-spike distance-to-threshold trajectories of neurones in monkey motor cortex. J Physiol 555: 831–850, 2004.
Young ED, Robert JM, and Shofner WP. Regularity and latency of units in ventral cochlear nucleus: implications for unit classification and generation of response properties. J Neurophysiol 60: 1–29, 1988.
This article has been cited by other articles:
|
|
 |
|
  E. R. Williams and S. N. Baker Renshaw Cell Recurrent Inhibition Improves Physiological Tremor by Reducing Corticomuscular Coupling at 10 Hz J. Neurosci., May 20, 2009; 29(20): 6616 - 6624. [Abstract] [Full Text] [PDF]
|
|
|
|
 |
|
 E. R. Williams and S. N. Baker Circuits Generating Corticomuscular Coherence Investigated Using a Biophysically Based Computational Model. I. Descending Systems J Neurophysiol, January 1, 2009; 101(1): 31 - 41. [Abstract] [Full Text] [PDF]
|
|
|
|
|
|
D. S. Soteropoulos and S. N. Baker Quantifying Neural Coding of Event Timing J Neurophysiol, January 1, 2009; 101(1): 402 - 417. [Abstract] [Full Text] [PDF]
|
|
|
|
|
|
K. Nakamura, M. Matsumoto, and O. Hikosaka Reward-Dependent Modulation of Neuronal Activity in the Primate Dorsal Raphe Nucleus J. Neurosci., May 14, 2008; 28(20): 5331 - 5343. [Abstract] [Full Text] [PDF]
|
|
|
|
|
|
D. S. Soteropoulos and S. N. Baker Different Contributions of the Corpus Callosum and Cerebellum to Motor Coordination in Monkey J Neurophysiol, November 1, 2007; 98(5): 2962 - 2973. [Abstract] [Full Text] [PDF]
|
|
|
|
 |
|
 W. K. Mitchell, M. R. Baker, and S. N. Baker Muscle responses to transcranial stimulation in man depend on background oscillatory activity J. Physiol., September 1, 2007; 583(2): 567 - 579. [Abstract] [Full Text] [PDF]
|
|
|
|
|
|
C. L. Witham and S. N. Baker Network oscillations and intrinsic spiking rhythmicity do not covary in monkey sensorimotor areas J. Physiol., May 1, 2007; 580(3): 801 - 814. [Abstract] [Full Text] [PDF]
|
|
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
| HOME | HELP | FEEDBACK | SUBSCRIPTIONS | ARCHIVE | SEARCH | TABLE OF CONTENTS |
| Visit Other APS Journals Online |
