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/5/3763    most recent 00482.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 HighWire Citing Articles via ISI Web of Science (3) Citing Articles via Google Scholar Google Scholar Articles by Yu, B. M. Articles by Shenoy, K. V. Search for Related Content PubMed PubMed Citation Articles by Yu, B. M. Articles by Shenoy, K. V.

INNOVATIVE METHODOLOGY

Mixture of Trajectory Models for Neural Decoding of Goal-Directed Movements

¹Department of Electrical Engineering, ²Medical Scientist Training Program, ³Department of Neurosurgery, and ⁴Neurosciences Program, Stanford University, Stanford, California; and ⁵Gatsby Computational Neuroscience Unit, University College London, London, United Kingdom

Submitted 6 May 2006; accepted in final form 11 February 2007

	ABSTRACT

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION APPENDIX GRANTS ACKNOWLEDGMENTS REFERENCES

 
Probabilistic decoding techniques have been used successfully to infer time-evolving physical state, such as arm trajectory or the path of a foraging rat, from neural data. A vital element of such decoders is the trajectory model, expressing knowledge about the statistical regularities of the movements. Unfortunately, trajectory models that both 1) accurately describe the movement statistics and 2) admit decoders with relatively low computational demands can be hard to construct. Simple models are computationally inexpensive, but often inaccurate. More complex models may gain accuracy, but at the expense of higher computational cost, hindering their use for real-time decoding. Here, we present a new general approach to defining trajectory models that simultaneously meets both requirements. The core idea is to combine simple trajectory models, each accurate within a limited regime of movement, in a probabilistic mixture of trajectory models (MTM). We demonstrate the utility of the approach by using an MTM decoder to infer goal-directed reaching movements to multiple discrete goals from multi-electrode neural data recorded in monkey motor and premotor cortex. Compared with decoders using simpler trajectory models, the MTM decoder reduced the decoding error by 38 (48) percent in two monkeys using 98 (99) units, without a necessary increase in running time. When available, prior information about the identity of the upcoming reach goal can be incorporated in a principled way, further reducing the decoding error by 20 (11) percent. Taken together, these advances should allow prosthetic cursors or limbs to be moved more accurately toward intended reach goals.

	INTRODUCTION

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION APPENDIX GRANTS ACKNOWLEDGMENTS REFERENCES

 
Neural activity from motor cortical areas has been shown in a number of studies to be related to various aspects of the corresponding arm reach, including movement direction (Ashe and Georgopoulos 1994; Georgopoulos et al. 1986; Moran and Schwartz 1999; Riehle and Requin 1989; Tanji and Evarts 1976), movement extent (Riehle and Requin 1989), position (Ashe and Georgopoulos 1994; Paninski et al. 2004), velocity (Ashe and Georgopoulos 1994; Paninski et al. 2004), acceleration (Ashe and Georgopoulos 1994), posture (Caminiti et al. 1991; Scott and Kalaska 1997), speed (Moran and Schwartz 1999), joint angular velocity (Reina et al. 2001), force (Evarts 1968; Sergio and Kalaska 1998), and intended reach goal (Messier and Kalaska 2000; Shen and Alexander 1997). Although the coding scheme used by motor cortical areas is still incompletely understood (Fetz 1992; Moran and Schwartz 2000; Scott 2004), the regularities in the relationship between neural activity and aspects of the arm reach have enabled the development of direct brain-controlled prosthetic devices (Carmena et al. 2003; Chapin et al. 1999; Hochberg et al. 2006; Kennedy and Bakay 1998; Musallam et al. 2004; Santhanam et al. 2006; Serruya et al. 2002; Taylor et al. 2002). These devices are designed to allow disabled patients to regain motor function through the use of prosthetic limbs, or computer cursors, that are controlled by neural activity.

One of the key components of a neural prosthetic device is its decoding algorithm, which translates neural activity into desired movements. Examples of decoding algorithms that translate neural activity around the time of the movement (termed peri-movement activity) into continuous arm trajectories include population vectors (Taylor et al. 2002) and linear filters (Carmena et al. 2003; Hochberg et al. 2006; Serruya et al. 2002). Both of these decoding algorithms assume a linear relationship between the neural activity and arm state. In general, the arm state may include, but is not limited to, arm position, velocity, and acceleration.

Although these linear decoding algorithms are effective, recursive Bayesian decoders have been shown to provide more accurate trajectory estimates (Brockwell et al. 2004; Brown et al. 1998; Wu et al. 2004, 2006). Recursive Bayesian decoders are based on the specification of a probabilistic model comprising 1) a trajectory model, which describes how the arm state changes from one time step to the next, and 2) an observation model, which describes how the observed neural activity relates to the time-evolving arm state. If the modeling assumptions are satisfied, then Bayesian estimation makes optimal use of the observed data. Unlike the aforementioned linear decoding algorithms, recursive Bayesian decoders provide confidence regions for the arm state estimates and allow for nonlinear relationships between the neural activity and arm state.

