## Abstract

Irregularity of firing in spike trains has been associated with coding processes and information transfer or alternatively treated as noise. Previous studies of irregularity have mainly used the coefficient of variation (CV) of the interspike interval distribution. Proper estimation of CV requires a constant underlying firing rate, a condition that most experimental situations do not fulfill either within or across trials. Here we introduce a novel irregularity metric based on the ratio of adjacent intervals in the spike train. The new metric is not affected by firing rate and is very localized in time so that it can be used to examine the time course of irregularity relative to an alignment marker. We characterized properties of the new metric with simulated spike trains of known characteristics and then applied it to data recorded from 108 single neurons in the motor cortex of two monkeys during performance of a precision grip task. Fifty-six cells were antidromically identified as pyramidal tract neurons (PTNs). Sixty-one cells (30 PTNs) exhibited significant temporal modulation of their irregularity during task performance with the contralateral hand. The irregularity modulations generally differed in sign and latency from the modulations of firing rate. High irregularity tended to occur during the task phases requiring the most detailed control of movement, whereas neural firing became more regular during the steady hold phase. Such irregularity modulation could have important consequences for the response of downstream neurons and may provide insight into the nature of the cortical code.

## INTRODUCTION

When quantifying features of an experimentally recorded spike train, the firing rate is the most commonly used measure. However, other properties may also be of interest; one of these is the regularity of the discharge—how variable the interspike intervals are during periods of stationary rate. Irregularity of timing in spike trains has been has recently been the subject of considerable discussion in the literature in the context of rate coding versus temporal coding of the represented information (Shadlen and Newsome 1998; 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

### Definition and evaluation of irregularity metrics

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)/d*t*.

A metric that meets the above requirements and uses two adjacent intervals (*I*_{i} and *I*_{i}_{+1}) is (1) or alternatively for computational convenience (2) where | · | denotes absolute value, and log is the natural logarithm. The value of the metric applies over the entire time span of the two adjacent intervals; however, for computational convenience, we will treat it as a point event at the time of the middle spike. To avoid loss of data by binning, we first calculate the metric for each (overlapping) pair (i.e., intervals 1,2 then 2,3, etc) along each line of the raster, assigning each value to the time of the corresponding middle spike. After this is done for each trial of the raster, time bins can be imposed. In this way, even though a given interval pair straddles a time bin boundary, its irregularity metric value will be assigned to the appropriate bin, and all available interval pairs in the data will have been used. At this stage, we may examine the probability distribution function, average and SD of the metric within any time bin and across all trials. We will abbreviate the average metric as IR, where IR = mean(*m*). Characterization of the performance of IR as a measure of irregularity follows in results. Although we have performed numerical investigations of the properties of IR using intervals chosen from a gamma distribution, IR is independent of firing rate for interspike intervals that follow any arbitrary probability distribution function *f*(*I*). Formally, the expected value of IR for intervals distributed according to *f*(*I*) is the same as for intervals distributed according to *f*(*kI*).

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 (*I _{i}*) and local mean taken over

*n*adjacent intervals on either side (3) or alternatively for computational convenience (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 *L*_{V} (5) Like the metric of *Eq. 2,* this relies only on an interval pair and is thus well localized. It was not the aim of this paper to compare the performance of different metrics but rather to characterize one in detail. The remainder of this paper will therefore use the metric defined in *Eq. 2.*

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 *k _{i}* as the indices of those intervals that are smaller than the median, and

*j*as the indices of those that are larger than or equal to the median. The intervals were then reordered according to (6) The remaining (

_{i}*N*−

*M*) intervals were left in their original order. This method has the advantage that the interspike interval histogram is unchanged, whereas a negative serial correlation is introduced into the interval sequence. The magnitude of this correlation was varied by changing the fraction

*M*/

*N*from 0 to 50% in 10% steps.

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 O_{2}:N_{2}O) 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

We first characterize properties of the two-interval irregularity metric *m*; we will abbreviate its mean value over some appropriate set of intervals as IR [IR = mean(*m*); see methods]. Using intervals simulated using gamma distributions of constant rate and order, the dependence of IR and CV (coefficient of variation) on gamma order is shown in Fig. 1*A*. Exactly the same curves are obtained for any value of constant rate, demonstrating that IR and CV are both independent of constant firing rate value. Both IR and CV values change most rapidly between low gamma orders and change very little between gamma orders >15. The direct relation of IR to CV values is shown in Fig. 1*B* along with a diagonal line that would denote equality. Note that for gamma distributed intervals the CV is 1√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.

Analytically or by simulation we can obtain the probability density function (PDF) of the *m* metric. This is shown for two values of gamma order in Fig. 1*C*; the curve is one-sided because the *m* metric is defined by the absolute logarithm of the adjacent interval ratio. The highest probability occurs at an *m* value of 0, which corresponds to identical adjacent intervals. The probability rapidly decreases as adjacent intervals differ. Sample means from such a PDF are not normally distributed and have a long tail. This is shown in Fig. 1*D*, which shows the IR distribution over many samples of 10 values drawn from the order 4 PDF in Fig. 1*C*. Vertical lines indicated the mean and 2.5 and 97.5% points.

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 1*E* 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. 1*E* 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. 1*F*. The asymmetry of the PDF mean distribution in Fig. 1*D* underlies the difference between upper and lower confidence limits. Figure 1*F* 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. 1*A*), because it measures a ratio of two adjacent intervals, it is possible that its value could depend on the derivative of rate d(rate)/d*t*. 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 2*A* shows a step change in rate from 20 to 200 Hz. The IR measure (Fig. 2*B*) 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 2*C* illustrates a ramp rise in firing rate (from 20 to 400 Hz in 0.5 s). The time resolved IR (Fig. 2*D,* 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. 2*E*). Averaging many such step events produces Fig. 2*F,* 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. 2*E*), and the IR which is measured from pairs of intervals generated with stationary rate and the same gamma order (as in Fig. 1*A*). 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. 2*F* 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 (*I*_{2} in Fig. 2*E*) 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 (*I*_{1} and *I*_{3} in Fig. 2*E*) 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 *I _{i}* and

