## Abstract

Brain-driven interfaces depend on estimation procedures to convert neural signals to inputs for prosthetic devices that can assist individuals with severe motor deficits. Previous estimation procedures were developed on an application-specific basis. Here we report a coherent estimation framework that unifies these procedures and motivates new applications of prosthetic devices driven by action potentials, local field potentials (LFPs), electrocorticography (ECoG), electroencephalography (EEG), electromyography (EMG), or optical methods. The brain-driven interface is described as a probabilistic relationship between neural activity and components of a prosthetic device that may take on discrete or continuous values. A new estimation procedure is developed for action potentials, and a corresponding procedure is described for field potentials and optical measurements. We test our framework against dominant approaches in an arm reaching task using simulated traces of ensemble spiking activity from primary motor cortex (MI) and a wheelchair navigation task using simulated traces of EEG-band power. Adaptive filtering is incorporated to demonstrate performance under neuron death and discovery. Finally, we characterize performance under model misspecification using physiologically realistic history dependence in MI spiking. These simulated results predict that the unified framework outperforms previous approaches under various conditions, in the control of position and velocity, based on trajectory and endpoint mean squared errors.

## INTRODUCTION

Amyotrophic lateral sclerosis, spinal cord injury, brain stem infarcts, advanced-stage muscular dystrophies, and diseases of the neuromuscular junction profoundly disrupt voluntary muscle control. New technologies, variously called brain-computer interfaces (Leuthardt et al. 2004; Wolpaw and McFarland 2004), motor neural prostheses (Carmena et al. 2003; Hochberg et al. 2006; Serruya et al. 2002; Taylor et al. 2002), and cognitive prostheses (Musallam et al. 2004; Shenoy et al. 2003), represent a communication link that bypasses affected channels of motor output. Manually actuated devices, eye tracking, and other approaches (Frey et al. 1995) represent practical solutions for many patients but may not be feasible for individuals with profound motor deficits. Moreover, brain-driven interfaces, based on a variety of actuation technologies including functional electrical stimulation, active orthotics, exoskeletons, robotic arms, and wireless consumer electronics, have the potential to provide dexterous and natural control without muscle fatigue.

A brain-driven interface includes a method to monitor neural activity, an algorithm to map neural activity into intended device behavior, a device to be controlled, and a feedback mechanism from the device to the user (Andersen et al. 2004; Kubler et al. 2006; Lebedev and Nicolelis 2006; Schwartz et al. 2006; Wolpaw et al. 2002). This paper relates to the optimal mapping between preprocessed neural activity and estimates of the user's intention that determine control signals. The method presented here unifies four canonical approaches, demonstrates new applications, and suggests one path to further algorithm development.

In prosthetics literature, the optimal mapping is predominantly described as an estimation (filtering) problem followed by a control problem. First, an estimate of the user-intended prosthetic device state is calculated based on neural activity that serves as a noisy observation of that intention. Second, a control law determines inputs to the device that achieve this estimate of the user-intended device state. This optimization ignores feedback to the user but provides a practical approach that is accommodated within the existing framework of estimation theory or similarly, a tracking problem in stochastic control. The word “intention” is used to describe an objective behavior of the prosthetic device that the subject is driven to exhibit through mechanisms such as reward or gratification that may be intrinsic and extrinsic to the brain.

Previous approaches to the estimation problem include: manually adjusted linear combinations of power spectral band energies (Wolpaw and McFarland 1994), population vectors for automated but suboptimal linear mappings (Georgopoulos et al. 1986), linear regression for optimized linear mappings (Wessberg et al. 2000), support vector machines (Olson et al. 2005), and recursive Bayesian estimation procedures, including the Kalman filter (Wu et al. 2006), particle filter (Brockwell et al. 2004), and point process filter (Brown et al. 1998; Eden et al. 2004a).

Bayesian estimation allows better tracking than linear regression (including population vector based decoding) in off-line data analyses in various neural decoding algorithm studies (Barbieri et al. 2005; Brockwell et al. 2004; Brown et al. 1998; Eden et al. 2004a), although this improvement may not be universally observed in practice (Kim et al. 2006). This approach describes the intended state of a prosthetic device and observed neural activity as a sequence of random variables indexed by time. The trajectory model defines a prior on the sequence of intended device states. The observation model defines the relationship between neural activity and intended device states. Actual device states are determined from neural activity based on these trajectory and observation models.

Neural prosthetic device algorithms are assessed a variety of ways, including simulated off-line decoding (Eden et al. 2004a,b; Kemere and Meng 2005; Kemere et al. 2004; Srinivasan and Brown 2007; Srinivasan et al. 2005, 2006), off-line decoding of actual neural activity recorded during movements in laboratory animals (Brockwell et al. 2004; Hatsopoulos et al. 2004; Wessberg et al. 2000), and on-line prototypes involving laboratory animals or humans (Carmena et al. 2003; Hochberg et al. 2006; Leuthardt et al. 2004; Musallam et al. 2004; Serruya et al. 2002; Shenoy et al. 2003; Taylor et al. 2002; Wahnoun et al. 2006; Wolpaw and McFarland 2004; Santhanam et al. 2006). One ultimate gold standard in comparing methods is a double-blinded longitudinal study of device performance in activities of daily living for a population of individuals that experience well-defined stages and pathological mechanisms of a target motor deficit. Among these various methods of assessment, simulation provides an important first step in clarifying essential design concepts, demonstrating mathematical approaches, providing intuition about performance characteristics, and preliminarily assessing relative performance of competing algorithms.

In the following sections, we present an estimation framework for brain-driven interfaces that explicitly allows the designer to span a full range of device capabilities by employing a hybrid state space composed of interacting discrete and continuous valued random processes. This method is shown to generalize previous Bayesian approaches to prosthesis design, including finite state machines (Shenoy et al. 2003), free-arm-movement models (Wu et al. 2006), reaching-movement-trajectory models (Cowan and Taylor 2005; Kemere and Meng 2005; Kemere et al. 2003; Srinivasan and Brown 2006; Srinivasan et al. 2006a; Wu et al. 2004), the switching observation model (Wu et al. 2004), and the mixture-of-trajectories model (Yu et al. 2007). One possible filtering procedure is derived for point process observations on the hybrid state space, and connections are drawn to existing literature on hybrid estimation for Gaussian observation processes (switching Kalman filters).

To demonstrate the versatility of this framework, two prosthetic device applications are described in simulation. The first application uses ensemble spiking activity from point process models of primary motor cortex (MI) to coordinate reaching movements to a target that may change within a discrete set of targets over the course of the movement. The second application uses power in EEG bands to navigate a wheelchair with definitive stopping. The reaching example is further expanded into a third application that demonstrates adaptive filtering under neuron death and discovery. In the final application, the effect of model misspecification is examined by incorporating physiologically realistic history dependence into simulated MI spiking activity.

## METHODS

### Elements of the hybrid framework

In the formulation of the neural prosthesis estimation problem, the user communicates the intended state of the prosthetic device via neural signals. The optimal brain-driven interface must convert these neural signals into an estimate of the intended device state that minimizes some distance metric (cost) to the intended sequence of device states. The cost is commonly assumed to be some form of mean squared error for continuous-valued device states and frequency of proper classification for discrete-valued states (Andersen et al. 2004; Kubler et al. 2006; Lebedev and Nicolelis 2006; Schwartz 2004; Wolpaw et al. 2002). Implicit in this formulation, a controller is subsequently expected to receive the estimate and drive the device to the corresponding state with the required precision and response time.

To maintain generality, we describe the user-intended device state at time step *k* by a vector of discrete random variables *s _{k}* and continuous random variables

*x*. The user drives neural activity

_{k}*n*

_{k}

^{1:C}from

*C*channels at time step

*k*based on the desired device states

*s*and

_{k}*x*. Although we refer to

_{k}*n*

_{k}

^{c}as the activity of the

*c*th neuron in the specific context of multiple single-unit recording,

*n*

_{k}

^{c}may correspond more generally to the

*c*th signal of any measurement of activity at time step

*k*, including single neuron spiking, multiunit activity, continuous electric field measurements, and even eye movements. Neural activity that leads or lags the desired device states can be incorporated by supplying the appropriate advanced or delayed neural activity in place of the current neural activity. As an example of discrete and continuous-valued states in the context of a classical center-out reaching task, delay period spiking activity in dorsal premotor cortex can be described as neural activity driven by the upcoming target [represented by a discrete set of states (Hatsopoulos et al. 2004; Shenoy et al. 2003; Srinivasan et al. 2006b)], and peri-movement spiking activity in primary motor cortex can be described as neural activity driven by arm velocity [represented by a continuous set of states (Georgopoulos et al. 1986; Hatsopoulos et al. 2004; Moran and Schwartz 1999; Serruya et al. 2002; Wessberg et al. 2000)]. However, neural observations of both discrete and continuous states are not necessary. The history of activity at this time step,

*H*= (

_{k}*n*

_{1}

^{1:C},

*n*

_{2}

^{1:C}, …,

*n*

_{k−1}

^{1:C}), may also affect

*n*

_{k}

^{1:C}due to recurrent neural connections and other sources of history dependence.

As an illustration of these variables, consider driving a car with your EEG instead of with your arms. Your intention to increment or decrement the gear as well as the current gear position can be captured by a discrete variable *s _{k}*, whereas the desired wheel or gas pedal angles can be further described by the continuous variable

