HOME HELP FEEDBACK SUBSCRIPTIONS ARCHIVE SEARCH TABLE OF CONTENTS		QUICK SEARCH:   [advanced] Author: Keyword(s): Year:  Vol:  Page:

This Article Abstract Full Text (PDF) All Versions of this Article: 97/2/1221    most recent 00448.2006v1 Alert me when this article is cited Alert me if a correction is posted Citation Map Services Email this article to a friend Similar articles in this journal Similar articles in ISI Web of Science Similar articles in PubMed Alert me to new issues of the journal Download to citation manager Citing Articles Citing Articles via ISI Web of Science (1) Citing Articles via Google Scholar Google Scholar Articles by Chi, Z. Articles by Margoliash, D. Search for Related Content PubMed PubMed Citation Articles by Chi, Z. Articles by Margoliash, D.

Template-Based Spike Pattern Identification With Linear Convolution and Dynamic Time Warping

¹Department of Statistics, University of Connecticut, Storrs, Connecticut; ²Department of Statistics, Florida State University, Tallahassee, Florida; and ³Department of Organismal Biology and ⁴Anatomy, University of Chicago, Chicago, Illinois

Submitted 28 April 2006; accepted in final form 15 November 2006

	ABSTRACT

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION APPENDIX GRANTS ACKNOWLEDGMENTS REFERENCES

 
Pattern identification for spiking activity, which is central to neurophysiological analysis, is complicated by variability in spiking at multiple timescales. Incorporating likelihood tests on the variability at two timescales, we developed an approach to identifying segments from continuous neurophysiological recordings that match preselected spike "templates." At smaller timescales, each component of the preselected pattern is represented by a linear filter. Local scores to measure the similarities between short data segments and the pattern components are computed as filter responses. At larger timescales, overall scores to measure the similarities between relatively long data segments and the entire pattern are computed by dynamic time warping, which combines the local similarity scores associated with the pattern components, optimizing over a range of intercomponent time intervals. Occurrences of the pattern are identified by local peaks in the overall similarity scores. This approach is developed for point process representations and binary representations of spiking activity, both deriving from a single underlying statistical model. Point process representations are suitable for highly reliable single-unit responses, whereas binary representations are preferred for more variable single-unit responses and multiunit responses. Testing with single units recorded from individual electrodes within the robust nucleus of the arcopallium of zebra finches and with recordings from an array placed within the motor cortex of macaque monkeys demonstrates that the approach can identify occurrences of specified patterns with good time precision in a broad range of neurophysiological data.

	INTRODUCTION

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION APPENDIX GRANTS ACKNOWLEDGMENTS REFERENCES

 
Neurons exhibit patterned activity at multiple scales of time (Reich et al. 2000, 2001; Zugaro et al. 2004). Analysis of such data, often obtained as extracellular recordings, is central to the goals of neurophysiology. A subset of this general problem is identification of spiking activity patterns in one record that match previously identified patterns in another record. For example, the analysis of sleep mechanisms of learning and memory has in part involved comparison of ongoing discharge during sleep with discharges recorded during wakeful behavior (Buzsáki 1989). Based on the identified patterns, some inferences can be made on the information processing of the brain during sleep (Buzsáki 1998; Dave and Margoliash 2000; Hahnloser et al. 2006; Louie and Wilson 2001, 2002; Nádasdy 2000; Skaggs and McNaughton 1996; Wilson and McNaughton 1994). This is a difficult problem because of the variability of neuronal activity in sleep relative to wakeful behavior and because there is no time referent during sleep. As another example, from a neural decoding perspective it is of interest to estimate the times of certain landmark events during a course of behavior based on the associated spiking activity (Chi and Margoliash 2001). The estimation of behavioral timing can be useful when it is combined with procedures designed to estimate spatial trajectories of motor output (Georgopoulos et al. 1986; Hatsopoulos et al. 2004; Smith and Brown 2003; Wu et al. 2004, 2006); however, it appears that this type of estimation has not received much attention in the literature. We show herein that pattern identification for spiking activity can provide a good estimate for behavioral timing.

Template matching has been the basis of many pattern-identification procedures (Abeles et al. 1988, 1993; Chi et al. 2003; Dave and Margoliash 2000; Louie and Wilson 2001; Nádasdy et al. 1999). In these procedures, a template is constructed based on exemplar spike trains that exhibit a certain pattern of interest. During data analysis, each segment of the data is scored in terms of its similarity to the template. Segments with high similarity scores are interpreted as exhibiting the pattern of interest and are thus extracted for further analysis.

In most of the current procedures, the template is treated as nearly time locked to stimuli or behavioral events. Except for some local variability, such as random perturbation in spike timing, global variability is expressed as uniform timescaling, which multiplies all the interspike intervals (ISIs) by a single factor (Chi et al. 2003; Louie and Wilson 2001; Nádasdy et al. 1999). However, in some systems spiking activity can exhibit "time warping," i.e., nonuniform time compression or expansion on a relatively large timescale (Dave and Margoliash 2000; Nádasdy et al. 1999). As far as we know, the procedure introduced by Victor and Purpura (1996) is the only one that explicitly addresses time warping. The procedure was developed to analyze the precision of temporal coding during specific short time intervals of sensory stimulation. In the procedure, the ISIs are assumed to vary independently and similarity scores are computed with dynamic time warping (DTW). Such scores are known as "text edit distances" in the field of computer science (Wagner and Fischer 1974). Similar procedures based on text edit distances are extremely useful in bioinformatics (Aach and Church 2001; Bishop and Thompson 1986; Needleman and Wunsch 1970; Waterman et al. 1987). Sequences of discrete behaviors such as syllables in birdsongs are also occasionally analyzed using the text edit distance (Margoliash et al. 1991; Sankoff and Kruskal 1983). DTW is also commonly used in speech recognition (Huang and Gray 1988; O’Shaughnessy 1999; Rabiner and Juang 1993), birdsong recognition (Anderson et al. 1996; Ito and Mori 1999; Kogan and Margoliash 1998), and other applications.