The function of the trajectory model is to build into the recursive Bayesian decoder prior knowledge about the form of the reaches. The model may reflect 1) the hard, physical constraints of the limb (for example, the elbow cannot bend backward), 2) the soft, control constraints imposed by neural mechanisms (for example, the arm is more likely to move smoothly than in a jerky motion), and 3) the physical surroundings of the patient and his/her objectives in that environment. These statistical regularities and constraints can be learned by observing the reaching behavior and fitting the trajectory model to the actual reaches, which is the approach adopted in this paper. The degree to which the trajectory model reflects the kinematics of the actual reaches directly affects the accuracy with which trajectories can be decoded from neural data. A commonly used trajectory model is the random-walk model (Brockwell et al. 2004; Brown et al. 1998), which captures the fact that arm trajectories tend to be smooth. In other words, small changes in arm state from one time step to the next are more likely than large changes. The random-walk model is part of a family of trajectory models based on linear dynamics perturbed by Gaussian noise, which we refer to collectively as linear-Gaussian models. Linear-Gaussian models have been successfully applied to decoding the path of a foraging rat (Brown et al. 1998; Zhang et al. 1998), as well as arm trajectories in ellipse tracing (Brockwell et al. 2004), pursuit tracking (Shoham et al. 2005; Wu et al. 2004, 2006), and "pinball" tasks (Wu et al. 2004, 2006).

When selecting a trajectory model, one is typically faced with a trade-off between how accurately the trajectory model captures the movement statistics and the computational demands of its corresponding decoder. For example, for real-time applications, we may decide to use a relatively simple trajectory model because of its low computational cost, even if it fails to capture some of the salient properties of the observed movements. Even though we may be able to identify a more complex trajectory model that could yield more accurate decoded trajectories, the computational demands of the corresponding decoder may be prohibitive in a real-time setting. In this work, we present a general approach to constructing trajectory models that can exhibit rather complex dynamical behaviors, whose decoder can be implemented to have the same running time as simpler trajectory models. The core idea is to combine simple trajectory models, each accurate within a limited regime of movement, in a probabilistic mixture of trajectory models (MTM) (Kemere et al. 2003, 2004a,b).

We demonstrate the utility of this approach by developing a mixture of trajectory models suitable for goal-directed movements in settings with multiple goals. A common usage mode of real-time prosthetic systems involves guiding a computer cursor to acquire discrete goals in the subject's virtual workspace. This design approach was adopted in studies using multi-electrode neural recordings in monkeys (Carmena et al. 2003; Serruya et al. 2002; Taylor et al. 2002) and humans (Hochberg et al. 2006) and both electroencephalographic (EEG) (Wolpaw and McFarland 2004) and electrocorticographic (ECoG) (Leuthardt et al. 2004) recordings in humans. As shown by Donoghue and colleagues (Hochberg et al. 2006), the goal-directed cursor control design can allow a paralyzed patient to operate a computer interface controlling a variety of useful functions, including television and simulated email, thus illustrating its immediate clinical benefits. Because of the prevalence of goal-directed reaching in everyday life, this goal-directed design is likely to continue to be fruitful in systems involving prosthetic limbs, in addition to computer cursors, that are driven by the brain's activity. Indeed, many of the basic movements a paralyzed patient would desire are goal directed, such as reaching for a cup, picking up the phone, and feeding oneself. The utility of prosthetic systems based on goal-directed movements has fueled the development of statistical models and decoding algorithms tailored for goal-directed movements (Cowan and Taylor 2005; Kemere and Meng 2005; Kemere et al. 2002, 2003, 2004a,b; Srinivasan and Brown 2006; Srinivasan et al. 2005, 2006; Yu et al. 2005).

Goal-directed movements can be characterized by the following three properties. First, each movement is typically directed toward one of a (possibly large) number of discrete goals available in the subject's workspace. These goals may be visual targets presented on a computer screen or physical objects located near the subject. Second, repeated movements to the same goal are not identical. For example, there may be variability in movement speed or curvature. Third, the trajectories generally start at rest, proceed out to the desired goal, and end at rest. Previous studies considered goal-directed movements toward a single stationary (Kemere and Meng 2005; Srinivasan et al. 2005, 2006) or dynamic (Srinivasan and Brown 2006) goal with a known arrival time or stereotyped movements to multiple goals (Kemere et al. 2002, 2004b). In this work, we develop a mixture of trajectory models that captures all three properties of goal-directed movements and show how it can be used to decode movements from neural activity.

While the peri-movement neural activity is informative of the moment-by-moment details of the desired movement, there may be additional information available about the identity of the desired reach goal well before the desired time of movement onset. For example, if the phone rings, there is a greater chance that the goal of the upcoming reach will be the phone rather than the light switch. Even without such an external clue, the upcoming goal identity can often be inferred from neural activity related to motor preparation (termed delay activity because motor preparation is typically studied using a delayed-reach behavioral task) (e.g., Churchland et al. 2006b,c; Kurata 1993; Messier and Kalaska 2000; Riehle and Requin 1989; Shen and Alexander 1997; Weinrich and Wise 1982). The type of information conveyed by delay activity is categorically different from that provided by peri-movement activity. Whereas peri-movement activity specifies the moment-by-moment details of the arm trajectory (e.g., Ashe and Georgopoulos 1994; Moran and Schwartz 1999; Paninski et al. 2004; Schwartz 1992), delay activity has been shown to indicate the upcoming reach goal (Hatsopoulos et al. 2004; Musallam et al. 2004; Santhanam et al. 2006; Shenoy et al. 2003; Yu et al. 2004). It should be possible to use this goal information, when available, to improve the accuracy of the decoded trajectory. Brown and colleagues (Srinivasan et al. 2005, 2006) showed how to constrain free movement trajectories, given goal information that takes the form of a continuous distribution around a single goal. For multiple goals, the information usually takes the form of a discrete distribution across the possible goals. As with decoders we previously proposed (Kemere et al. 2002, 2003, 2004a,b), the MTM framework can naturally incorporate this goal information across multiple goals to improve the accuracy of the decoded trajectory. In contrast to a decoder that selects among a set of canonical trajectories (Kemere et al. 2002, 2004b), the MTM decoder can take into account behavioral variability across reaches to the same goal. Furthermore, the MTM decoder does not require the use of a linear filter, which was used in tandem with a mixture of hidden Markov models (Kemere et al. 2003) and a set of canonical trajectories with independent Gaussian noise at each time point (Kemere et al. 2004a).