*x*that evolves depending on the resulting gear position recorded in

_{k}*s*

_{k}_{+1}. The EEG amplitudes on

*C*different channels, indicated by

*n*

_{k}

^{1:C}, may depend on your discrete and continuous-valued intentions for the car but also the history of previous amplitudes

*H*because EEG typically exhibits a power spectral density that is shaped rather than flat. Note that the intended device state need not correspond literally to parts of the car as with the intention to increment or decrement the gear. The choice of states can dramatically impact ease of use, just as with the design of an interface to a consumer electronic device such as the MP3 player.

_{k}We define the hybrid state space as a joint probability density on the entire temporal sequence (trajectory) of intended states and neural activity *p*(*x*_{1}, *s*_{1}, *n*_{1}^{1:C}, *x*_{2}, *s*_{2}, *n*_{2}^{1:C}, …). Graphical models on acyclic graphs are a pictorial description of this joint density (Jordan 2004). The graphical model allows us to constrain the form of the joint density with a simple and consistent prescription. Consider our specific graphical model of the hybrid state space (Fig. 1*A*) in which we have illustrated only one segment of the full graphical model that corresponds to the entire trajectory. The circles, called nodes, denote random variables corresponding to the intended states and neural activity. The arrows specify interdependencies between the random variables. A consistent prior distribution on the entire set of nodes is provided by specifying distributions for each node conditioned on its parents, which are all nodes at the base of arrows that point to that node. Nodes without parents require unconditioned priors. The graphical model imposes a Markov structure, where any node is independent of all other nodes when conditioned on its parents. The hanging arrows directed toward *n*_{k}^{1:C} and *n*_{k+1}^{1:C} represent history dependence in the neural activity.

The probability distribution *p*(*n*_{k}^{1:C}|*x _{k}*,

*s*,

_{k}*H*) associated with this hybrid state space corresponds to the observation model because it relates the present measurement of neural activity to the present intention of the user and the history of neural activity. The probability distributions

_{k}*p*(

*x*

_{k}_{+1}|

*x*,

_{k}*s*

_{k}_{+1}) and

*p*(

*s*

_{k}_{+1}|

*s*) comprise the trajectory model; they describe the frequency and types of transitions in user intent for which the prosthetic device is designed.

_{k}The discrete probability distribution *p*(*s _{k}*

_{+1}|

*s*) in the trajectory model can be conveniently summarized as a discrete state transition matrix

_{k}**M**(1) This notation means that the entry in the

*i*th row and

*j*th column of

**M**corresponds to

*p*(

*s*

_{k}_{+1}=

*i*|

*s*=

_{k}*j*), the probability of desiring state

*i*at time step

*k*+ 1 given that state

*j*was desired at time step

*k*.

In principle, the observation model should properly describe the relation between neural observations and user intent, and the trajectory model should accurately reflect the distribution of user intents. Model mismatch describes errors that accumulate from an incorrect model specification. Whereas continuous field potentials (LFP, ECoG, EEG) are typically described by Gaussian observation models (Tarvainen et al. 2004), spiking activity at millisecond resolution is better described by point process observation models (Brown 2005; Brown et al. 2002; Daley and Vere-Jones 2003; Snyder and Miller 1991). The continuous component *p*(*x _{k}*

_{+1}|

*x*,

_{k}*s*

_{k}_{+1}) of the trajectory model can often be reasonably approximated as Gaussian to anticipate smooth changes in the user's continuous state intent when conditioned on a particular discrete state. The discrete component of the trajectory model

*p*(

*s*

_{k}_{+1}|

*s*), also called the state transition density, is generally defined by specifying a probability between 0 and 1 for each possible pair of intentions (

_{k}*s*

_{k}_{+1}, s

_{k}), although parameterization may be relevant to dealing with a large set of discrete intentions.

Alternate connections could have been used to describe the hybrid state space. The specific choice of connections made (Fig. 1*A*) draws on a standard form used in hybrid filtering on Gaussian observations (Murphy 1998; citeseer.ist.psu.edu/article/murphy98switching.html) but extends it to accommodate arbitrary history dependence. This imposed structure on the state space makes it easy to obtain estimates of the intended device state in a recursive fashion based on the latest set of neural activity. Moreover, the connections are sufficiently general as to accommodate a diverse set of applications.

Five previously described Bayesian approaches to neural prosthetics fall within this single framework (Fig. 1, *B–E*). A finite state machine description of the prosthesis (Shenoy et al. 2003) consists of a sequence of discrete user-intended states, rules for transitions between those states, and a relationship between states and neural signals (Fig. 1*B*). Free-arm-movement models (Wu et al. 2006) and reaching-movement-trajectory models (Cowan and Taylor 2005; Kemere and Meng 2005; Kemere et al. 2003; Srinivasan and Brown 2007; Srinivasan et al. 2005; Srinivasan et al.; Yu et al. 2005b), both describe continuous-valued arm-movement intentions that drive neural activity (Fig. 1*C*). The switching observation model (Wu et al. 2004) accommodates poorly sorted neural activity that might be better described by combinations of single-cell receptive fields (Fig. 1*D*). For example, multimodal multiunit activity can be described as a probabilistic combination of two or more unimodal tuning curves. The mixture-of-trajectories model (Yu et al. 2007) was designed for continuous-valued reaching movements to a stationary target drawn from a discrete set (Fig. 1*E*).

In subsequent discussion, we use some of these model names more generally to describe their use in situations for which they were not originally applied. We use “free movement” in reference to free-arm-movement models, but where the continuous state is not necessarily related to an arm, but rather some application-specific state, like wheelchair velocity. Similarly, we use “mixture of trajectories” more generally to denote any model with a fixed discrete state regardless of the specific application.

Because the hybrid framework is sufficiently general, these various applications can be implemented by appropriately defining states, the trajectory model, and the observation model. For example, the hybrid framework is equivalent to the free-movement model when defining only one possible discrete state, and applying a free-movement model to define the continuous component *p*(*x _{k}*

_{+1}|

*x*,

_{k}*s*

_{k}_{+1}) of the trajectory model. Similarly, the hybrid framework is equivalent to the mixture-of-trajectories model when the discrete state transition matrix

**M**is assigned to be the identity matrix

**I**. With

**M**=

**I**, the probability of transitioning to a different state from the current state is zero, capturing the notion of a static discrete state.

Although the hybrid state space depicted (Fig. 1) unifies these previous conceptions of neural prosthesis design, an estimation procedure (filter) must still be specified to generate probability densities of intended device states given neural activity from which average cost measures can be minimized. The resulting hybrid filter procedure can then be used to describe filtering in any of these previous conceptions. In the following sections, we develop a point process filter for spiking observations in the hybrid state space and review the corresponding Gaussian process filter for continuous field potentials.

### Point process models of ensemble spiking activity

Neural activity in the form of action potentials is a sequence of transient spiking events. We first specify an observation model that captures the quality of temporally localized events as well as possible dependencies between neurons in an ensemble.

Signals of this nature are naturally described by point process observation models (Brown 2005; Brown et al. 2002; Daley and Vere-Jones 2003; Snyder and Miller 1991). The crux of the point process description of a single neuron is its conditional intensity function. This is the instantaneous probability of firing as a function of elapsed time *t* and generally conditioned on continuous-valued signals *x*(*t*), discrete-valued signals *s*(*t*), and spiking history *H*(*t*), where *N*(*t*) denotes the total number of spikes generated by the neuron since some arbitrary starting time (2) We introduce additional notation to accommodate a population of neurons in a discrete time setting. For the *k*th discrete time step of length δ_{k} seconds, the conditional intensity of neuron *c* is represented as λ_{k}^{c} in units of spikes per second. Spiking activity at the *k*th time step is summarized by a vector *n*_{k}^{1:C} = (*n*_{k}^{1}, *n*_{k}^{2}, …, *n*_{k}^{C}) of binned spike counts. The *c*th element of *n*_{k}^{1:C} contains the total number of spikes generated by the *c*th neuron in the respective δ_{k}-second interval. The observation model for the spiking activity *n*_{k}^{1:C} of an ensemble of *C* neurons binned at δ_{k} second intervals is approximated (Eden et al. 2004a) as follows (3) where *H _{k}* = (

*n*

_{1}

^{1:C},

*n*

_{2}

^{1:C}, …,

*n*

_{k−1}

^{1:C}) is the history of spiking activity at step

*k*for the ensemble. This is an approximation in two regards. First, neurons are assumed to be independent conditioned on the history of population activity and current intended device state. This assumption still captures a causal notion of probabilistic dependence among neurons, that for example, the spiking history of one neuron might affect the present spiking probability of another neuron. Second, the discrete time observation model in 3 approximates the exact continuous-time observation model for a point process (Truccolo et al. 2005).

### Filtering spikes with the hybrid framework

To develop an estimation procedure that maps spikes to hybrid device states, we looked to the switching Kalman filter (Bar-Shalom et al. 2001; Murphy 1998; citeseer.ist.psu.edu/article/murphy98switching.html) which maps Gaussian signals to hybrid device states. We could possibly bin the spike trains (lump them into sequential intervals of time) and then apply a standard switching Kalman filter. However, spike trains that have been binned (lumped into sequential intervals of time) may only satisfy the Gaussian assumption as the bin size grows. This results in a tradeoff between the user's control of when an action is supposed to happen versus how it is supposed to happen.