The variability of spiking activity at local versus relatively large timescales can be roughly described as follows (Fig. 1A). Think of the occurrences of a pattern of interest as "texts." The "codewords" of the texts are short intervals with spikes. The definition of codeword should be problem dependent. For different systems, a codeword may contain a single spike, multiple spikes, or a high-frequency burst. Within the codewords, the detectable variability can be reasonably characterized as random local perturbation in spike timing and spike count. Between the codewords, on the other hand, the gaps often have quite visible nonuniform variability, resulting in time warping. The gaps need not be activity free. Instead, they can contain "noise" spikes that do not repeat in different occurrences of the pattern.

View larger version (12K): [in this window] [in a new window]   FIG. 1. A schematic of template matching. A, top: a sketch of a template registered in [0, D0]. Template may be a single-unit (SU) spike train, an array of firing rates of multiple single units, or other suitable representations of spiking activity. Shaded blocks B1, ..., B4 represent "codewords" of the template that are almost "rigid" and can be perturbed only locally, whereas the blank blocks I1, ..., I5 represent intervals that can be compressed or expanded in time. Bottom: matching the template with a data segment at time x. Arrows define how different parts of the template and those of the segment are matched. Onset of the template is mapped to the onset x of the data segment. By compressing or expanding Ii, the codewords are aligned and compared with the data in the shaded blocks indicated by the arrows. Overall similarity score between the template and the data segment is the sum of the similarity scores between the codewords and the corresponding data intervals, minus a weighted sum of the modifications on the intervals Ii. Depending on problems, the spiking activity in the intervals Ii of the data are either ignored or treated as noise and included in the overall similarity score as a minus term. B: template matching for single-unit spike trains, as for the nucleus robustus arcopallium (RA) data. Top: each Bi is a burst in the template and Ii is an interburst interval (IBI). Bottom: bursts in the data segment are matched with the template bursts. Spike s in the data segment does not correspond to any burst in the template and is counted as a "noise" spike. Amount of noise in Ii is evaluated by the number of spikes in them.

The local comparison should be based on specific representations of spiking activity. We shall consider two types of representations. The point process representation registers spike times continuously. As a result, the neural activity is characterized as processes of points in time. The binary representation divides the neural activity into small time bins, and assigns 1 to a bin if there is at least one spike in it and 0 otherwise. As a result, the neural activity is characterized as sequences of 0 and 1 s. Although the local comparison algorithms for these two representations look quite different, they derive from nearly identical arguments based on likelihood tests.

By using DTW, the template-based approach automatically "parses" each identified data segment into codewords and gaps in between, thus allowing a more detailed interpretation of the identification results (Fig. 1A). Such parsing can yield additional insight. If the spiking activity pattern is preselected in association with a certain behavior, and each of its codewords is reliably associated with a landmark event during the behavior, then the parsing provides a way to infer from the neural activity when the landmark events occur within the overall behavior. Presently, the dominant procedures for inference of the time course of behavior are based on regression or state-space models (Brockwell et al. 2004; Carmena et al. 2003; Serruya et al. 2002; Smith and Brown 2003; Wu et al. 2004, 2006). They are designed to estimate spatial trajectories of motor output with minimized average error over time. The estimation of the timing of events during motor output can therefore provide a useful new perspective to the decoding of neuronal activity.

In what follows, under METHODS, we describe the identification procedures and present pseudocode for their implementations. The procedure for the point process representation was implemented in C++, and that for the binary representation in MATLAB (The MathWorks, Natick, MA). Under RESULTS, we illustrate the utility of the procedures with two brief neurobiological examples. The procedure for the point process representation is applied to a single-unit (SU) recording collected during sleep from the nucleus robustus arcopallium (RA) of zebra finches (Dave and Margoliash 2000). The procedure for the binary representation is applied to multielectrode recordings from the primary motor cortex (MI) of a behaving macaque monkey during a target pursuit task (Hatsopoulos et al. 2004; Serruya et al. 2002). In this example, the times when a target was reached and the following target appeared in each successful trial can be estimated with error <0.2 s. This is somewhat surprising in view of the relatively high level of spike-timing variability in MI. For each procedure, we conduct a sensitivity analysis, demonstrating that the procedures are robust.

	METHODS

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION APPENDIX GRANTS ACKNOWLEDGMENTS REFERENCES

 
The neuronal recordings analyzed in RESULTS were obtained from adult male zebra finches (Taeniopygia guttata) and Macaca mulatta monkeys. The experiment procedures were conducted under protocols approved by the University of Chicago IACUC for birds (Margoliash) and monkeys (Hatsopoulos).

Because of the nature of the recording technology and cortical activity, the levels of accuracy of these two recordings were very different. First, in terms of signal-to-noise ratio (SNR), the accuracy of separating action potentials from background was much higher for the recordings from RA. The SNRs of spikes from SU recordings in the RA of zebra finches are generally >10, whereas those from the array recordings in macaque monkeys can be as low as 3. The difference in SNRs result from the different experimental designs and thus different recording technologies are used: movable but small numbers (one to four) of electrodes in the bird versus a fixed array of large numbers (>100) of electrodes in the monkey. Second, in terms of neurophysiological characteristics, the spike-timing variability of zebra finch RA units is much lower as well (Yu and Margoliash 1996; cf. Hatsopoulos et al. 2004). The differences are likely to influence the performance of pattern identification. The sampling rate of recording is 20 kHz for the RA data and 30 kHz for the MI data.