We first present the MTM framework in its general form. Then, we construct an MTM that is appropriate for goal-directed reaches in settings with multiple goals and show how it can be used to decode arm trajectories from neural data. Next, we detail the behavioral task and neural recordings, along with how goal information can be extracted from delay activity. Finally, we compare the decoding accuracy of the MTM decoder with that using simpler trajectory models.

	METHODS

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION APPENDIX GRANTS ACKNOWLEDGMENTS REFERENCES

 
Mixture of trajectory models framework

Recursive Bayesian decoders require the specification of a trajectory model that describes the statistics of arm trajectories we expect to observe. Ideally, we seek to construct a complete model of neural motor control that captures the hard physical constraints of the limb (Chan and Moran 2006), the soft control constraints imposed by neural mechanisms, and the physical surroundings and context. One way to approximate such a complete model is to probabilistically combine trajectory models each of which is accurate within a limited regime of movement (Kemere et al. 2002, 2003, 2004a,b). Examples of movement regimes include different parts of the workspace, different reach speeds, and different reach curvatures. For the particular implementation tested here, each movement regime will correspond to movements heading toward a particular reach goal. At the onset of a new movement, the movement regime is unknown, or imperfectly known, and so the full trajectory model is composed of a mixture of the individual, regime-specific trajectory models. Here, we develop a recursive Bayesian decoder based on a mixture of trajectory models.

The task of decoding a continuous arm trajectory involves finding the likely sequences of arm states corresponding to the observed neural activity. At each time step t, we seek to compute the distribution of the arm state x_t given the peri-movement neural activity y₁, y₂, ... , y_t (denoted by {y}₁^t) observed up to that time. This distribution is written P(x_t|{y}₁^t) and termed the state posterior. Here, y_t is a vector of binned spike counts across the neural population at time step t, and t = 1 corresponds to the time at which we begin to decode movement. If the actual movement regime m is perfectly known before the reach begins, then we can compute the state posterior based only on the individual trajectory model corresponding to that regime. This distribution is written P(x_t|{y}₁^t, m) and termed the conditional state posterior. However, in general, the actual movement regime is unknown or imperfectly known, so we need to compute P(x_t|{y}₁^t, m) for each m {1, ... , M}, where M is the number of movement regimes (also referred to as mixture components).

To combine the M conditional state posteriors, we can simply expand P(x_t|{y}₁^t) by conditioning on the movement regime m

(1)

In other words, the state posterior is a weighted sum of the conditional state posteriors. The weights P(m|{y}₁^t) represent the probability that the actual movement regime is m, given the observed spike counts up to time t. Bayes' rule can then be applied to these weights in Eq. 1, yielding the key equation for the MTM framework

(2)

The conditional state posteriors P(x_t|{y}₁^t, m) and likelihood terms P({y}₁^t|m) in Eq. 2 can be computed or approximated using any of a number of different recursive Bayesian decoding techniques, including Bayes' filter (Brown et al. 1998), particle filters (Brockwell et al. 2004; Shoham et al. 2005), and Kalman filter variants (Wu et al. 2004, 2006). If available, prior information about the identity of the movement regime can be incorporated naturally into the MTM framework using P(m) in Eq. 2. This information must be available before the reach begins and may differ from trial to trial. If no such information is available, the same P(m) (e.g., a uniform distribution) can be used across all trials.

The computational complexity of the MTM decoder is M times that of computing P(x_t|{y}₁^t, m) and P({y}₁^t|m) for a particular mixture component m. Because the computations for each mixture component can theoretically be carried out in parallel, it is possible to set up the MTM decoder so that its running time remains constant, regardless of the number of mixture components M. In other words, the MTM approach enables the use of more flexible—and potentially more accurate—trajectory models without a necessary penalty in decoder running time. Furthermore, the MTM decoder preserves the real-time properties of its constituent estimators and is thus suitable for real-time prosthetic applications.

MIXTURE OF TRAJECTORY MODELS FOR GOAL-DIRECTED REACHES. The particular probabilistic model explored in this work is

(3)

(4)

(5)

where m {1, ... , M} indexes the reach goal and M is the number of reach goals. The dynamical arm state at time step t {1, ... , T} is x_t ^px1, which includes position, velocity, and acceleration terms, as specified in the APPENDIX. The corresponding observation, s_t–lagiⁱ {0, 1, 2,...}, is a peri-movement spike count for unit i {1, ... , q} taken in a time bin of width , where lag_i is the time lag (in time steps) between the neural firing of the ith unit and the associated arm state. For notational convenience, the spike counts across the q simultaneously recorded units are assembled into a q x 1 vector y_t, whose ith element is s_t–lagiⁱ. This is the y_t that appears in Eqs. 1 and 2. The parameters A_m ^pxp, b_m ^px1, Q_m ^pxp, _m ^px1, V_m ^pxp, lag_i , c_i ^px1, d_i do not depend on time and are fit to training data, as subsequently described.

Equations 3 and 4 define the trajectory model, which describes how the arm state x_t changes from one time step to the next. In this case, the full trajectory model is a mixture of linear-Gaussian trajectory models, each describing the trajectories toward a particular reach goal indexed by m. By this definition, each movement regime corresponds to movements heading toward a particular reach goal. Conditioned on the reach goal m, the trajectory model is a linear-Gaussian dynamical model.¹ Although the MTM framework will be illustrated in this work using Eqs. 3 and 4, mixtures of other trajectory models can also be used. For example, it is possible to define, for each reach goal, a linear-Gaussian model with a time-varying forcing term (Kemere and Meng 2005; Srinivasan et al. 2005, 2006) or a canonical trajectory (Kemere et al. 2002, 2004b).