To avoid this tradeoff, we used the point process observation model (previous section) as a statistical description of spiking that is accurate on a millisecond-by-millisecond time scale (Brown et al. 1998; Truccolo et al. 2005). This necessitated the development of a point process filter for hybrid states. Just as there are several approaches to the switching Kalman filter that balance computational complexity and accuracy (Bar-Shalom et al. 2001), there are several possible ways to filter spikes for the hybrid framework.

Our point process filter is adapted from a mixture-of-Gaussians switching Kalman filter called the Interacting Multiple Model (IMM) (Bar-Shalom et al. 2001). The IMM has been a popular choice to balance complexity and accuracy for a wide array of Gaussian applications that track multiple targets with multiple sensors, such as camera-based human motion tracking (Farmer et al. 2002), and radar-based airplane tracking (Mazor et al. 1998). The specific constraints of an application dictates the proper balance between complexity and accuracy, which is commonly determined through empirical testing. In the prosthesis application, demands for hardware memory, energy consumption, and computation speed must be weighed against the accuracy needed for the user to achieve ease of use and reliability. Although the IMM cannot be applied directly to point process observations, it informs our approach to hybrid filtering for spikes because prosthesis applications share the same essential structure as previous multi-target multi-sensor applications. We summarize the point process hybrid filtering procedure in the following section.

By combining this observation model and the hybrid state space defined in the previous section, we derive a specific filter to estimate hybrid device states from ensemble activity that has been modeled as a point process (appendix c). The exact, continuous-time solution to this filtering problem is a stochastic partial differential equation known as point process filtering with jump-Markov processes (McGarty 1974). The method in this section is one possible approximation.

Our point process filter derivation first manipulates probability densities without specifying their functional form and later introduces the functional form of the point process observation model given in *Eq. 3.* The resulting point process filter retains the same flavor as the IMM filter. Just as the IMM filter involves a bank of Kalman filters that run in parallel, our hybrid filter employs a bank of stochastic state point process filters (Eden et al. 2004a) that run in parallel, one for each possible value of the discrete state at a particular time step. A practical note on numerical issues for implementation is given in appendix f.

### Spike filtering with the hybrid framework in nine steps

Each iteration of the point process hybrid filter involves nine basic steps. The quantities *p*(*s _{k}*|

*H*

_{k}_{+1}) and

*p*(

*x*|

_{k}*s*,

_{k}*H*

_{k}_{+1}) come from the previous iteration, where

*p*(

*s*

_{1}) and

*p*(

*x*

_{1}) are used instead for the first iteration. A Gaussian approximation to a probability density on the continuous state

*x*is specified by a mean and covariance matrix. A probability mass function on the discrete state

_{t}*s*is specified by a list of probabilities for each possible value of

_{k}*s*. See appendix f for a practical note on numerical issues.

_{k}### Step 1

Compute

### Step 2

Compute

### Step 3

Approximate with the Gaussian approximation to mixtures of Gaussians (see methods).

### Step 4

Calculate the Gaussian approximation to *p*(*x _{k}*

_{+1}|

*s*

_{k}_{+1},

*n*

_{k+1}

^{1:C},

*H*

_{k}_{+1}). Specifically, for each value that

*s*

_{k}_{+1}can take on, send

*p*(

*x*|

_{k}*s*

_{k}_{+1},

*H*

_{k}_{+1}) through one full iteration of a point process filter (see Appendix A) with observation equation

*p*(

*n*

_{k+1}

^{1:C}|

*x*

_{k}_{+1},

*s*

_{k}_{+1},

*H*

_{k}_{+1}) and state equation

*p*(

*x*

_{k}_{+1}|

*x*,

_{k}*s*

_{k}_{+1}). Retain these densities (1 Gaussian for each possible value of

*s*

_{k}_{+1}) for the next iteration.

### Step 5

Calculate . Note that and are posterior and prediction covariance terms from step 4 (see also appendix a). This is the Laplace approximation (see appendix e). Here, λ_{k+1}^{c} denotes the conditional intensity of neuron *c* at time step *k* + 1, which may be a function of *s _{k}*

_{+1},

*x*

_{k}_{+1}, and

*H*

_{k}_{+1}.

### Step 6

Calculate . Retain this density for the next iteration.

### Step 7

Calculate using the results from steps 4 and 6.

### Step 8

Choose the discrete and continuous device states for step *k* +1 based on steps 6 and 7 and your cost function. For example, to approximately minimize average classification error, choose the value of *s _{k}*

_{+1}that maximizes step 6. To approximately minimize mean squared error, choose the average value of

*x*

_{k}_{+1}under the density calculated in step 7.

### Step 9

Return to step 1 for the next time step.

### Filtering continuous field potentials with the hybrid framework

Continuous field potentials are also viable sources for the control of prosthetic devices, such as with EEG (Wolpaw and McFarland 2004), ECoG (Leuthardt et al. 2004), and LFP (Pesaran et al. 2006). An EEG-based device has the potential for wide application because it is completely noninvasive. The ECoG and LFP approaches may allow cheaper and more robust hardware solutions than spike-driven interfaces because skull screws and coarse electrodes may suffice for these signals where micromachined multiunit arrays are needed to record stable ensemble spiking activity.

The physiological basis of these continuous field potentials is varied and different from that of ensemble spiking activity. Additional research is needed to understand effective training paradigms and hardware design as they pertain to each of these signal sources. However, existing filtering procedures are sufficient to incorporate these signals into the hybrid framework (Bar-Shalom et al. 2001). This is because continuous field potentials have been extensively modeled as Gaussian observation processes, including autoregressive moving average (ARMA) models (Tarvainen et al. 2004). As a result, the many types of switching Kalman filter can be applied directly to accommodate these signals into the hybrid framework.

The interacting multiple model (IMM) (Bar-Shalom et al. 2001) is the switching Kalman filter that is analogous to the point process filter presented in the previous section. The IMM derivation can be written in almost the same fashion, except that the observation model 3 is now Gaussian. To specify the Gaussian observation model, represent the *C* continuous field potentials at time step *k* in a vector *n*_{k}^{1:C} = [*n*_{k}^{1}, *n*_{k}^{2}, …, *n*_{k}^{C}]′. The entries in this vector could be field potential amplitudes, power in specific frequency bands, or some more complicated transformation of the raw signals. The Gaussian observation model relates *n*_{k}^{1:C} to the intended device states, waveform history *H _{k}*, and zero-mean Gaussian observation noise ε

_{k}as follows (4) The covariance matrix of ε

_{k}is denoted

*R*. Here,

_{k}*x*is a

_{k}*J*× 1 vector of continuous states,

*D*is a

_{k}*C*×

*J*matrix that may depend on

*s*, and

_{k}*f*(

_{k}*s*,

_{k}*H*) is a function that maps neural history to a

_{k}*C*× 1 vector, such as with ARMA models. The functional form of

*f*and the remaining parameters of the observation model may in general depend on the discrete intended state

_{k}*s*.

_{k}The IMM procedure is simply the point process hybrid filter procedure of the previous section, with modifications to steps 4 and 5. In step 4, Kalman filters are applied instead of point process filters. These are the Gaussian filter *Eqs. A6* and *A7* in appendix a. In step 5, *p*(*n*_{k+1}^{1:C}|*s _{k}*

_{+1},

*H*

_{k}_{+1}) is calculated based on 4 and the prediction density

*p*(

*x*

_{k}_{+1}|

*s*

_{k}_{+1},

*H*

_{k}_{+1}) from step 4. Specifically,

*p*(

*n*

_{k+1}

^{1:C}|

*s*

_{k}_{+1},

*H*

_{k}_{+1}) is a separate Gaussian for each possible value of

*s*

_{k}_{+1}with distribution given in

*Eq. A8.*

These two modifications, together with the hybrid filter procedure described in the previous section, complete the description of the IMM filter. The IMM filter is illustrated in the following text (results, *Application 2*) with a wheelchair navigation task with definitive stopping based on EEG-band energy.

## RESULTS

We now illustrate four applications of the hybrid filtering framework to brain-driven interfaces: *1*) reaching to discrete targets that switch during movement based on MI ensemble spiking activity; *2*) wheelchair navigation with definitive stopping based on EEG; *3*) adaptive filtering of MI ensemble spiking with ongoing neuron death and discovery; and *4*) filtering under model mis-specification due to physiologically realistic history dependence in MI ensemble spiking.

In each section, we illustrate the relative performance of three filters described in methods: free-movement estimation (Fig. 1*C*) (Wu et al. 2006) mixture of trajectory estimation (Fig. 1*E*) (Yu et al. 2007), and the hybrid filtering framework (Fig. 1*A*).

### Application 1: reaching to discrete targets that switch during movement

The first application employs the hybrid framework to track a reaching movement with a mid-flight change in the desired target. In the switching-target-reach task, the subject is required to reach to one of a discrete set of targets with a prosthetic arm driven by ensemble spiking activity from MI. Each reach must be completed within 2 s in a two-dimensional plane from the origin to one of eight targets arranged evenly on a circle of 0.25 m radius. In addition, the target changes once during the course of the movement, requiring the user to make a corrective maneuver to the new target location. The initial target and final target are chosen at random with equal probability over the target set. The switch time, unknown to the user, is drawn uniformly between 0.2- and 1.2-s postmovement onset. These parameters are chosen to explore reaching movements at a realistic spatial scale for humans while maintaining peak arm velocities that are comparable to those studied in related primate electrophysiology experiments of MI (Hatsopoulos et al. 2004; Li et al. 2001; Moran and Schwartz 1999).