Part one: point process representation

The procedure described in this part is suitable when 1) the pattern of interest consists of well-defined spike bursts; 2) across occurrences of the pattern, the bursts exhibit low variability; and 3) the data being analyzed constitute a single spike train, consisting of time registry of spikes accurately estimated from the SU recording. Under these conditions, an exemplar spike train expressing the pattern may be used as a template. Under RESULTS, more detail will be given on how templates are selected for SU activity in RA.

The procedure directly compares the spike times. Figure 1B illustrates how to compare a template with data registered around a time location. In general, the template need not start or end with a spike. Together with its bursts, the intervals at the beginning and at the end of the template and those between the bursts are all parts of the template and will be referred to as interburst intervals (IBIs). The bursts are used as the codewords of the template pattern. The IBIs can be thought of as independent "springs," and the bursts as relatively incompressible "blocks" interconnected by the springs. By compressing or stretching the IBIs, the template bursts are shifted accordingly and compared with the data registered in the corresponding time intervals. Given a set of modified IBIs, the score of similarity between the warped template and the data segment aligned with it is a sum of the similarity scores associated with the bursts, minus a weighted sum of the amount of changes in the IBIs and the number of spikes in the data segment that do not match any template spikes. Optimizing the weighted sum over a range of modified IBIs then yields the optimal score at the time location. We next consider appropriate similarity scores and their computation.

TEMPLATE. A template spike train has to be registered in an interval [0, D] (Fig. 1B). If the template has n bursts, then in general it has n + 1 IBIs, including the interval between its onset (t = 0) and first spike, and the interval between its last spike and offset (t = D). A subset of template spikes is a burst if all its ISIs are less than a specified threshold. Fix parameter > 0, which controls the precision for comparing spike times. Let H_n+1 = D, T₀ = 0, and for 1 i n

The ith burst interval is defined as [H_i, T_i] with duration B_i = T_i – H_i; the ith IBI is defined as [T_i–1, H_i] with duration I_i = H_i – T_i–1 (H stands for "head"; T for "tail").

LOCAL SCORES FOR TEMPLATE BURSTS. First, we fix a kernel function K(x). This will be used to measure the temporal discrepancy between a data spike and a template spike. Some commonly used kernel functions defined on [–1, 1] are constant 1 ("square" kernel), 1 – |x| ("triangular"), 1 – x² ("Epanechnikov"), and (1 – x²)² ("biweight"). To evaluate the similarity between the ith template burst and a data burst s₁, s₂,... in [x, x + B_i], the idea is to translate both to [0, B_i] so they can be compared directly. If the ith template burst consists of spikes t₁, t₂,..., then the similarity score between the two bursts is

(1)

with 0 the maximum absolute negative contribution a nonmatching data spike can make to the local score. From Chi et al. (2003), F_i(x) can be computed as a linear convolution of the entire data spike train with _i.

GLOBAL SCORE FOR THE ENTIRE TEMPLATE. Denote by N(a, b) the number of data spikes between a and b. For the ith IBI, let G_i(x) 0 be a cost function for its warping. Suppose the warped template has its IBIs changed by v₁, ..., v_n+1. Denote

Then the total similarity score between the warped template and the data at x is

(2)

Equation 2 is derived in the appendix and can be explained by Fig. 1B. By changing the IBIs, the ith template burst is aligned to the data in [x + H_i + V₁ⁱ, x + T_i + V₁ⁱ], yielding the first term on the right-hand side of Eq. 2. The second term penalizes the changes on the template IBIs. In the warped template, the IBIs are [x + T_i–1 + V₁^i–1, x + H_i + V₁ⁱ], 1 i n + 1. Data spikes in these intervals are counted as "noise" and reduce the similarity. This gives rise to the third term in Eq. 3.

By the probabilistic interpretation given in the APPENDIX, the overall similarity score assigned to the time location x is defined as

(3)

In general, there can be multiple solutions to the optimization in Eq. 3, each one being a set of changes in the IBIs. Let ₁(x), ..., _n+1(x) be a solution. By Fig. 1B, we can establish a one-to-one correspondence between the intervals I₁, ..., I_n+1 in the template and I₁, ..., I_n+1 in the data segment and identify the latter as a set of optimal matching IBIs. Then we can identify the intervals B₁, ..., B_n between I₁, ..., I_n+1 as a set of optimal matching burst intervals. The interlacing I₁, B₁, I₂, ..., B_n, I_n+1 form a partition of the data segment between x and the right endpoint of I_n+1, which is an optimal match to the template that one can possibly find at time x under the similarity score. The optimal matching IBIs I_i and burst intervals B_i are computed as follows

(4)

Because of the one-to-one correspondence between their subintervals, we can inspect in more detail how similar the optimal matching data segment and the template are.

DTW FOR THE GLOBAL COMPARISON. To compute L(x) and ₁(x), ..., _n+1(x), let

Clearly, l(x, v₁, ..., v_n+1) = l₁(x, v₁, ..., v_n+1) and thus L(x) = L₁(x). Because

it follows that L_k(x) = max_v [M_k(x, v) + L_k+1(x + v)]. Therefore L(x) and ₁(x), ..., _n+1(x) can be computed by DTW. Table 1 summarizes the procedure.

1. )L(x) , where is a given threshold;
   
2. )L(x) = max {L(y) :|y – x| r}, with r the radius for local comparison; and
   