Equation 5 defines the observation model, which describes how the observed peri-movement spike counts s_t–lagiⁱ relate to the arm state x_t. In Eq. 5, the linear mapping c'_ix_t + d_i is a cosine tuning model (Georgopoulos et al. 1982), where c_i is the "preferred state vector." This linear mapping is then passed through an exponential to ensure that the mean firing rate of the ith unit at time t – lag_i, exp(c'_ix_t + d_i), is non-negative. Note that, whereas each mixture component indexed by m in the trajectory model (Eqs. 3 and 4) can have different parameters leading to different arm state dynamics, the observation model (Eq. 5) is the same for all m.

Although the neural activity is known to be physically driving the trajectories, the probabilistic model Eqs. 3–5 specifies that the neural activity y_t is dependent on the arm state x_t. This model incorrectly implies, for example, that noise arising from the mechanical properties of the muscles that corrupts the arm trajectory should also show up in the neural activity in motor cortical areas. Nevertheless, models with this "inverted" structure have been shown to effectively decode arm trajectories (Brockwell et al. 2004; Shoham et al. 2005; Wu et al. 2004, 2006). The primary motivation for adopting such a structure is that there are established techniques for efficiently estimating an unobserved time series with known dynamics (in this case, the arm trajectory) from noisy observations (in this case, the neural spike counts). These techniques are detailed in the next section.

RECURSIVE BAYESIAN DECODING. Arm trajectories can be decoded from neural activity by applying Bayes' rule to the statistical relationships of Eqs. 3–5. Having observed the neural data, we seek the likely sequences of arm states that could have led to those neural observations. For each m and t, we need to compute the following two terms that appear in Eq. 2: the conditional state posteriors P(x_t|{y}₁^t, m) and the likelihood terms P({y}₁^t|m).

The conditional state posteriors can be obtained by iterating the following two updates. First, the one-step prediction is found by applying Eq. 3 to the conditional state posterior at the previous time step

(6)

Second, the conditional state posterior at the current time step is computed using Bayes' rule

(7)

Note that P(y_t|x_t, {y}₁^t–1, m) has been replaced by P(y_t|x_t) to obtain Eq. 7 because, given the current arm state x_t, the current observation y_t does not depend on the previous observations {y}₁^t–1 nor the reach goal m (cf. Eq. 5). The terms in the numerator of Eq. 7 are the observation model from Eq. 5 and the one-step prediction from Eq. 6. The denominator of Eq. 7 can be obtained by integrating the numerator over x_t, as shown later in Eq. 9.

When the trajectory and observation models are both linear-Gaussian, all of the relevant distributions are Gaussian and the integral in Eq. 6 can be computed exactly. Taking the iterations defined by Eqs. 6 and 7 is identical to applying the standard Kalman filter.

However, the particular observation model here (Eq. 5) is not linear-Gaussian. This leads to distributions that are difficult to manipulate and the integral in Eq. 6 cannot be computed analytically. We instead use a modified Kalman filter that uses a Gaussian approximation during the measurement update step (Eq. 7). We approximate the conditional state posterior as a Gaussian matched to the location and curvature of the mode of P(x_t|{y}₁^t, m), as detailed in the APPENDIX. This Gaussian approximation then allows the integral in Eq. 6 to be computed analytically because each mixture component of the full trajectory model (Eq. 3) is linear-Gaussian. This yields a Gaussian one-step prediction, which is fed back into Eq. 7.

The likelihood terms P({y}₁^t|m) can be expressed as

(8)

where

(9)

The integral in Eq. 9, which cannot be computed analytically, is approximated using Laplace's method (MacKay 2003). Note that this involves the same Gaussian approximation in x_t (i.e., the same mean and covariance) as made above for P(x_t|{y}₁^t, m).

EVALUATING PERFORMANCE. The state posterior P(x_t|{y}₁^t) in Eq. 1 is a multimodal distribution. To compare the performance of different decoders and to control a prosthetic cursor or arm, we need to collapse this multimodal distribution into a single decoded trajectory. In other words, we need to summarize the belief embodied in the state posterior with a single value _t at each time point. This can be done by first defining a loss function L, which specifies the loss incurred by the summary _t for each possible value of x_t. The single decoded trajectory is then the _t that minimizes the average loss under the state posterior

(10)

Here, we choose to use the instantaneous sum of squared distance loss function

(11)

in which case the _t that minimizes the average loss (Eq. 10) is the mean of the state posterior

(12)

In particular, the mean of the state elements corresponding to arm position is taken to be the decoded position trajectory. To compare different decoders, we first compute the root-mean-square position error (E_rms) between the actual and decoded trajectories on a per-trajectory basis. This yields a distribution of E_rms values for a given decoder. The E_rms distribution of different decoders can then be compared and statistics of each distribution (such as mean and SE) can be computed.

The expectation in Eq. 12 can be expanded by conditioning on the reach goal m as in Eq. 1, yielding

(13)