How can the hybrid framework solve the switching-target-reach task? We begin by defining the continuous variable *x _{k}* as the arm state and the discrete variable

*s*as the target identity from a set of

_{k}*R*targets on a two-dimensional workspace (5) For this example, consider neural observations

*n*described as binned spikes of a point process (see

_{k}*Point process models of ensemble spiking activity*). The essential structure of this hybrid state space is depicted by the mixture-of-trajectories model (Fig. 1

*E*). To support switching targets, this static target diagram is altered by indexing the target

*s*with time as in the switching observation model (Fig. 1

*D*) (Wu et al. 2004).

Next we specify the conditional densities corresponding to each edge. Recall from methods that the density *p*(*s _{k}*

_{+1}|

*s*) is defined by a state transition matrix

_{k}**M**(6) This notation means that the entry in the

*i*th row and

*j*th column of

**M**corresponds to

*p*(

*s*

_{k}_{+1}=

*i*|

*s*=

_{k}*j*).

The density *p*(*x _{k}*

_{+1}|

*x*, s

_{k}_{k}

_{+1}) constrains the path of a reaching movement for any given target

*s*

_{k}_{+1}. This conditional density can be obtained from any of several reaching movement trajectory models (Kemere and Meng 2005; Kemere et al. 2003; Srinivasan et al. 2005; Srinivasan et al.; Yu et al. 2005b), directed specifically to the target location corresponding to

*s*

_{k}_{+1}. In this example, we use the standard (unaugmented) reach state equation detailed in (Srinivasan et al. 2005; Srinivasan et al.) based on the following free-movement state equation (7) The trajectory model described by the hybrid framework is not equivalent to the procedure used to simulate intended movements (described in the opening paragraph of this section). For example, whereas the hybrid framework assigns finite probability to the set of trajectories where there is more than one target switch, the reaching-movement simulation procedure allows only one switch that occurs randomly during the 0.2- to 1.2-s interval of the movement. Nevertheless, as demonstrated in the following text, the hybrid framework trajectory model of this problem will be a sufficiently appropriate description of the intended frequency and type of movement to allow neural observations to push the filter estimates along the simulated trajectory.

The point process observation model (*Eq. 3*) that describes *p*(*n*_{k}^{1:C}|*x _{k}*) is specified by a discrete-time conditional intensity function λ

_{k}

^{c}for each neuron

*c*. In this example, we choose a conditional intensity adapted from a model of primary motor cortex (Moran and Schwartz 1999) (8) (9) where

*v*and

_{x}*v*are the velocity components of

_{y}*x*in

_{k}*Eq. 5.*Here we assume that any lag between neural activity and the user's intentions is known and has been corrected to allow for the zero-lag indexing used in the preceding text. In practice, this lag can be estimated as another model parameter.

The parameters of the observation model *p*(*n*_{k}^{1:C}|*x _{k}*) can be tuned using point process adaptive filtering (Eden et al. 2004a) that also tracks changes due to neural plasticity. The parameters of the trajectory model in

*p*(

*x*

_{k}_{+1}|

*x*,

_{k}*s*

_{k}_{+1}) and

*p*(

*s*

_{k}_{+1}|

*s*) can be optimized a priori to reflect the types and frequency of behaviors that the neural prosthesis expects to support. Alternatively, adaptive methods will need to be developed to track the usage statistics of the device and adjust the trajectory model accordingly. In this example, we give our various competing filters an equal footing by providing them the actual trajectory and observation model parameters where applicable. A caveat is the state transition matrix

_{k}**M**. In free-movement estimation, this parameter is nonexistent. In mixture-of-trajectories estimation,

**M**is equal to the |

*S*| × |

*S*| identity matrix

**I**, where |

*S*| is the number of possible discrete states. For the hybrid framework, this parameter can be tuned to the expected frequency of target switches.

In these applications, we parameterized **M** as follows (10) where *b* = 1−*a*|*S*| − 1. This implies that transitions to states other than the current state are equiprobable. We found that performance was relatively insensitive to a range of choices between *a* = 0.9 and *a* = 0.99 but changed substantially for *a* = 1. (Although 0.99 is close to 1, this difference is geometrically magnified by successively multiplying 0.99 over multiple time steps.)

With the conditional densities specified, we can now use the hybrid point process filtering framework to drive the prosthetic device with ensemble spiking activity from motor cortex. We compare the performance of the hybrid framework against free-movement estimation and the mixture-of-trajectories model in a simulated analysis of the switching-target-reach task. The free-movement estimation procedure is implemented using a standard point process filter. This is mathematically equivalent to our hybrid framework where each target is given infinite uncertainty. The mixture-of-trajectories estimation method reported in Yu et al. (2007) is modified to implement the same reach state equation that our hybrid filter uses to provide equal grounds for comparison. This is mathematically equivalent to the hybrid filter with state transition matrix **M** = 1.

We also examine the effect of premovement instructed delay period activity that may be available to the prosthetic device. Such activity is known to provide information about the desired target in posterior parietal cortex (Andersen and Buneo 2002), premotor cortex (Weinrich and Wise 1982), frontal cortex (Schall 2001), and other brain regions. Premovement target information is easily incorporated into the mixture-of-trajectories model and hybrid filter by specifying a nonuniform initial posterior density on the target states *p*(*s*_{l}). We use one fixed moderately informative nonuniform posterior density (see parameter table) to simulate this premovement target information.

These filtering procedures were compared in a simulated version of the switching-target-reach task. The simulation comprised two stages. First, the subject's desired arm movement was generated based on the reach state equation (Srinivasan et al. 2005, 2006a) which is related to the stochastic optimal-control model (Bertsekas 2005; Kemere and Meng 2005; Todorov 2004). Optimal control has previously been described with relation to arm movement (Flash and Hogan 1985; Uno et al. 1989), including more recently, stochastic optimal control (Todorov and Jordan 2002). Second, the corresponding ensemble spiking activity from primary motor cortex (MI) was simulated based on *Eq. 8*, a velocity-tuned point process model of MI spiking activity (Moran and Schwartz 1999; Truccolo et al. 2005), using the time rescaling procedure that has previously been described in detail (Brown et al. 2002).

A summary of parameter choices is given in Table 1. The reach state equations require specification of the noise increment covariance and target uncertainty for each discrete target (Srinivasan et al. 2006). With states corresponding to position and velocity as described in *Eq. 5*, the original noise covariance of ε_{k} in *Eq. 7* was chosen to be nonzero in the entries corresponding to velocity increment variances (11) The uncertainty in each target state Π_{T} was also diagonal with (12) where *q*, *r*, and *p* are specified in Table 1.

For this example and application 3, movements and spikes were simulated at 0.01-s resolution. A 0.001-s resolution is typically used to simulate spikes with time rescaling, as with application 4. The coarser time resolution results in some distortion of the K-S plots based on time rescaling such as when multiple spikes are collapsed into one spike within a 10-ms interval, but speeds the simulation 10-fold.

The subject's arm movement was governed by the same unaugmented reach state equation (Srinivasan et al. 2006) used in the preceding text to define *p*(*x _{k}*

_{+1}|

*x*,

_{k}*s*

_{k}_{+1}) in the hybrid framework. Arm movement at any given time step followed the reach state equation corresponding to the current target with low target uncertainty (see parameter table). The constants

*q*,

*r*, and

*p*in that table refer to specific entries of the diagonal matrices for noise covariance and target uncertainty specified previously (Srinivasan et al. 2006). The target itself was allowed to switch once during the course of the movement. The target switch time was assigned at random, uniformly from a discrete set of possible times between 0.2 and 1.2 s postmovement onset, spaced at 0.2 s.

Because ensemble spiking is governed by conditional independence (see *Eq. 3*), the spiking activity of each cell could be generated separately. To generate the spike train of a given cell, the arm trajectory was first passed through the point process model in *Eq. 8.* The conditional intensity generated by the point process model of each neuron was then used to produce ensemble spiking activity based on the time rescaling theorem (Brown et al. 2002).

For each neuron *c*, model parameters β_{0}^{c} and β_{1}^{c} were chosen (see parameter table) to reflect typical background firing rate and depth of modulation for primate MI neurons during instructed-delay center-out reaching movements (Truccolo et al. 2005). The model parameter β_{2}^{c} was drawn randomly from a uniform distribution over [−π, π] to ensure that preferred directions were represented over all angles over the course of multiple trials. Note, however, that preferred directions were not necessarily evenly distributed on any given trial because of the random variation in drawing from the uniform distribution just once. Neurons in this simulated MI ensemble exhibited background firing rates of 10 spike/s and firing rates of 24.9 spike/s at a speed of 0.2 m/s in the preferred direction.