3. )L(x) > L(y) for at least one y with |y – x| r. This is to make sure that L is not a constant in [x – r, x + r].
   

A reasonable choice for r is D, the duration of the template. When some of the extracted segments overlap in time, a disjoint subset of the segments may be selected.

The selection of the threshold value in general is a difficult problem. One important reason is that the stochastic properties of the data can be difficult to grasp. In RESULTS we will describe an ad hoc method for the RA activity that tries to simulate the noise in the data. A different ad hoc method was used for the data collected from the MI of monkey, which will also be described in RESULTS.

ADAPTIVE SELECTION OF PARAMETER VALUES. The parameters and can be set based on estimation. For , we first compute the mean value d of the template ISIs, which can be regarded as the minimum resolution to distinguish template spikes. If the squarekernel is used in Eq. 1 to compute the local scores, then set = d/2. To use a different kernel K, has to be adjusted. The reason is that except for x = 0, K(x) is strictly less than the square kernel. Without adjusting , the kernel puts less total weight to data spikes that match the template spikes, resulting in artificially lower similarity scores. We set

(5)

where c_K > 1 such that the area under K(x/c_K) is 2. The value of c_K is 2 for the triangular kernel, 1.5 for the Epanechnikov kernel, and 1.875 for the biweight kernel.

From Eq. A1 in the APPENDIX, log (₀/)/log ('/₀), where ' is the firing rate within template bursts, is the firing rate within template IBI, and ₀ is the mean firing rate in nonmatching spike trains. Because for a Poisson process the density is the reciprocal of the mean ISI, we set the value of based on the approximation

(6)

where d', d, and d₀ are the mean ISIs in template bursts, template IBIs, and the entire data spike train, respectively. In estimating d, if a template IBI has no spike, then its duration substitutes for a sample ISI. Whereas d' and d are fixed once templates are selected, d₀ is dependent on data. Therefore the selection of is adaptive.

NUMERICAL ISSUES. In data analysis, L(x) is computed on a discrete time grid jd, with j = 0, ±1,..., and d > 0 a fixed step size. Also, the warping of each IBI has to be bounded. Table 2 contains a version of the algorithm that takes into account these constraints. It requires that F_i(x₀ + jd) be computed and stored for a given x₀. We apply the subroutine below to compute F_i(x₀ + jd). Let a denote the smallest integer no less than a, and a the largest integer no greater than a.

1. )Let F_i(x₀ + jd) = 0 for all j.
   
2. )For each data spike s and j = j₀, ..., j₁, increase F_i(x₀ + jd) by _i(s – x₀ – jd), where j₀ = (s – x₀ – )/d and j₁ = (s – x₀)/d.
   
3. )Return all F_i(x₀ + jd).
   

The procedure described in this part is suitable when 1) the spiking activity is recorded from multiple single units and exhibits significant amount of variability in spike timing across trials and/or 2) the neural signals have lower SNRs compared with SU recordings; typically, the spike trains are obtained from multielectrode recordings by spike sorting. Under these conditions, the spike train of a single unit typically cannot represent the variability in a pattern of interest and therefore it alone cannot be directly used as a template for pattern identification. Instead, we shall construct a template as a linear filter that is made up of log-likelihood ratios estimated from sample spike trains of all the recorded single units. The statistical model involved in the construction is an integral part of the pattern identification procedure and is described below.

BASIC SETUP. Suppose the spiking activity pattern we are interested in is associated with a sequence of n landmark events during the course of a certain type of behavior. For example, if a behavioral task requires a monkey to reach n targets in sequence, then for i = 1, ..., n, the ith event may be defined as the cursor reaching the ith target (see more detail in RESULTS). The n – 1 intervals between the times of consecutive events during an occurrence of the event sequence will be referred to as interevent intervals (IEIs).

We suppose that the multielectrode spike trains are recorded from the same set of C units and each data spike has been associated with one of the units by spike sorting. Fix parameter > 0 as the duration of time bin. For c = 1, ..., C, the binary representation for thespiking activity of the cth unit is a binary sequence x_c(1), x_c(2), x_c(3), ..., such that

When is small enough, almost all time bins can have at most one spike, and the binary representation induces little loss of information on spike count.

To start, one needs to collect a sample of multielectrode spike trains associated with multiple occurrences of the event sequence. For each sample spike train, if T₁, ..., T_n are the associated event times, then the segments of the spike train in surrounding intervals [T_i – A, T_i + B] are registered as codewords, where A, B 0 are fixed integers. As in the case of point process representation, the codewords cannot be warped in time.

STATISTICAL MODEL. The procedure is based on the following assumptions.

• Conditional on the occurrence of an event at time t, the spiking activities of the units in the time interval [t – A, t + B] are independent, such that the activity of each unit forms a Poisson process.
  
• Conditional on the occurrence of the entire event sequence, the joint spiking activities of the units associated with the individual events and the IEIs are all independent of each other.
  

For i = 1, ..., n, c = 1, ..., C, and j = –A, ..., B, let p_i,c(j) be the probability that the cth unit generates spikes in the jth time bin away from the ith event. The probability can be estimated using a peristimulus time histogram method. That is, for the ith event and cth unit, we first align the sample spike trains of the unit at the event, then estimate p_i,c(j) as the fraction of those sample spike trains that have spikes in the jth time bin relative to the time of the event.

For i = 1, ..., n – 1, we use q(t, _i) to model the probability density of the duration of the ith IEI. The parameter _i can be estimated based on the sample as well.