The interpretation of Eq. 13 is similar to that of Eq. 1. If the desired reach goal m is perfectly known before the reach begins, the decoded trajectory (_t) is computed based only on the individual trajectory model (i.e., the mixture component) corresponding to that reach goal. The decoded trajectory, in this case, is simply the mean of the conditional state posterior corresponding to the known reach goal, written E [x_t|{y}₁^t, m] and termed the component trajectory estimate for m. However, in general, the desired reach goal is unknown or imperfectly known, so we need to compute a component trajectory estimate E [x_t|{y}₁^t, m] for each m {1, ... , M}. The final decoded trajectory (_t) is a weighted sum of these component trajectory estimates, as shown in Eq. 13. As in Eq. 1, the weights P(m|{y}₁^t) represent the probability that the desired reach goal is m, given the observed spike counts up to time t.

In this work, we compare the performance of four decoders. The first is a state-of-the-art decoder presented by Kass and colleagues (Brockwell et al. 2004) based on a random-walk trajectory model (RWM) in acceleration. The trajectory (Eqs. A7 and A8) and observation (Eq. A9) models are defined in the APPENDIX. The second decoder is based on a single linear-Gaussian trajectory model (STM) shared across reaches to all goals. It is defined by Eqs. 3 and 4 for the special case of M = 1. The STM decoder uses the observation model shown in Eq. 5. The RWM and STM decoders provide points of comparison for the following two MTM decoders, both of which are based on Eqs. 3–5. Whereas the MTM_M decoder uses only peri-movement activity, the MTM_DM decoder uses both delay and peri-movement activity. In Eq. 2, the same P(m) (in this case, a uniform distribution) is used across all trials for MTM_M. In contrast, a different P(m) is used on each trial for MTM_DM based on the prior goal information extracted from delay activity.

Goal-directed reach task and neural recordings

Animal protocols were approved by the Stanford University Institutional Animal Care and Use Committee. We trained two adult male monkeys (Macaca mulatta, monkeys G and H) to perform delayed center-out reaches for juice rewards. As illustrated in Fig. 1A, visual targets were back-projected onto a frontoparallel screen 30 cm in front of the monkey. The monkey touched a central target and fixated his eyes on a crosshair at the upper right corner of the central target. After a center hold period of 500 (monkey G) or 400–600 ms (monkey H, selected randomly and uniformly in this range), a pseudorandomly chosen reach goal was presented at one of eight possible radial locations (30, 70, 110, 150, 190, 230, 310, 350°)² 10 cm away. After a pseudorandomly chosen instructed delay period of 200, 750, or 1,000 ms, the "go" cue (signaled by both the enlargement of the reach goal and the disappearance of the central target) was given and the monkey reached to the goal. After a hold time of 250 (monkey G) or 200 ms (monkey H) at the reach goal, the monkey received a liquid reward. The next trial started 200 (monkey G) or 100 ms (monkey H) later.

View larger version (30K): [in this window] [in a new window]   FIG. 1. Delayed reach task and neural recordings. A: task timeline (top), simultaneously recorded spike trains (middle), and arm and eye position traces (bottom) are shown for a single trial. Blue and red lines correspond to horizontal and vertical position, respectively. Full range of scale for the arm and eye position is ±15 cm from the center target. Trial from experiment H20041106.1. B: spatial arrangement of the 8 reach goals and corresponding spike histograms for one representative unit (Unit H20041217.23.0). Bars below histograms indicate delay (hatched) and peri-movement (gray) activity. Dotted lines denote reach goal onset and movement onset.

During experiments, monkeys sat in a custom chair (Crist Instruments, Hagerstown, MD) with the head braced and the nonreaching arm strapped to the chair. The presentation of the visual targets was controlled using the Tempo software package (Reflective Computing, St. Louis, MO). A custom photodetector recorded the timing of the video frames with 1-ms resolution. The position of the hand was measured in three dimensions using the Polaris optical tracking system (Northern Digital, Waterloo, Ontario, Canada; 60 Hz, 0.35-mm accuracy), whereby a passive marker taped to the monkey's fingertip reflected infrared light back to the position sensor. Eye position was tracked using an overhead infrared camera (Iscan, Burlington, MA; 240 Hz, estimated accuracy of 1°).

A 96-channel silicon electrode array (Cyberkinetics, Foxborough, MA) was implanted straddling dorsal premotor (PMd) and motor (M1) cortex (monkey G, right hemisphere; monkey H, left hemisphere), as estimated visually from local landmarks, contralateral to the reaching arm. Surgical procedures were described previously (Churchland et al. 2006c; Santhanam et al. 2006). Spike sorting was performed off-line using techniques described in detail elsewhere (Sahani 1999; Santhanam et al. 2004; Zumsteg et al. 2005). Briefly, neural signals were monitored on each channel during a 2-min period at the start of each recording session while the monkey performed the behavioral task. Data were high-pass filtered and a threshold level of three times the RMS voltage was established for each channel. The portions of the signals that did not exceed threshold were used to characterize the noise on each channel. During experiments, snippets of the voltage waveform containing threshold crossings (0.3 ms precrossing to 1.3 ms postcrossing) were saved with 30-kHz sampling. After each experiment, the snippets were clustered as follows. First, they were noise-whitened using the noise estimate made at the start of the experiment. Second, the snippets were trough-aligned and projected into a four-dimensional space using a modified principal components analysis. Next, unsupervised techniques determined the optimal number and locations of the clusters in the principal components space. We then visually inspected each cluster, along with the distribution of waveforms assigned to it, and assigned a score based on how well separated it was from the other clusters. This score determined whether a cluster was labeled a single-neuron unit or a multineuron unit.