In total, five filtering procedures were compared in the simulated switching-target-reach task: free-movement estimation, mixture of trajectory estimation, and hybrid filtering, the last two methods being evaluated with and without premovement target information. Figures 2–5 provide a comprehensive view of the ability of these filtering procedures to convert MI spiking activity into reaching movements to switching targets. Figures 2 and 3 show sample trajectories driven by ensemble spiking activity under the various estimation procedures for a population of 25 neurons with a target switch at 1-s postmovement-onset. Figures 4 and 5 characterize how filter performance scales with ensemble size and target switch time for each of these procedures.

Sample decoding results from one trial without premovement target information (Fig. 2) show that the hybrid framework combines the strengths of free-movement estimation and the mixture-of-trajectories model. By incorporating target information, the hybrid framework and mixture-of-trajectories estimates drive the prosthetic arm to rest at the desired target location, whereas free-movement estimation leaves the arm displaced from the target and still moving at the 2-s mark. However, this same target information also causes the mixture-of-trajectories estimate to pull toward each passing target late in the reach (Fig. 2*A*). This “gravity effect” is reflected in the target probabilities under the mixture-of-trajectories model (Fig. 1*C*). In the second half of the reach, the current heading causes the passing targets (black and red lines) to quickly become highly likely, drawing the trajectory estimate toward those corresponding target locations. The hybrid framework overcomes this problem because it anticipates that targets may switch. By choosing *a* = 0.99 in the state transition matrix, the target densities (Fig. 1*D*) decay with time, and additional supporting neural activity is required to drive the probability of any given target to dominate the others. This mollifies the gravity effect of the mixture-of-trajectories model.

The hybrid filter also handles premovement target information differently from the mixture-of-trajectories model (Fig. 3). With the premovement information, the first target's probability under the mixture-of-trajectories model (blue line, Fig. 3*C*) approaches certainty faster than before (Fig. 2*C*). However, single trial decoding results (Fig. 3, *A* and *B*) show that the mixture-of-trajectories estimate appears to persist to the original target location even when the desired trajectory has begun to reorient to the new target. This is also seen in the target probabilities (Fig. 3*C*) where the first target (blue) dominates 200 ms beyond the time of the target switch. In contrast, the hybrid framework incorporates the target information early in the reach but progressively “forgets” or downweights its influence because it anticipates the possibility of a target switch, again by using *a* = 0.99. The free-movement estimate does not incorporate premovement information.

We next examined filter performance over a wide range of ensemble sizes, ranging from 15 to 80 neurons. Root mean squared (RMS) error was evaluated in two ways: averaged over the entire trajectory (Fig. 4, *A* and *B*) and over the endpoint at the 2-s mark (Fig. 4, *C* and *D*). Additionally, we examined the fidelity of position tracking (Fig. 4, *A* and *C*) and velocity tracking (*B* and *D*) separately. RMS error decreases for all five methods with larger ensemble sizes.

Trajectory RMS errors are typically smaller than endpoint RMS errors because trajectories begin with the accurate initial condition and accumulate error with time. All methods appear to perform equally well in endpoint error except free-movement estimation which does not incorporate the discrete target locations. Moreover, endpoint error appears to level out faster than trajectory error. This is likely due to the fact that just a few MI neurons are needed to make an accurate target classification, and once the accurate classification is made, the mixture-of-trajectories and hybrid framework methods will drive the prosthetic arm to rest at that target.

Premovement target information appears to provide a slight or insignificant improvement, but this is largely due to the moderate information provided by our choice of initial target prior. Higher fidelity premovement target information will likely make overshooting more pronounced in mixture-of-trajectories estimation and decrease RMS error in the first half of the reach generated by hybrid estimation.

Earlier target switches are easier to track for all methods than later target switches (Fig. 5) for a population of 25 neurons. Later switch times require faster velocity corrections, causing trajectory RMS error to rise across all methods (Figs. 5, *A* and *B*). Trajectory RMS error accumulates rapidly with later switch time for the mixture-of-trajectories model that lags in reorienting the arm movement unlike hybrid estimation, which anticipates switching and reorients quickly.

Endpoint errors (Fig. 5, *C* and *D*) under mixture-of-trajectories and hybrid estimation are largely insensitive to switch time in contrast to free-movement estimation. For mixture-of-trajectories and hybrid estimation, neural observations after the switch are sufficient to classify the target correctly, and because these latter methods incorporate the set of target locations, the prosthesis movement can reliably converge to the target.

Receiving information from the premovement activity that results in an incorrect maximum likelihood target classification is comparable to a zero second switch time because in both cases, premovement activity initially pushes path estimates toward the wrong final target. This represents the easiest case for tracking switching movements because subsequent neural activity over the full interval of reach time is available to correct estimates toward the final target. In later switches, shorter intervals of neural activity are available to redirect the arm movement. This is the regime where hybrid estimation shows marked improvement over the mixture-of-trajectories approach.

The simulation predicts that performance breaks down for all methods under moderate ensemble sizes for very late switches, where the target can no longer be reliably identified and high-velocity corrective movements must be tracked. A more subtle trend (Fig. 5, *A* and *B*) shows that free-movement estimation performs substantially worse than the mixture-of-trajectories model for early switch times but slightly better in very late switch times. These very late switch times make the overshoot and gravity effects of the mixture-of-trajectories model so pronounced that resulting trajectory estimates accumulate more RMS error than even the free-movement model.

### Application 2: EEG-based wheelchair navigation with definitive stopping

Motorized wheelchairs allow users with severe motor deficits to navigate within and outside the home. EEG has previously been suggested as one mechanism by which users could specify movements of their wheelchairs (Wolpaw et al. 2002). Although recordings from scalp-based EEG leads are susceptible to task-independent neural activity and a variety of artifacts, basic two dimensional cursor control was recently demonstrated using EEG in an on-screen center-out reaching task (Wolpaw and McFarland 2004). In the wheelchair application, the ability to reliably stop and start movements is critical to safety and functionality. Definitive stopping (and starting) is essential in a variety of practical settings, such as crossing streets, boarding a subway train, and navigating a grocery store aisle.

In this application, we simulated 10 user-intended consecutive point-to-point movements of a wheelchair in a 10 × 10-m workspace, punctuated by periods of rest (Fig. 6*B–G*). The intended movement kinematics were simulated based on a velocity-based linear quadratic control model used previously in the neural prosthetic device literature (Bertsekas 2005; Kemere and Meng 2005). Endpoints were drawn uniformly over the workspace, rest periods were drawn uniformly between 0 and 5 s, and travel times for each movement were calculated based on average velocities drawn uniformly between 0.5 and 2 m/s. Kinematics were simulated in discrete time increments of 100 ms.

Power in EEG-bands from 20 electrodes were simulated (Fig. 6*A*) at each time step *k* based on the user-intended wheelchair kinematics (13) This equation represents EEG-band energy as a linear function of user-intended wheelchair velocity, corrupted by zero-mean additive Gaussian noise ε_{k}. A very similar relationship was described in recent EEG (Wolpaw and McFarland 2004) and ECoG (Leuthardt et al. 2004) implementations. Although drawn independently at each time step, the additive noise is correlated across leads, as given by the following 20 × 20 noise covariance matrix (14) Three filters discussed in methods were designed to reconstruct intended wheelchair movement from EEG-band energy: free-movement estimation (Fig. 1*C*) (Wu et al. 2006), mixture of trajectory estimation (Fig. 1*E*), and the hybrid filter (Fig. 1*A*). With continuous-valued observations, free-movement estimation is equivalent to the standard Kalman filter (see Kailath et al. 2000 or appendix a), whereas mixture of trajectory estimation and hybrid filtering are akin to the interacting multiple model (IMM) framework (see Bar-Shalom et al. 2001) (see also methods). All three filters were provided exact parameters corresponding to the EEG model (*Eq. 13*) as equal grounds for comparison, although these parameters can be learned in practice with a training set based on standard approaches such as multiple linear regression (Kailath et al. 2000). Mixture of trajectory estimation and the hybrid filter require that we define discrete state variables, continuous state variables, and their interdependence at each time step *k*. The starting and stopping behaviors are naturally described as discrete states: *s _{k}* ∈ {stopped, moving}. The intended two dimensional position and velocity of the wheelchair corresponds to a continuous state

*x*, expressed as the kinematic vector on the right hand side of

_{k}*Eq. 13*in the preceding text. In both mixture of trajectory and hybrid filter implementations of the wheelchair system, the intended kinematics

*x*evolves depending on the discrete movement state

_{k}*s*(15) Here, each

_{k}*v*is a zero mean Gaussian random variable with the following covariance (16) The probability of being in either discrete state given the past discrete state,

_{k}*p*(

*s*

_{k}_{+1}|

*s*), evolves according to the following equation (17) In the mixture-of-trajectories model (Fig. 1

_{k}*E*), the discrete state is static, corresponding to a state transition matrix . In contrast, the hybrid filter accommodates nonzero probability of switching between moving and stopped states. Here, we use , although the exact choice of

**M**is in general based on the frequency with which states are expected to alternate. The probabilities of each discrete state are initialized uniformly at 0.5 in both filters.

Filters based on these design choices were implemented as described in methods. Filter reconstructions were graphed against the actual trajectory in position and velocity (Fig. 6, *B–G*). The free-movement and mixture-of-trajectory reconstructions almost completely overlap with each other (black) because the posterior probability of discrete states in the mixture-of-trajectories filter approaches unit probability in the *s _{k}* =