LOCAL SCORES FOR SPIKING ACTIVITY PATTERNS AROUND EVENTS. To identify the occurrences of the whole event sequence, the first step is to evaluate the likelihood of the occurrence of a single event in each time bin. The derivation of the likelihood is given in the APPENDIX. Here we consider the computational procedure only. For the ith event and the cth unit, we introduce _i,c(j), such that for j = –A, ..., B, it is evaluated as

It follows that

(7)

where _* denotes the mathematical convolution and _i,c*x_c(t) is the tth entry of the convolution between _i,c and x_c. The derivation of the formula is given in the APPENDIX. We therefore use

as the similarity score between the joint spiking activity starting in the tth time bin and the pattern associated with the ith event. The sequence _i,c acts as a discrete linear filter. By convolving _i,c with the entire binary sequence x_c, the score is obtained for all time bins.

GLOBAL SCORE FOR THE EVENT SEQUENCE. For the ith IBI, i = 1, ..., n – 1, and t 1, let G_i(t) = –log q(t, _i). Then for any time bin t and I₁, ..., I_n–1, the logarithm of the probability that the tth time bin contains the onset of an occurrence of the event sequence with IEIs I₁, ..., I_n–1 can be expressed as l(t, I₁, ..., I_n–1) plus a constant, where

(8)

We therefore define the score assigned to the time bin t as

(9)

Up to an additive constant, L(t) is the maximum likelihood that an occurrence of the event sequence starts in the tth time bin.

Denote by Î₁(t), ..., Î_n–1(t) the maximizers associated with L(t). As in the case of the point process representation, L(t), Î₁(t), ..., Î_n–1(t) can be computed using DTW. The algorithm is sketched in Table 3.

The actual times of the events are estimated as t_i. In practice, because the score function L(t) can be quite noisy, it may need to be smoothed to extract meaningful local maximum points.

To summarize, the numerical procedure for the binary representation starts with estimating p_i,c(j) and q(x, _i). Next the procedure constructs filters _i,c for i = 1, ..., n, c = 1, ..., C, and also constructs G_i(x) for i = 1, ..., n – 1. The scores and potential occurrences of the event sequence are then computed by the algorithm in Table 3.

	RESULTS

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION APPENDIX GRANTS ACKNOWLEDGMENTS REFERENCES

 
Application to RA of zebra finch

We applied the procedure for point process representations to the SU spiking activity in RA of adult male zebra finch. The data were high signal-to-noise ratio recordings from chronically implanted movable electrodes in RA, collected in a previous study (Dave and Margoliash 2000). Male birds directed their songs toward females housed in adjacent half-cages.

For each RA unit, the goal was to identify segments in its spontaneous spiking activity during sleep that exhibited temporal patterns similar to its premotor activities associated with daytime singing. Such segments have been interpreted as replay of premotor activity during sleep that may be involved in learning and maintaining song (Margoliash 2002). The female-directed songs of adult male zebra finches are highly stereotyped intense behaviors (Sossinka and Böhner 1980), acquired through learning (Brenowitz et al. 1997), and are dynamically maintained in adults by auditory feedback (Nordeen and Nordeen 1992). The songs constitute sequences of "syllables" repeated in a fixed pattern to form one or more repeated "motifs." Premotor spiking activity in RA associated with song renditions are highly stereotyped and exhibit precise temporal patterning within bursts, with the SD in spike timing as low as 0.2 ms (Chi and Margoliash 2001; Yu and Margoliash 1996). Therefore the point process representation is preferred to the binary representation for RA activity.

For these data, the spontaneous activity of a unit during sleep can be probed with one or more templates. Once a signature motif of the bird was identified, any spike train of the unit associated with a rendition of the motif could serve as a template (Fig. 2, A and B). To partition a template into bursts and IBIs, we aligned the spike trains associated with all the identified renditions of the motif by their onset times (Fig. 2C). The raster plot revealed common temporal structures of the spike trains that had a many-to-many correspondence to the acoustic structures of the motif. Based on the structures, we partitioned the spike trains into bursts and IBIs by the criteria that adjacent bursts be at least 20 ms apart. Within the bursts of these spike trains, the mean ISI was d' = 3.06 ± 2.68 ms, and within their IBIs the mean ISI was d = 73.02 ± 55.97 ms. The same partition was obtained when we applied a more refined alignment procedure which minimized the total L₁ distances of the spike trains (Chi and Margoliash 2001).

View larger version (46K): [in this window] [in a new window]   FIG. 2. A: spectrograph and waveform of a rendition of zebra finch song motif, which has a duration of 660 ms. B: a single-unit neuronal trace associated with the motif rendition. Trace has the same duration as the motif rendition, but starts 40 ms ahead. Lead time is taken out in the plot for the 2 to be aligned. C: raster plot of spike trains associated with different renditions of the same motif, aligned by the onset times of the renditions. Neuronal trace in B, in the 3rd row of the plot, is partitioned into 6 burst intervals marked by horizontal bars and 7 IBIs based on the raster plot.

To implement the procedure in Table 2, we used the biweight kernel to compute local similarity scores by Eq. 1. The parameter to control the temporal precision for spike matching was = d'/2 x 1.875 = 3.06/2 x 1.875 = 2.869 ms by Eq. 5. To calculate the global similarity scores by Eq. 2, we set G_i(x) = 0 and fixed the maximum proportion of warp at 1/5 for each template IBI. The score function was computed on a grid with step size 0.5 ms.

During the spontaneous activity of the unit in sleep, the mean ISI was d₀ = 49.05 ± 20.85 ms. By Eq. 6, the penalty on each nonmatching spike was = 0.1434. We then calculated the similarity scores across the data to identify segments that exhibited similar patterns to the template.