Figure 1A shows the delayed reach task timeline, along with neural and behavioral data for a single trial with a lower-right reach goal. We later refer to the time between reach goal onset and the "go" cue as the delay period. Figure 1B illustrates the spatial arrangement of the eight reach goals, as well as the corresponding spike histograms for one representative unit across the eight reach goals. Each spike histogram was obtained by averaging the spike trains across multiple trials with the same reach goal. In broad terms, delay activity occurs during the delay period (always before movement onset), whereas peri-movement activity occurs around the time of the reach. The precise windows of delay and peri-movement activity used in this work are defined in later sections.

The monkeys were trained over several months and multiple data sets of the same behavioral task were collected. Each data set was collected in one day's recording session. For each monkey, we chose to analyze a data set with a large number of successful trials. The selected data sets consisted of 1,456 successful trials for monkey G (experiment G20040508) and 1,072 successful trials for monkey H (experiment H20041217), not including trials with 200-ms delay periods. The data set for monkey G (H) included 30 (56) single-neuron units and 68 (143) multineuron units, for a total of 98 (199) units.

The results reported in this work are cross-validated by randomly splitting the entire data set by trials into J roughly equal-sized parts. For each j {1, ... , J}, the jth part served as test data for a model trained on the other J – 1 parts. We used J = 9 (11) for the data set for monkey G (H). To evaluate decoder performance at different numbers of neural units, we further randomly split each part by units into K equal subparts. Each subpart contained the same number of trials and identical behavioral data as its parent, but with only 1/K of the neural data. To meaningfully compare the two data sets, we roughly equalized the number of units in each subpart. Unless otherwise specified, the results presented here assume K = 1 (98 units) for monkey G and K = 2 (99 units) for monkey H.

Incorporating goal information from delay activity

Up to this point, the neural activity discussed has been peri-movement activity, which takes place around the time of movement and specifies the moment-by-moment details of the arm trajectory. In the delayed-reach task, there is also neural activity present during an instructed delay period that directly precedes the "go" cue (termed delay activity). As shown in Crammond and Kalaska (2000) and Churchland et al. (2006c), neurons with delay activity are typically also active in the absence of an instructed delay during the reaction time period. Rather than specifying the moment-by-moment details of the trajectory, delay activity has been shown to reliably indicate the upcoming reach goal (Hatsopoulos et al. 2004; Musallam et al. 2004; Santhanam et al. 2006; Shenoy et al. 2003; Yu et al. 2004). The data sets for both monkeys G and H contain both delay and peri-movement activity on each trial. Furthermore, both types of activity may be emitted by the same unit on a single trial, as can be seen in Fig. 1.

The following describes how the reach goal can be decoded from delay activity by applying Bayes' rule. Let z be a q x 1 vector of spike counts across the q simultaneously recorded units in a prespecified time window during the delay period on a single trial. The distribution of spike counts (from training data) for each reach goal m can be fit to either a product of Gaussians (Maynard et al. 1999; Yu et al. 2004)

(14)

or a product of Poissons (Hatsopoulos et al. 2004; Shenoy et al. 2003)

(15)

where µ_i,m, _i,m², and _i,m are the parameters of the ith unit for the mth reach goal. The z_i notation in Eqs. 14 and 15 specifies that the distribution is describing the ith element of the vector z. In both models, the units are assumed to be independent given the reach goal. It would be natural to introduce conditional dependencies between the units using a general multivariate Gaussian, but there are often difficulties in estimating an invertible covariance matrix for tens to hundreds of units with a limited number of training trials (Maynard et al. 1999).

For any test trial, the probability that the upcoming reach goal is m given the delay activity z can be computed by applying Bayes' rule

(16)

where P(m) in Eq. 16 is assumed to be uniform. The most likely reach goal [i.e., the one with the largest P(m|z)] is usually taken to be the decoded reach goal (Hatsopoulos et al. 2004; Musallam et al. 2004; Santhanam et al. 2006; Shenoy et al. 2003; Yu et al. 2004).

The accuracy of the goal decoder (Eq. 16) varies with the duration and placement of the time window in which spikes are counted, as well as the precise spike count model P(z|m) that is used (Hatsopoulos et al. 2004; Santhanam et al. 2006). Optimizing the goal decoder is beyond the scope of this work and is treated in detail in the aforementioned references. Instead, we focus here on how to incorporate this goal information, if available, when decoding continuous arm trajectories. For this purpose, we choose to use the Gaussian model (Eq. 14) with a 200-ms spike count window starting 150 ms after the appearance of the reach goal.

The goal information from delay activity, P(m|z), can be incorporated naturally in the MTM framework in the place of P(m) in Eq. 2. The distribution P(m) in Eq. 2 represents the prior knowledge (i.e., before movement onset) that the upcoming reach goal is m. Because the delay activity entirely precedes movement onset and provides information about the upcoming reach goal, it can be used to set P(m) in Eq. 2 on a per-trial basis.

It is important to note that the most likely goal from Eq. 16 is not simply assumed here to be the goal of the upcoming reach. On a given trial, the delay activity may not definitively indicate the goal of the upcoming reach (e.g., two different reach goals may have significant probability) or it may indicate an incorrect goal for the upcoming reach. In this case, we would like to allow the subsequent peri-movement activity to determine the goal of the reach, or even correct the mistake, "in-flight." Instead of making a hard goal decode based on delay activity, the entire distribution P(m|z) is retained and passed to the MTM framework. For simplicity, we make the approximation that delay activity is informative only of the upcoming reach goal and is independent of the peri-movement activity; in other words, we assume that z is not directly coupled with x_t or y_t.

	RESULTS

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION APPENDIX GRANTS ACKNOWLEDGMENTS REFERENCES

 
In this section, we evaluate and compare the performance of decoding goal-directed movements using the RWM, STM, MTM_M, and MTM_DM decoders. For all decoders, we first fit the model parameters to training data, as detailed in the APPENDIX. The test data for a single trial consisted of 1) the arm trajectory, taken from 50 ms before movement onset to 50 ms after movement end at dt = 10-ms time steps; 2) the peri-movement spike counts, taken in nonoverlapping = 10-ms bins and temporally offset from the arm trajectory by the optimal lag found for each unit; and 3) the delay-period spike counts, taken in a single 200-ms bin starting 150 ms after the appearance of the reach goal. Arm trajectories in the test phase were used to evaluate the accuracy of the trajectories estimated from neural data. Because neural data collection ended shortly after movement end, the arm trajectories were not padded as in the training phase.