*moving*state within the first second of movement. Consequently, the free-movement and mixture-of-trajectory filters become virtually identical. For the hybrid filter, the probability of

*s*=

_{k}*moving*and

*s*=

_{k}*stopped*alternate in alignment with periods of intended movement and rest (Fig. 6

*H*). The posterior density still shows oscillations during the resting periods because the EEG-band signals during zero-intended velocity are consistent with both the movement (i.e., movement with 0 velocity) and resting dynamics (

*Eq. 15*).

Although very little error can be seen at the meter scale (Fig. 6, *B–G*), by magnifying the *x* velocity about four rest periods (Fig. 6*I*), we see that the free-movement and mixture-of-trajectory implementations exhibit a “resting tremor” of ∼10 cm/s amplitude, whereas the hybrid filter x-velocity remains closer to 1 cm/s. A normalized histogram of reconstructed speed during the simulated rest periods (Fig. 6*J*), shows that the hybrid filter is better able to hold still compared with the free-movement and mixture of trajectory filters. Importantly, definitive stopping does not compromise the ability to start movements from rest, or track movements during *s _{k}* =

*moving*periods (Fig. 6

*I*).

### Application 3: adaptive spike filtering under ongoing neuron death and discovery

In the current brain-driven interface prototypes, parameters of the observation model that relate neural activity to the user's intentions are learned by the algorithm during a training session consisting of real, imagined, or viewed actions (Carmena et al. 2003; Leuthardt et al. 2004; Musallam et al. 2004; Serruya et al. 2002; Shenoy et al. 2003; Taylor et al. 2002; Wahnoun et al. 2006; Wolpaw and McFarland 2004). In practice, this facilitates rapid (also termed “instant”) improvements in device performance because the device learns to match the user, whereas the user can maintain their operating assumptions about the device. Moreover, these parameters may need to be adjusted based on ongoing changes in brain function such as neuron death. Finally, parameters may need to be learned anew without additional training, such as when new neurons appear on recording electrodes. This can be achieved by employing observations of neural activity from existing neurons to learn parameters of the new neurons (Eden et al. 2004b). Neuron discovery arises in the use of adjustable-depth MEMS electrodes (Muthuswamy et al. 2005) and automated electrode positioning algorithms (Cham et al. 2005).

The concept of adaptive estimation is central to all of these parameter learning problems in neural prosthetic device applications. Although neural activity parameters are commonly fit in prosthesis applications using separate intervals of training data (Carmena et al. 2003; Hochberg et al. 2006; Leuthardt et al. 2004; Musallam et al. 2004; Santhanam et al. 2006; Serruya et al. 2002; Taylor et al. 2002; Wahnoun et al. 2006; Wolpaw and McFarland 2004), these parameters can also be adjusted on the fly by simultaneously estimating the user's intentions and neural activity parameters (Eden et al. 2004a; Eden et al. 2004b). In a more modular “lock-step” implementation, the user's intentions are estimated, followed by refinements of the neural activity parameter estimates, but these actions are alternated at every time step. This recursive procedure can be applied to neural activity in training periods, as well as afterward, in the face of neuronal death and discovery.

The hybrid framework, and our implementations of the free-movement and mixture-of-trajectories filters, are all compatible with recursive adaptive estimation procedures, such as the stochastic state point process filter designed to process spiking activity (Eden et al. 2004a) that was previously demonstrated in MI decoding of arm kinematics during neuron death and discovery (Eden et al. 2004b).

In application 3, we use MI-based coordination of reaching movements (application 1) to demonstrate that the hybrid framework is still consistent with this adaptive filtering framework under the “lockstep” procedure. This simply alters the existing implementation of each filter by adding a new point process filter (appendix a) (see also Eden et al. 2004a) that generates new estimates of MI tuning curve parameters at each time step. This filter uses current estimates of kinematics as if these estimates were exact [also called certainty equivalence (Bertsekas 2005)].

We now simulated a sequence of 300 intended reaching movements with target switching and repeated this simulation 50 times to construct percentiles on this performance (Fig. 7). Similarly, intended reaching movements were reconstructed with a population of 25 neurons. However, in this application, neurons were sequentially removed and replaced every 30 trials (once per minute of simulated movement time) with a new, randomly generated tuning curve unknown to the decoding algorithm. In response to death and discovery, the algorithm simply re-initializes its uncertainty (estimate covariance) of the new neuron's parameters to some large value. The initial estimates of the new neuron's parameters are assigned to the average antipreferred direction and average background firing rate of the estimated population tuning curves, which has been found to be an appropriate initialization in a previous study of adaptive filtering with MI neurons (Eden et al. 2004b). Following this initialization, the filter continues the lockstep learning process.

In total, 10 neurons were replaced over each run of 300 simulated movements. These conditions represent, in some ways, a more difficult scenario than encountered in practice: death and discovery are not expected to occur simultaneously, and this rate of death and discovery (1/min) far exceeds the anticipated rate. However, the tuning curve parameters were initialized to their exact values at the beginning of the simulation. This allowed us to dissociate performance degradation from initial parameter uncertainty versus ongoing neuronal death and discovery. In practice, these initial parameter estimates can be obtained (independent of the filter implementation) from training data to arbitrary precision (by using larger data sets) using standard parameter optimization procedures such as Matlab's glmfit routine if the model class is chosen to result in a convex parameter optimization (Truccolo et al. 2005). This is the case with the MI and EEG-band neural activity models described in applications 1 and 2, respectively.

Lockstep adaptive versions of all three filters were able to learn new neuron parameters and produce estimates of intended reaching kinematics throughout the simulation (Fig. 7). Estimated parameter values approximate actual values for the first replaced MI neuron in one 300-run sequence (Fig. 7, *B–D*), and the lockstep adaptive hybrid filter appears to converge faster and more reliably for α_{1} (Fig. 7*C*) and α_{2} (Fig. 7*D*) parameters in this example. Actual tuning curves over the entire population are shown for one run sequence at trial 1 (Fig. 7*G*) and trial 300 (Fig. 7*H*), where the first 10 neurons have been replaced with different, randomly drawn tuning curves. At trial 1, the tuning curve parameter estimates generated by the lockstep adaptive hybrid filter are initialized to equal the true simulated parameters (Fig. 7*I*). These tuning curve estimates still approximate the actual tuning curve at trial 300 (Fig. 7*J*), despite ongoing death and discovery of neurons in the recorded ensemble. The estimated tuning curve for neuron 10 at trial 300 (Fig. 7*J*) looks washed out. It has just been replaced on this trial, and its parameter estimates have been initialized to correspond to the anti-preferred direction of the estimated ensemble tuning curves, at the average background firing rate.

Sample trajectory decodes from the first (Fig. 7*E*) and last (*F*) trials of a 300-run sequence show qualitatively similar performance. However, the 90% confidence intervals in velocity trajectory RMS error (Fig. 7*A*, constructed with 5th and 95th percentiles over the 50 trials) were substantially improved for the hybrid filter versus the other two methods, on the order of 10 cm/s reduction in the 95th percentile. Moreover, the mixture-of-trajectories approach fared slightly worse in this measure than the free-movement approach, perhaps because of its pull toward passing targets (Fig. 7*F*) and as discussed under application 1 (Figs. 2 and 3). The filtering methods, especially the hybrid filter, show transient spikes in the RMS error at every 30 trials, corresponding to the time point at which a neuron has been replaced. These transient effects could potentially be removed by incorporating parameter uncertainty into the trajectory estimate, a consideration that was ignored by using certainty equivalence in the lockstep filtering procedure. Together, these results emphasize that adaptive estimation is compatible with all three filter types, and that hybrid filtering can facilitate more accurate on-line parameter estimation.

### Application 4: filtering under model mis-specification due to physiologically realistic history dependence in MI ensemble spiking

Inevitably, the mathematical models we use to describe the relationship between neural activity and the underlying intention will be incorrect [as from George Box (1979), “all models are wrong, but some are useful”]. Although application 3 showed that adaptive estimation could compensate unknown parameter assignments, we now quantify performance degradation in the face of model mis-specification, again in the context of MI spiking activity and reaching movements to switching targets.

Experimentally recorded spiking activity from intact primates can be filtered in an off-line fashion to examine filtering under model mis-specification in an attempt to predict closed-loop performance in a human patient population with representative movement disorders. In this application, we instead simulate spikes from empirical statistical models that have been fit to experimentally recorded data from intact primates. These history-dependent point process generalized linear models have recently applied to capture unprecedented realism in MI spike timing based on the Akaike Information Criterion, a model quality index that balances goodness-of-fit against model complexity to avoid overfitting (Truccolo et al. 2005). As a result, the spiking activity generated by these models is qualitatively and quantitatively similar to the experimentally recorded primate MI neurons from which they were derived. The advantage of applying empirical statistical models is that we can systematically evaluate the effect of including or removing structure in MI spiking (such as history dependence) on the performance degradation of competing algorithms that are not perfectly matched to this structure.

In this application, the model form used to generate spiking activity is different from the model form used to decode these spikes. Specifically, MI spiking activity is now simulated with physiologically realistic history dependence, in addition to modulation based on arm velocity. The strength of this history dependence is described in Fig. 8*A*, where the parameter values have been approximately extracted from Fig. 5*B* of Truccolo et al. (2005), which describes MI spiking activity recorded from a monkey performing a visuomotor pursuit-tracking task with a manipulandum that controls an on-screen cursor.