The results from one session of the data are shown in Fig. 3; plot A displays the similarity scores across. To identify matching segments, we set the threshold for similarity score = N/3, where N = 41 was the total number of spikes in the template. Four above-threshold local peaks were identified in this session, with scores 16.6, 19.2, 15.9, and 15.4, respectively. The second and fourth peaks are marked by x and z in Fig. 3A. By Eq. 4, we extracted the optimal matching segments at the peak locations and partitioned them into burst intervals and IBIs that were mapped one to one to those in the template (Fig. 3, B and C). The identified segments exhibited a high degree of similarity to the template and also showed nonuniform time warping of the IBIs.

View larger version (23K): [in this window] [in a new window]   FIG. 3. Detection for RA activity. A: similarity scores over 5 min of the SU spontaneous activity in RA during sleep, using the spike train in Fig. 2B as the template. Three local peaks in the similarity scores are marked as x, y, and z. B, top row: enlarged view of the scores in an interval of 11 s around x. 2nd row: RA activity in the same interval. Horizontal bar indicates a segment of RA activity that may have a pattern similar to that of the template. Onset of the segment is the location of the peak. 3rd row: an enlarged view of the segment of neuronal trace. Bursts in the segment are labeled with the indices of the corresponding optimal matching template bursts. A match to the 3rd template burst is missing. 4th row: same spike train in Fig. 2B shown for comparison. C: neuronal trace in the interval marked by z. Spike between the bursts labeled 4 and 5 is counted as a "noise" spike. D: neuronal trace in the interval marked by y. Relatively low score at y is indicative of the low level of similarity between the template and the activity in the interval.

To see whether DTW was necessary in this example, we collected all the identified optimal matching segments from the data. There were 15 of them, with mean score 16.6 ± 1.9. To see how much warping of the template was involved to match these segments, we analyzed the two longest IBIs in the template, which were the fourth and fifth IBIs. The time compressions/expansions of the two IBIs, –4.2 ± 21.5 and –2.5 ± 9.6 ms per segment, showed no significant correlation with each other (r = 0.34, P = 0.22, 95% CI = [–0.21, 0.73], t-test). For instance, to match the segment at x in Fig. 3A, the two IBIs were compressed by 33.5 ms and expanded by 6.5 ms, respectively. This suggests that nonuniform, wide-ranging time warping was present in the data, supporting the use of DTW for the pattern identification.

Another question is whether the probability distribution of the time warping in the premotor activity could be used in the pattern identification for the sleep data. The answer seems to be negative. There is evidence that the variability of IBI in the spontaneous activity during sleep differed from that during wakeful behavior. In the preceding data, the variance of the fourth IBI in the identified matching segments was significantly larger than the one in the exemplar spike trains associated with motif renditions. The ratio of the two was r = 18.61 (P < 10^–5, 95% CI = [6.25, 55.44], F-test). The variances of the fifth IBI on the other hand did not exhibit significant difference (r = 2.20, P = 0.15, 95% CI = [0.74, 6.54]). The change in the statistical properties of neuronal activity across different behavioral states poses difficulty for the selection of parameter values, as we tried to use the templates obtained in one state to analyze the activity in another.

Sensitivity to temporal precision of spike matching

We had used = 1.875, where = 1.53 ms was the half mean ISI within the bursts of the exemplar spike trains associated with the motif renditions, and the coefficient 1.875 "standardizes" the biweight kernel. To examine how sensitive the pattern identification was to , we compared the results to those obtained with different preselected values of = 1.875 while keeping the other parameters unchanged. For larger , the temporal precision required for template matching was lower and more segments with above-threshold similarity scores were expected. The question is how sensitive the number of identified segments was to the change in . Denote by n the number of segments that were identified using a preselected value of but missed by the original procedure that used = 1.53 ms and n' the number of those the other way around

Recall that the original procedure identified a total of 15 segments. Therefore the summary shows that the number of identified segments was reasonably stable for around the value selected based on the exemplar spike trains. Except for = 0.5, all the extra segments identified using the preselected values of were not completely overlooked by the original procedure. They actually showed up as local peaks with relatively high scores (>10). Therefore the increasing trend in n could be offset by lowering the threshold in the original procedure. The extra segment identified with = 0.5 was also not completely overlooked by the original procedure. It had roughly 50% overlap in time with a segment identified by the procedure.

Sensitivity to the penalty on mismatching spikes

In the procedure, the penalty = 0.1434 was selected adaptively, based on the firing rates of the premotor activities associated with motif renditions and the spontaneous activity during sleep. We compared the results to those obtained with different fixed values of while keeping the other parameters unchanged. Let n be the number of segments identified by using a preselected value of but missed by using the adaptively selected value and n' the number of those the other way around:

Recall that with the adaptively selected = 0.1434, 15 optimal matching segments were identified. The summary shows that for values of around the adaptively selected value, the identification results were quite similar.

Sensitivity to kernel

We evaluated how sensitive the pattern identification procedure was to the kernel. We implemented the procedure using the square, triangular, and Epanechnikov kernels, respectively. By Eq. 5, had to be adjusted accordingly. Denote by n the number of segments that were identified using a different kernel but missed by using the biweight kernel and n' the number of those the other way around:

Recall that with the biweight kernel, 15 optimal matching segments were identified. Therefore the summary shows that the procedure was reasonably stable by using different kernels. The differences in the number of identified segments arose from the small changes in similarity scores. For example, all four extra segments identified using the square kernel actually generated local peaks in the score function obtained with the biweight kernel. The scores were just slightly below the threshold = 41/3, with differences 0.54, 1.20, 1.47, and 0.50, respectively.

Simulations for threshold selection and performance evaluation