*I*

_{i}_{+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 2

*G*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 2

*H*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 3*A* 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.

To understand the source of the IR errors better, we introduce the signed metric *ms* = log(*I _{i}*/

*I*

_{i}_{+1}), omitting the absolute value in the definition of the normal irregularity metric

*m*(see methods). Define a new measure by

*S*= mean[log(

*I*/

_{i}*I*

_{i}_{+1})] over an appropriate set of interval pairs as in a time bin. The PDF of

*ms*is strongly modified by presence of a rate step. Figure 3

*B*compares this PDF for a gamma 8 process with and without a 5× rate step. Without the rate step the PDF is symmetric about zero, with the rate step, there is considerable asymmetry that is the source of the IR error. Note, however, that the distorted PDF values are on the negative side in this example. The PDF on the positive side is similar in shape to that with no step change. The mean

*S*will differ from zero due to the asymmetry caused by the rate step. Absolute values of such deviations from zero are shown (for a typical 100-ms bin and many trials) in Fig. 3

*C*. As in Fig. 3

*A*, after the initial increase up to a rate ratio of 5, there is a decrease for larger rate steps; this again is simply the dilution by normal intervals falling within the bin. Figure 3

*D*shows the upper 95% confidence limits for measure

*S*, for the situation where there is no change in rate, as a function of the number of metrics averaged. This indicates the limits within which we can expect

*S*to remain if there is no rate change.

The one-sided distortion of the PDF of *ms* shown in Fig. 3*B* 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 5× step condition shown in Fig. 3*B* 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. 3*E*. Note that the maximum span of the error is about 1/10 of that for the uncorrected IR values in Fig. 3*A*. The choice of threshold value on *S* to trigger the correction procedure is arbitrary, but Fig. 3*F* 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. 4*A*), 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. 4*D* 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. 4*C*) was quite different from the time course of the IR (Fig. 4*E*). 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. 4*E*). Figure 4*F* 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. 4*E,* where IR is low, are more narrowly distributed than those from the region with high IR (b).

Further examples of applying the IR metric to experimental data are given in Figs. 5 and 6 using the same layout as in Fig. 4. Figure 5 contrasts task execution with hands ipsi- and contralateral to the recording site. There was a profound rate modulation during contralateral performance and a much smaller but nevertheless clear change in firing rate when the ipsilateral hand was used (Fig. 5, *B* and *C*). The IR measure modulated substantially during contralateral performance but was constant during ipsilateral movement. The normalized interval distributions from selected regions of interest (Fig. 5*F*) show that this unit exhibited short interval bursting behavior; the changes in IR during contralateral performance predominantly reflected the altering proportions of bursts versus longer intervals. Note that the IR values of Fig. 5*E* lie mainly between 1 and 2.5. This cell's discharge was therefore at some times even more irregular than a Poisson process (expected IR value 1.39), reflecting the presence of bursts.

Contralateral data from two further neurons is shown in the same format in Fig. 6. Both neurons exhibited strong IR modulation with decreased IR during the hold phase. This similarity of IR modulation is of interest because these neurons showed opposite firing rate changes during the task. More generally, comparison of Figs. 4–6 shows that the IR modulation can be dissociated from the firing rate changes. During the hold phase, we show situations with high rate and low IR (Figs 4 and 6, *right*), low rate and high IR (Fig. 5, contralateral), and low rate with low IR (Fig. 6, *left*).

A population analysis of the available dataset is shown in Fig. 7. Figure 7*A* 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 7*B* 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. 7*C*). 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.

To investigate the variation of IR modulation around the population mean, a principal component analysis was carried out. This used the IR modulation during contralateral task performance of the significantly modulated cells. The first three principal components are shown in Fig. 7, *D—F;* together, these accounted for 75% of the variance in the IR modulation. Figure 7*G* shows the projection of the IR modulation data onto the second and third principal component axes. Both identified PT neurons (•) 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 7*G.* 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. 7*G*. The firing rate patterns from Fig. 7*I*, *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. 7*H*, *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. 7*G.* 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. 7*J*) and mean firing rate (Fig. 7*K*) 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

### Instantaneous irregularity measure

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. 2*G*). 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. 7*I*), 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

This work was funded by The Wellcome Trust to S. N. Baker and National Institutes of Health Grants MH-46428 and DC-01249 to G. L. Gerstein.

## Acknowledgments

The authors thank D. Soteropoulos, S. Rhodes, D. Wetmore, and R. Pyper for assistance in gathering the experimental data described in this report.

## Footnotes

The costs of publication of this article were defrayed in part by the payment of page charges. The article must therefore be hereby marked “

*advertisement*” in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.

- Copyright © 2006 by the American Physiological Society