Simulated traces of kinematics (Fig. 8*B*, *1* and *2* from the top) and neural activity (Fig. 8*B*, *3–6* from the top) show that spiking activity that results from the same movement is qualitatively different under the physiologically realistic history-dependent model (Truccolo et al. 2005) (Fig. 8*B*, 5 and 6 from the top) versus the velocity modulated model (*Eq. 9*) used in applications 1 and 3 (Fig. 8*B*, 3 and 4 from the top). Note that the conditional intensity, or instantaneous firing rate, of the history-dependent model, well exceeds 200 spike/s, but these deflections are momentary, so that they generate physiologically observed bursts rather than producing an unrealistic, sustained elevation in mean firing rate.

How well do each of the three filters perform when they assume the basic velocity modulated model (*Eq. 9*) while they decode spiking activity that was actually generated with physiologically realistic history dependence (Truccolo et al. 2005)? We simulated 1,000 intended reaching movements to switching targets and reconstructed the intended movements with the three filter types to generate empirical probability densities on four performance metrics: position trajectory RMS error (Fig. 8*C*), position endpoint RMS error (Fig. 8*D*), velocity trajectory RMS error (Fig. 8*D*), and velocity endpoint RMS error (Fig. 8*E*). Each of these figures (8, *C* and *D*) includes panels for (from top to bottom) free-movement, mixture-of-trajectories, and hybrid filters. Each panel displays the distribution of performance when spikes are simulated with physiologically realistic history dependence (red) and with the original velocity modulated model (black). Position and velocity trajectory RMS errors degrade the least for the hybrid filter, followed by the mixture-of-trajectories, and free-movement filters. Both hybrid and mixture-of-trajectories filters maintain low endpoint RMS errors, while the free-movement filter shows relatively strong degradation in this performance measure.

Why does the hybrid filter show the least degradation under model mis-specification, followed by the mixture-of-trajectories and free-movement filters? The model mis-specification tends to push trajectory estimates away from their true values, whereas priors on the movement (embodied by the trajectory model) help restore estimates closer to their true values. Because the hybrid filter includes a switching discrete state, it is better able to describe the typical set of movements than the mixture-of-trajectories filter. Similarly, both the hybrid filter and mixture-of-trajectories filter constrain their priors on trajectory endpoint based on a known set of possible target locations. This not only improves endpoint RMS error but mitigates degradation under model mismatch that otherwise significantly affects the free-movement filter.

## DISCUSSION

We have introduced a unified approach for the design of filters for prosthetic devices. By using this technique, we can map spikes and continuous field potentials to estimates of the user's intention for a wide array of neural prosthetic device applications. The technique draws on Bayesian filter theory to generalize the dominant approaches to filter design in neural prosthetic devices (Brockwell et al. 2004; Carmena et al. 2003; Cowan and Taylor 2005; Eden et al. 2004b; Hochberg et al. 2006; Kemere and Meng 2005; Kemere et al. 2003; Musallam et al. 2004; Santhanam et al. 2006; Serruya et al. 2002; Shenoy et al. 2003; Srinivasan and Brown 2007; Srinivasan et al. 2006a; Taylor et al. 2002; Wessberg et al. 2000; Wolpaw and McFarland 2004; Wu et al. 2004, 2006; Yu et al. 2005, 2007).

The hybrid framework and filtering technique captures these previous methods completely but also allows for more flexibility. For this reason, device performance with this technique is expected to be equal or better than these previous methods. The versatility of this technique is illustrated in four applications using simulated ensemble spiking activity or multi-channel EEG: reaching to discrete targets that switch during movement, wheelchair navigation with definitive stopping, adaptive spike filtering under ongoing neuron death and discovery, and filtering under model mis-specification.

For both the hybrid point process filter and the IMM Switching Kalman Filter for Gaussian observation models, the number of operations at time step *k* scales with |*S _{k}*|, the number of values that a discrete state variable can take on. This is because the posterior density on the discrete state is nonparametric, and the posterior density on the continuous state is represented as a mixture of |

*S*| Gaussians. The particle filter (Brockwell et al. 2004; Doucet et al. 2001; Ergun et al. 2007), a Monte Carlo approach, would increase the fidelity the posterior density at the expense of increased computational cost. Ultimately, the way in which the posterior density is represented will depend on the cost of computation versus device performance in any specific application.

_{k}As shown in the previous section, the hybrid framework accommodates multiple discrete random processes by condensing them into one. Unfortunately, *n* discrete random variables, each with *p* possible values at step *k*, results in a condensed random variable with |*S*| = *n*^{p}. Fortunately, filtering on the hybrid framework can be parallelized fairly directly. This means that even with large |*S _{k}*|, the device can be controlled in real time if the hardware supports parallel computations. Parallelizing a digital hardware implementation may not necessarily save energy but could allow a slower clock speed. Parallelized algorithms do not necessitate the development of new types of processors [although power considerations are motivating new hardware architecture (Sarpeshkar et al. 2005)]; multithreading emulates parallelization on a single CPU, and even the prototyping language Matlab can now be parallelized with standard desktop computers (Kepner 2007). Formal analysis of hardware implementation is a key topic for subsequent investigation.

This real-time implementation is also possible because our algorithms are the same order of complexity as hybrid filtering based on the extended Kalman filter, which has been used extensively in mission-critical real-time military applications (Bar-Shalom et al. 2001). The algorithms are recursive rather than batch-processed, meaning that estimates of the user's intention can be updated using the latest neural observations without re-analyzing previously recorded activity.

In many applications, however, the number of discrete states can be kept small by using context. Context means that the space of device states is restricted at any given step in a way that still allows the user to eventually reach the desired device state. Consider how you organize files on your computer. By arranging your files in a sequence of subdirectories, you make it easy to scan through the list of files at each step. By placing all your files on the desktop, you are forced to select your file from a very large list, even though the file is just one mouse click away.

Looking forward, we expect to draw extensively on the rich field of dynamic Bayesian networks to address future applications. Prototyping is needed to determine the best computation/accuracy tradeoff for specific prosthetic devices. Learning and real-time sensory feedback (visual, somatosensory, auditory) must also be considered in developing algorithms that define the prosthetic interface. Associated technologies like computer vision and robotic control can be integrated with the hybrid framework to enhance real-world performance measures.

Ultimately, these approaches must be compared in a longitudinal study of device performance in activities of daily living for a population of individuals that experience well-defined stages and pathological mechanisms of a target motor deficit. Off- and on-line decoding studies in the laboratory will be important intermediate steps between simulated algorithm performance and end-stage clinical evaluation.

Finally, estimation with a minimum average cost criterion is not the only approach to formally describing the prosthetics problem. Future work will explore stochastic control, hierarchical design architectures, robustness, power consumption, hardware architecture, monetary cost, and other themes in systems design to achieve the type of performance in practical tasks that is necessary to benefit the full spectrum of limited motor function, from locked-in syndrome to single-arm amputation.

## APPENDIX

### Appendix A: Approximate point process filter for Gauss-Markov process (discrete time)

The Gaussian approximation to the posterior density with a Taylor series expansion about the prediction mean is employed in the following filter equations (Eden et al. 2004a). Consider a Gauss-Markov trajectory model (A1) A point process observation model is specified for an ensemble of *C* neurons. The conditional intensity function of the *c*th neuron, denoted λ_{k}^{c}, may depend on *H _{k}* and

*x*. For the

_{k}*k*th time step and

*c*th neuron,

*n*

_{k}

^{c}spikes arrive in a δ

_{k}time interval.

The prediction density mean *x _{k}*

_{+1|k}and covariance

*W*

_{k}_{+1|k}are (A2) (A3) The posterior density covariance

*W*

_{k}_{+1|k+1}and mean

*x*

_{k}_{+1|k+1}are (A4) (A5) Consider instead, an array of

*c*neural signals

*n*

_{k}

^{1:C}= [

*n*

_{k}

^{1},

*n*

_{k}

^{2}, …,

*n*

_{k}

^{C}]′ described by a Gaussian observation model (such as EEG) with mean

*D*

_{k}*x*+

_{k}*f*(s

_{k}_{k},

*H*) and variance

_{k}*R*. Here

_{k}*x*is a

_{k}*J*× 1 vector of continuous states,

*D*is a

_{k}*C*×

*J*matrix that may depend on

*s*, and

_{k}*f*(

_{k}*s*,

_{k}*H*) is a function that maps neural history to a

_{k}*C*× 1 vector, such as with ARMA models. The posterior density covariance and mean are then given by the standard Kalman filter equations (Kailath et al. 2000) (A6) (A7) The probability density in step 5 of the point process hybrid filter (see preceding text) is then replaced by (A8) to correspond to the Interacting Multiple Model (IMM) approach to switching Kalman filters (Bar-Shalom et al. 2001).

### Appendix B: Gaussian approximation to mixture of Gaussians

Consider a distribution composed of the weighted average of *R* multidimensional Gaussians (B1) with weights *d _{i}*, and where

*N*(

*x*; μ

_{i},

*W*) denotes the Gaussian probability density function with mean μ

_{i}_{i}and covariance

*W*.

_{i}The following standard approximation (Bar-Shalom et al. 2001) is obtained by moment matching [calculating the mean and covariance of *p*(*x*)] (B2) where (B3) (B4)