For the pattern identification procedure, the issues of threshold selection and sensitivity to noise are complicated by the fact that the statistical properties of the sleep data are likely to be different in some unknown but important aspects from the premotor activity that generated the templates. The simulation method is based on the idea that by testing the pattern identification procedure on artificial data that incorporate certain statistical properties of the sleep data, some reasonable inference might obtain.

We generated simulated data using the exemplar spike trains associated with motif renditions (Fig. 2A) because they had unambiguous interpretation. The pattern identification procedure was then applied to the simulated data to yield score functions. Because the peak values in these score functions were associated with "true targets" (i.e., spike trains known to exhibit the pattern of interest), we could use their distributions to draw inferences on the ability of the procedure to find true targets corrupted by noise.

Given an exemplar spike train T registered in [0, D], a simulated spike train was generated as follows.

• Each spike in T was randomly deleted with probability q and each remaining spike was randomly shifted by s N(0, ).
  
• The modified T was embedded into a larger interval I = [0, D + 2B] and a sample from the Poisson process with density was added to I as "noise" spikes.
  
• Finally, some of the spikes were removed from the resulting spike train so that all remaining ISIs were at least .
  

Because DTW automatically accounts for time warping of IBIs, we did not randomly modify the IBIs. In all the simulations, = 1 ms, which is a plausible lower bound for the ISI in RA. We also fixed B = 500 ms in all the simulations. The value of was selected based on the average firing rate of the spontaneous activity during sleep, which was estimated at 21 Hz. We set = 20 Hz. The parameters q and are unobservable from the spontaneous activity during sleep. We set around or above the half-mean (=1.53 ms) of the ISI within the bursts of the exemplar spike trains, so that there was a good chance for a spike in the exemplars to have a large "jitter" relatively to the ISIs. We conducted four simulations, with q = 1/4, 1/3 and = 1.5, 2 ms, respectively.

For each of the exemplar spike trains, 100 simulated spike trains were generated. The raster plot of a total of 1500 simulated spike trains with = 1.5 ms and q = 1/3 is shown in Fig. 4. Despite the relatively high chance of spike deletion in the true targets, the simulated spike trains maintained the overall temporal structure of the exemplar spike trains. However, not all of them exhibited enough similarity to the exemplars. A small percentage a of the simulated spike trains either failed to yield similarity scores above = N/4, with N the number of spikes in the template, or failed to have peak locations in score function within 50 ms from the onsets of the true targets. Excluding these spike trains, denote by S the peak value in the score function of any remaining spike train. The distribution of S was close to a Normal distribution. As or q increased, a increased and the distribution of S shifted toward 0.

A threshold could be selected based on the distribution of S. A sensible choice would be µ_S – 2_S. For all the simulations, this value was a little less than the actually used threshold value N/3. Therefore a few more segments would have been identified, and the sum of a and the tail probability under N(µ_S, _S) (2.5%) could serve as a nominal false negative rate.

View larger version (76K): [in this window] [in a new window]   FIG. 4. Raster plot of simulated spike trains with "noise" added to randomly modified copies of the exemplar spike trains in Fig. 2A. Horizon bar represents the time interval of 660 ms where the random copies are inserted. Each exemplar spike train generated 100 random spike trains that appear in the same block in the raster plot.

Application to MI of macaque monkey

We applied the procedure for the binary representation of neuronal activity to recordings obtained from the arm area of MI in the left hemisphere of a trained adult macaque monkey. The data were collected using a chronically implanted silicon microelectrode array composed of 100 electrodes (1.0-mm electrode length; 400-µm interelectrode separation; Cyberkinetics Neurotechnology System, Foxboro, MA). Based on histological evidence from previous implants that we have performed, the electrode tips are likely to be in lower portions of layer 3 or in layer 5. To acquire extracellular action potentials, signals were amplified (gain = 5,000), band-pass filtered (250 Hz to 7.5 kHz) and sampled at 30 kHz per channel. Only waveforms that crossed a threshold were stored and spike-sorted using Off-line Sorter (Plexon, Dallas, TX). The monkey was trained to repeatedly perform a sequence of movements (see following text) to receive a juice reward. From a decoding perspective, it is of interest to have a means to estimate the times of certain landmark events in the behavior during each trial based on the associated neuronal activity. These estimates were compared with the actual event times recorded during each trial.

The monkey was operantly trained to perform a target pursuit task by moving a cursor to a target with its right arm. The cursor and the target were displayed on a horizontal projection surface in front of the monkey. The monkey’s arm rested on cushioned arm troughs secured to links of a two-joint robotic arm (KINARM system; Scott 1999) underneath the projection surface. The shoulder joint was abducted 90° such that shoulder and elbow flexion and extension movements were made in the horizontal plane. The target was programmed to appear on the four corners of a fixed square. A trial could start at any time after the end of the previous one. At the start of a trial, the target appeared at the upper right corner of the square and the monkey had to move the cursor to it. Once reached, the target jumped counterclockwise to the next corner and the monkey had to move the cursor again, and so on (Fig. 5). A trial was registered as successful only if the monkey reached the target at all of the four corners within 5 s. In the data we collected, there were 240 successful trials; see Fig. 5 for some sample traces of hand movement in successful trials. A total of up to C = 49 MI single units were simultaneously recorded.

View larger version (14K): [in this window] [in a new window]   FIG. 5. Trajectories of hand movement of a macaque monkey in 5 successful trials of the target pursuit task as described under RESULTS. In each of the trials, the monkey’s hand passed the 4 corners of a square in counterclockwise order to reach a target that jumped from one corner to another, starting from the top right corner. Time when the top right corner was reached was defined as the onset of the movement sequence in a trial. Corners of the square are denoted by circles.