Figure 2 details, for a particular test trial (monkey G, 98 units), how the MTM decoded trajectory was obtained and compares the trajectory estimates produced by the different decoders. From Eq. 13, the MTM decoded trajectory is a weighted sum of component trajectory estimates E [x_t|{y}₁^t, m], one for each reach goal indexed by m {1, ... , 8}. In Fig. 2, B and C, the component trajectory estimates are plotted in the top panels, whereas the middle panels show how the corresponding weights P(m|{y}₁^t) evolved during the course of the trial.

View larger version (22K): [in this window] [in a new window]   FIG. 2. A representative test trial in which the use of delay activity improved the mixture of trajectory models (MTM) decoded trajectory. Top panels compare the actual trajectory (all panels, black) with the decoded trajectories for the (A) random-walk trajectory model (RWM, dark green) and single linear-Gaussian trajectory model (STM, light green); and a mixture of linear-Gaussian trajectory models, (B) one that uses only peri-movement activity (MTMM, red), and (C) one that uses both delay and peri-movement activity (MTMDM, orange). Ellipses denote 95% confidence intervals at 3 different time steps. Yellow squares represent the visual reach goals presented to the monkey in actual dimensions. Top panels in B and C also show the 8 component trajectory estimates for the MTM (cyan, blue, and magenta for the 3 components with the largest weights; gray for the other 5 components). Their corresponding weights, as they evolve during the trial, are plotted in the middle panels. Bottom panels compare the actual and estimated single-trial speed profiles using the same color conventions as in the top panels. Time zero corresponds to 60 ms before movement onset (i.e., one time step before we begin to decode movement). Note that the red and orange traces in the top panels are overlaid with the cyan trace. For this trial, root-mean-square position error (Erms) was 11.8, 26.9, 10.9, and 10.5 mm for the RWM, STM, MTMM, and MTMDM, respectively. Monkey G, 98 units (Experiment G20040508, trial ID 474).

If delay activity is available, it can be used to set a nonuniform P(m) in Eq. 2 on a per-trial basis, as previously discussed. The only difference between B and C in Fig. 2 is that the MTM decoder used delay activity in the latter, but not the former. In Fig. 2C (middle), the weights at t = 0 represent the probabilities of each reach goal based only on delay activity, before any peri-movement activity was observed. In this case, the delay activity indicated that the actual reach goal (cyan) was more probable than the other goals. This prior knowledge of the identity of the upcoming reach goal was then taken into account when updating the weights P(m|{y}₁^t) during the course of the trial as more and more peri-movement activity was observed. Note that using delay activity affected only P(m) in Eq. 2; the conditional state posteriors P(x_t|{y}₁^t, m) and the likelihood terms P({y}₁^t|m) remained unchanged. As the means of the conditional state posteriors, the component trajectory estimates therefore also remained unchanged, as can be verified by comparing Fig. 2, B and C (top). For the trial shown in Fig. 2, the use of delay activity reduced the competition between the actual reach goal (cyan) and the neighboring goals (blue and magenta). Compared with Fig. 2B (middle), the weight for the actual reach goal (cyan) in Fig. 2C (middle) was higher at every time point, the clearest effect seen during the first 200 ms. In other words, by using delay activity, the decoder was more certain of the actual reach goal throughout the trial. In Fig. 2C, a weighted sum of the eight component trajectory estimates (top) using these weights (middle) yields the MTM decoded trajectory (top, orange; E_rms: 10.5 mm).

By comparing the MTM decoded trajectories with the actual trajectory in Fig. 2, B and C (top), we see that the use of delay activity decreased the decoding error and tightened the confidence ellipses for this trial. The derivation of the MTM confidence intervals are given in the APPENDIX. Both MTM decoded trajectories had lower decoding error than the RWM (E_rms: 11.8 mm) and STM (E_rms: 26.9 mm), whose decoded trajectories are plotted in Fig. 2A (top). On this trial, the RWM decoder produced a reasonably accurate decoded trajectory, whereas the STM decoded trajectory proceeded slowly outward with wide confidence intervals.

These decoders can also be used to estimate the bell-shaped speed profile of the actual reach (Fig. 2, bottom). For the STM and MTM, we computed speed using its exact nonlinear relationship with the velocity elements in the state vector, rather than directly taking the speed element in the state vector, which involves a linear approximation. Compared with the speed profiles estimated by the RWM (dark green) and STM (light green) decoders, those estimated by the MTM decoders (red and orange) seem to better track the actual bell-shaped speed profile (black); this observation is shown even more clearly in the subsequent example trials.