### Appendix C: Derivation of a point process hybrid filter to map spikes to hybrid prosthetic device states

For the *k*th discrete time step, define the user-intended continuous state *x _{k}*, discrete state

*s*, and the ensemble spiking activity of all

_{k}*C*neurons

*n*

_{k}

^{1:C}. The history of ensemble spiking at time step

*k*is given by

*H*= (

_{k}*n*

_{1}

^{1:C},

*n*

_{2}

^{1:C}, …,

*n*

_{k−1}

^{1:C}). Define the observation model

*p*(

*n*

_{k+1}

^{1:C}|

*x*

_{k}_{+1},

*s*

_{k}_{+1},

*H*

_{k}_{+1}) that represents the relationship between user intentions and spiking activity. Define the trajectory model

*p*(

*x*

_{k}_{+1}|

*x*,

_{k}*s*

_{k}_{+1}) and discrete state transition density

*p*(

*x*

_{k}_{+1}|

*x*,

_{k}*s*

_{k}_{+1}) that reflect the distribution of intentions that the user is expected to request over time.

In this section, we seek a recursive method to obtain *p*(*x _{k}*

_{+1},

*s*

_{k}_{+1}|

*n*

_{k+1}

^{1:C},

*H*

_{k}_{+1}) from

*p*(

*x*

_{k}*,s*|

_{k}*n*

_{k}

^{1:C},

*H*) and

_{k}*n*

_{k+1}

^{1:C}. This constitutes the point process hybrid filtering procedure.

For our specific hybrid state space in Fig. 1*A*, (C1) This implies that our problem is equivalent to obtaining *p*(*x _{k}*

_{+1}|

*s*

_{k}_{+1},

*n*

_{k+1}

^{1:C},

*H*

_{k}_{+1}) and

*p*(

*s*

_{k}_{+1}|

*n*

_{k+1}

^{1:C},

*H*

_{k}_{+1}) from

*p*(

*x*|

_{k}*s*,

_{k}*n*

_{k}

^{1:C},

*H*),

_{k}*p*(

*s*|

_{k}*n*

_{k}

^{1:C},

*H*), and

_{k}*n*

_{k+1}

^{1:C}.

Note that (C2) We now calculate *p*(*x _{k}*

_{+1}|

*s*

_{k}_{+1}

*n*

_{k+1}

^{1:C},

*H*

_{k}_{+1}) using

*Eqs. C3*–

*C8*and calculate

*p*(

*s*

_{k}_{+1}|

*n*

_{k+1}

^{1:C},

*H*) using

_{k}*Eq. C9.*

Observe that (C3) where *p*(*x _{k}*

_{+1}|

*s*

_{k}_{+1},

*H*

_{k}_{+1}) is the prediction density given by the Chapman-Kolmogorov equation (C4)

*Equations C3*and

*C4*comprise one step of a filter on

*p*(

*x*|

_{k}*s*

_{k}_{+1},

*H*

_{k}_{+1}) with the observation model

*p*(

*n*

_{k+1}

^{1:C}|

*x*

_{k}_{+1},

*s*

_{k}_{+1},

*H*

_{k}_{+1}) and trajectory model

*p*(

*x*

_{k}_{+1}|

*x*,

_{k}*s*

_{k}_{+1},

*H*

_{k}_{+1}) =

*p*(

*x*

_{k+1}|

*x*,

_{k}*s*

_{k+1}). For computational simplicity, we approximate both the trajectory model and posterior density

*p*(

*x*

_{k}_{+1}|

*s*

_{k}_{+1},

*n*

_{k+1}

^{1:C},

*H*

_{k}_{+1}) to be Gaussian. Such a filter (reproduced under appendix a) is developed in (Eden et al. 2004a) for point processes using a Taylor expansion about the prediction density mean rather than the posterior density mean employed in (Brown et al. 1998).

The density *p*(*x _{k}*|

*s*

_{k}_{+1},

*H*

_{k}_{+1}) is obtained by (C5) This density is a mixture of Gaussians that is approximated by one Gaussian density using a standard moment-matching formula given in appendix b.

The first density in the summation (*Eq. C5*) is calculated as follows (C6) where (C7) Here *p*(*s _{k}*

_{+1}|

*s*) is the discrete state transition density, and

_{k}*p*(

*s*|

_{k}*H*

_{k}_{+1}) is the posterior density on the discrete state, given in the previous iteration.

The second density in the summation (*Eq. C5*) is given by a quantity retained from the previous step (C8) This statement is verified in appendix d.

We now calculate *p*(*s _{k}*

_{+1}|

*n*

_{k+1}

^{1:C},

*H*

_{k}_{+1}) in

*Eq. C2*using the following relation (C9)

*Equation C7*calculates

*p*(

*s*

_{k}_{+1}|

*H*

_{k}_{+1}). The density

*p*(

*n*

_{k+1}

^{1:C}|

*s*

_{k}_{+1},

*H*

_{k}_{+1}) is given by the following integral. (C10) An approximation to this integral for point process observations is given by Laplace approximation as detailed in appendix e. Finally,

*p*(

*n*

_{k+1}

^{1:C}|

*H*

_{k}_{+1}) is a normalizing factor obtained by summing the numerator over all possible values of

*s*

_{k}_{+1}.

### Appendix D: Corollary

Verify *Eq. C8* that *p*(*x _{k}*|

*s*,

_{k}*s*

_{k}_{+1},

*H*

_{k}_{+1}) =

*p*(

*x*|

_{k}*s*,

_{k}*H*

_{k}_{+1}) (D1)

*From Fig. 1A*, observe that (D2) Thus

*Eqs. D1*and

*D2*imply that (D3)

### Appendix E: Laplace approximation of *p*(*n*_{k+ 1}^{1:C}|*s*_{k+1}, *H*_{k+1})

_{k+1}

This section derives the Laplace approximation of *Eq. C10*, repeated in the following text for convenience (E1) Define (E2) The Laplace approximation to *Eq. E1* is given by (E3) where the mode *x*_{0} maximizes *p*(*n*_{k+1}^{1:C}|*s _{k}*

_{+1},

*x*

_{k}_{+1},

*H*

_{k}_{+1})

*p*(

*x*

_{k}_{+1}|

*s*

_{k}_{+1},

*H*

_{k}_{+1}) for a given

*n*

_{k+1}

^{1:C}.

Approximate the mode as in (Eden et al. 2004a) using a prediction density, in this case given by (E4) Under this approximation, the following equalities hold (E5) (E6) where is precisely the variance of the Gaussian approximation to the posterior density given in (Eden et al. 2004a) and appendix a.

Using *Eqs. E5* and *E6*, the Laplace approximation (*Eq. E3*) simplifies to (E7) Express *p*(*n*_{k+1}^{1:C}|*x _{k}*

_{+1},

*s*

_{k}_{+1},

*H*

_{k}_{+1}) using a discrete-time approximation for point processes (Eden et al. 2004a) (E8) where λ

_{k+1}

^{c}denotes the conditional intensity of neuron

*c*at time step

*k*+ 1, which may be a function of

*s*

_{k}_{+1},

*x*

_{k}_{+1}, and

*H*

_{k}_{+1}.

Substituting this approximation into *Eq. E7*, we have the final approximate equation for *p*(*n*_{k+1}^{1:C}|*s _{k}*

_{+1},

*H*

_{k}_{+1}) (E9)

### Appendix F: Spike filtering with the hybrid framework: practical note on numerical issues

This section documents four points to consider when implementing the hybrid filter.

First, the spike filtering (hybrid point process) filter described in this paper uses a bank of stochastic state point process filters (SSPF), described in Eden et al. (2004a) and appendix a. As with the SSPF, the prediction or posterior covariance may become singular because of numerical implementation or badly conditioned if the values in certain matrix elements are dramatically smaller than others. In a practical implementation, it is useful to check that a covariance matrix is well-conditioned or invertible before taking the inverse operation required by the SSPF. If the posterior covariance is not invertible, perform a Fisher's scoring step instead of executing the posterior covariance equation, by removing the −(*n*_{k}^{c} ∂2 log λ*k*, *s**k*+1*c*∂*xk* term of the posterior covariance equation for just that time step. If the prediction covariance is badly conditioned, retain the prediction covariance as the posterior covariance.

Second, you may encounter divide-by-zero or floating-point errors if you incorrectly implement the nine step spike filtering procedure. Check that you are not dividing by a discrete state probability that has approached zero.

Third, to generate smoother continuous state trajectories, such as in application 1 of results, augment your state space to include acceleration terms, and introduce the nonzero diagonal term of increment covariance only in the acceleration dimensions.

Fourth, note that application 1 of results is a discrete-target version of problem of reaching to drifting targets that evolve over a continuum of positions. The discrete nature of the targets in application 1 necessitates the hybrid framework. Similarly, look for parallels between your application and discrete or continuous versions of it.

## GRANTS

This research was supported by the National Institutes of Health Medical Scientist Training Program Ruth Kirschstein National Research Service Award T32 GM-07753-28 to L. Srinivasan and Grant R01 DA-015644 to E. N. Brown and National Science Foundation Grant CCR-0325774 to S. K. Mitter.

## Footnotes

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

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

- Copyright © 2007 by the American Physiological Society