with parameter _i = (a_i, b_i). We shall mainly report results attained by incorporating the fitted Gamma distributions in the scoring. Later, we will show that similar results were attained as well without exploiting these distributions.

The temporal pattern of the spiking activity in MI has more variability than in RA. Taking this into account, the time bin duration was selected so that 99.9% of the time bins each contained at most one spike. Under this condition, the largest time bin duration was about 10 ms. We therefore chose = 10 ms. The binning reduced the size of the data while keeping most of its information. The hand position data were also down-sampled at 10-ms resolution to estimate the times when the corners were reached.

We used the movement sequences in the first 200 successful trials as the training data set to estimate the firing probabilities p_i,c(j) within 1 s around each event, with i = 1, ..., 4 indexing the events, c = 1, ..., 49 indexing the units, and j = –100, ..., 100 indexing the time bins within 1 s around a time of event. The parameters _i for the Gamma distributions were estimated as well. The estimates of p_i,c(j) and _i were then used to construct the template. The recording after the 200th successful trials was used to test the accuracy of the identification procedure.

The score function L(t) was computed according to Table 3 over the entire recording and then smoothed with a low-pass filter with cutoff frequency 0.5 Hz. Because successful trials occurred at a much lower frequency (0.21 Hz) than the cutoff frequency, the smoothing had no aliasing effect on the estimation of the onset times of the movement sequences (i.e., the times when the first corner was reached) in successful trials.

Across the entire recording period, 622 time locations t_i were identified where the smoothed score function (t) peaked locally (Fig. 6). Based on the assumption that the scores were higher near the onsets of the movement sequences in successful trials than elsewhere, only t_i with (t_i) > (t_i±1) were selected, resulting in 202 estimated onset times. For each of the selected t_i, the corresponding event times were estimated as in Table 3 (Fig. 7).

View larger version (38K): [in this window] [in a new window]   FIG. 6. Detection of movement sequences during successful trials in the target-pursuit task for the macaque monkey. Top: dark gray curve denotes the smooth score function using a low-pass filter on the original similarity scores over the entire recording period. Black dots mark all the estimated "max-maximum" points, which are local peak points of the smoothed score function with higher values than their adjacent local peaks. Gray vertical lines mark the actual onset times of movement sequences during successful trials (i.e., times of reaching the first corner). Bottom: sections extracted from the top panel showing details of the score function and identified onsets of movement sequences. Left panel: taken from the training data set (recording until the 200th successful trials). Right panel: test data set (the remaining 40 trials).

View larger version (37K): [in this window] [in a new window]   FIG. 7. Estimated times of events within identified movement sequences for the target-pursuit task experiment. In the top panel, the gray curve denotes the smoothed matching scores; the stars mark all the estimated starting points (first corner) of squared movements and the dots are the estimated points at the other 3 corners. Vertical lines mark the actual event times at the 4 corners of squared movements: the solid lines are for first corners and the dashed lines are for the other 3 corners. Bottom panels extracted from the top panel are 2 legible examples in the training (left) and test (right) data.

The power of the procedure was defined as the ratio of the number of true positives to that of successful trials. Evaluated over the entire data set, the power was high (173/240). The result was similar when the power was evaluated within the test data set only (29/40). The procedure also attained a high rate of true positives, defined as the ratio of the number of true positives to the number of selected t_i. For the test data set, the rate of true positives was 29/40. (Incidentally, within the test data set, the number of selected t_i was 40, the same as the number of successful trials.) The rate of true positives for the training data set was even higher (144/162).

To see how sensitive the above evaluation of power was to the threshold for the mean error, we varied the threshold value from 0.5 to 1.5 s, causing the power to change from 62 to 73% (Fig. 8A). The slight change in the power indicates that the evaluation was robust to the threshold selection thus justifies using 1 s as the threshold value.

View larger version (12K): [in this window] [in a new window]   FIG. 8. A: robustness of power with respect to varied thresholds ranging from 0.5 to 1.5 s. Power varied from 62 to 73% in the examined range. B: histogram of mean error of estimated event times associated with the true positives identified from the entire recording. For each true positive, the mean error is the average of the errors in the estimates of the times when the 4 corners were reached. Among a total of 173 true positives, 144 were identified from the training data set and the remainder from the test data set.

Finally, we tested the procedure without incorporating the estimated distributions of the IEIs in scoring. Equivalently, all G_i were set to 0. For the test data set, we again attained a high power (29/40) and a high rate of true positives (29/41). The mean error in the estimated event times was 0.226 ± 0.196 s and the errors associated with the four corners were 0.536 ± 0.607, 0.045 ± 0.062, 0.092 ± 0.189, and 0.231 ± 0.433 s. The results for the entire data set and those for the training data set were similar to these reported results. The results were only marginally inferior to those attained by exploiting the estimated distributions of the IEIs.

Comparisons with Kalman filter decoding method

A number of decoding methods have been proposed for MI of Macaque monkeys, including population vectors, linear filters, artificial neural networks, or Kalman filters (Carmena et al. 2003; Georgopoulos et al. 1986; Serruya et al. 2002; Wu et al. 2006). These methods have been demonstrated to be effective and accurate in estimating kinematic parameters such as hand movement direction or hand trajectory. However, none of them addresses the estimation of the timing of landmark events during a movement. This is largely attributed to the experiment paradigms these methods are designed for, which generate fairly straightforward hand movement of monkey, making it unnecessary to decode event timing during the movement. In contrast, in our squared movement paradigm, the monkey moved its hand following four corners of a square in each trial. In this paradigm, the times when the hand reached th