In contrast to Fig. 2, Fig. 3 shows a trial (monkey G, 98 units) where the peri-movement activity alone was able to quickly determine the actual reach goal without much competition from neighboring goals. This can be seen in Fig. 3B (middle), where the weight corresponding to the actual reach goal (cyan) rose to unity after about 100 ms and stayed there for the remainder of the trial. As a result, the resulting MTM decoded trajectory (top, red; E_rms: 11.2 mm) was quite accurate. As in Fig. 2, we can incorporate delay activity if available; however, in this case, the dominant weight at t = 0 (blue) did not correspond to the actual reach goal (cyan), as seen in Fig. 3C (middle). In other words, the delay activity incorrectly indicated the identity of the upcoming reach goal. However, as these weights were updated by the observation of peri-movement activity, this "error" was soon corrected (within about 100 ms). From that point on, the weight corresponding to the actual reach goal dominated. Despite this error at the beginning of the trial, the MTM decoded trajectory in Fig. 3C (top, orange; E_rms: 13.0 mm) still headed to the correct goal and provided a reasonably accurate estimate of the arm trajectory. The larger confidence ellipses for MTM_DM compared with MTM_M reflect the competition between the actual (cyan) and neighboring (blue) reach goals. The decoded trajectories for the RWM (dark green) and STM (light green) are shown in Fig. 3A for comparison. As in Fig. 2, both MTM decoded trajectories yielded lower decoding error than the RWM (E_rms: 27.7 mm) and STM (E_rms: 24.2 mm). Furthermore, as shown in the bottom panels, the speed profiles estimated by the MTM decoders (red and orange) tracked the actual bell-shaped speed profile (black) more closely than those estimated by the RWM (dark green) and STM (light green) decoders.

View larger version (22K): [in this window] [in a new window]   FIG. 3. Representative test trial in which the peri-movement activity corrected an incorrect goal identification from delay activity. Figure conventions are identical to those in Fig. 2. For this trial, Erms was 27.7, 24.2, 11.2, and 13.0 mm for RWM, STM, MTMM, and MTMDM, respectively. Monkey G, 98 units (Experiment G20040508, trial ID 676).

Having demonstrated how the MTM framework produces trajectory estimates on individual trials, we can quantify and compare the average performance of the various decoders (RWM, STM, MTM_M, MTM_DM) across entire data sets. Figure 4 illustrates the following two main results, which hold true across both monkeys. First, a mixture of linear-Gaussian trajectory models (MTM_M) provides lower decoding error than either of the nonmixture trajectory models (RWM and STM) (Wilcoxon paired-sample test, P < 0.01). Compared with the STM decoder, the MTM_M decoder reduced E_rms from 22.5 to 13.9 mm (22.8 to 11.8 mm) in monkey G (H). Second, the use of prior goal information P(m) in the MTM framework (MTM_DM) can further decrease decoding error (Wilcoxon paired-sample test, P < 0.01). Compared with the MTM_M decoder, the MTM_DM decoder reduced E_rms from 13.9 to 11.1 mm (11.8 to 10.5 mm) in monkey G (H). Because the MTM decoder is inherently parallelizable (as described in METHODS), these performance gains can be obtained without an associated increase in decoder running time. The superior performance of the MTM_M compared with the RWM and STM can be explained by the fact that the MTM better captures the kinematics of goal-directed reaches. This can be seen in both the generative (prior) speed profiles (Fig. A1, bottom panels), as well as the decoded (posterior) speed profiles (Figs. 2 and 3, bottom panels). If delay activity is available, this additional source of information can be naturally incorporated in the MTM framework to further improve decoding performance (MTM_DM). The relative performance of the RWM and STM decoders is subsequently addressed in the context of Fig. 8.

View larger version (9K): [in this window] [in a new window]   FIG. 4. Erms (mean ± SE) comparison for the RWM, STM, MTMM, and MTMDM decoders. A: monkey G (98 units). B: monkey H (99 units).

View larger version (22K): [in this window] [in a new window]   FIG. A1. Position trajectories (top panels) and speed profiles (bottom panels). A: collected empirically, B: generated by the RWM, C: generated by the STM, and D: generated by the MTM. Only 24 reaches (3 to each reach goal) are shown in each column for clarity.

View larger version (15K): [in this window] [in a new window]   FIG. 8. Erms (mean ± SE) comparison of RWM (dark green), STM (light green), MTMM (red), and MTMDM (orange) decoders at different numbers of units. Dashed curves: monkey G, solid curves: monkey H. Vertical gray bar indicates the number of units used for the performance reported in Fig. 4.

View larger version (15K): [in this window] [in a new window]   FIG. 5. Two-dimensional histogram of Erms differences between pairs of decoders for (A) monkey G (98 units) and (B) monkey H (99 units). Horizontal axis: Erms difference between MTMM and STM; vertical axis: Erms difference between MTMDM and STM; diagonal axis: Erms difference between MTMDM and MTMM. Grayscale intensity (log scale) indicates the number of trials lying in each bin. Dotted lines represent the means of the Erms differences along each axis. Letters a, b, c, and d show where the trials in Figs. 2, 3, 6, and 7 lie on the histogram, respectively.

View larger version (22K): [in this window] [in a new window]   FIG. 6. Outlying test trial in which the MTMM decoded trajectory exhibited a snap-to-component effect. Figure conventions are identical to those in Fig. 2. For this trial, Erms was 45.3, 18.0, 37.1, and 10.4 mm for RWM, STM, MTMM, and MTMDM, respectively. Monkey G, 98 units (Experiment G20040508, trial ID 1921).

View larger version (20K): [in this window] [in a new window]   FIG. 7. Outlying test trial in which the peri-movement activity was not able to correct an incorrect reach goal identified by the delay activity. Figure conventions are identical to those in Fig. 2. For this trial, Erms was 27.9, 25.1, 10.8, and 49.2 mm for RWM, STM, MTMM, and MTMDM, respectively. Monkey G, 98 units (Experiment G20040508, trial ID 1608).