1 Introduction

Brain-machine interfaces (BMIs) are a powerful class of assistive devices that may one day restore movement ability to paralyzed individuals (Schwartz et al. 2006). These devices act by creating a direct mapping between recorded neural activity and the movement of an external actuator, like a computer cursor or a robotic arm (Fig. 1) (Chapin 2004; Serruya et al. 2002; Taylor et al. 2002; Carmena et al. 2003; Musallam et al. 2004). Early clinical trials with intracortical BMIs, which use as their control signal the recorded activity of populations of single neurons, have recently shown that paralyzed individuals can effectively control computer cursors (Hochberg et al. 2006) and robotic arms of varying complexity (Hochberg et al. 2012; Collinger et al. 2013; Aflalo et al. 2015). However, much work yet needs to be done to give subjects control over artificial limbs that might rival control of the natural limb (Gilja et al. 2012). The design of BMI control algorithms that might enable stable, robust, and naturalistic control is an active area of research (Gilja et al. 2011).

Fig. 1
figure 1

Schematic of a BMI as a feedback control system. The major components of a feedback control system (namely, the controller, control signals, plant, and feedback) are laid out on top of a typical BMI cursor control schematic, where the brain is identified as the controller, the control signals are neural activity (often tapped out of primary motor cortex), the plant is the combination of the BMI decoder and the cursor, and feedback is accomplished by watching the cursor movements

Full size image

A BMI decoding algorithm specifies how recorded signals (like recordings from intracortical multielectrode arrays) get translated into movement of the prosthesis. Currently, nearly all BMI decoding algorithms are designed from an estimation standpoint: it is assumed that neurons are tuned to different intended movements, and recorded firing rates are treated as noisy observations of that underlying motor intent. Thus, differences in BMI decoding algorithms are typically interpreted to result from the different assumptions that they make about how neurons represent motor intent. Algorithms that fall into this class include linear estimators such as the population vector algorithm (PVA) (Georgopoulos et al. 1986; Georgopoulos et al. 1988; Velliste et al. 2008) and optimal linear estimator (OLE) (Kass et al. 2005; Salinas and Abbott 1994; Collinger et al. 2013), state-space models such as the Kalman filter (Wu et al. 2006; Hochberg et al. 2012), Laplace-Gaussian filter (Koyama et al. 2010a; Koyama et al. 2010a), and the unscented Kalman filter (Li et al. 2009), as well as a host of related variants (eg, particle filter (Brockwell et al. 2004), SSKF (Malik et al. 2011), VBR (Li et al. 2011), ReFIT-KF (Gilja et al. 2012), PT (Zhang Y and Chase SM 2013), SDKF (Golub et al. 2014), and others). Differences in the performance of various algorithms are typically assumed to relate to how well the algorithms’ assumptions about motor intent encoding match the true underlying encoding of the neurons.

However, another way to interpret BMI system design is from a control system perspective. Each prosthetic effector, in conjunction with its decoding algorithm, acts as a control system that the subject needs to learn how to use through trial-and-error. Robotic arms are, by definition, physical control systems, while brain-controlled computer cursors may or may not be programmed to follow physical laws. When interpreted this way, differences in the performance of two algorithms might instead stem from differences in the control systems themselves: some control systems may be more usable than others, due to their physical characteristics or ease of conceptualization. Here we review some common BMI cursor decoding algorithms and derive the physical systems with which they correspond. We then re-interpret findings from the literature on which decoding algorithms work best in light of the physical systems that they represent. Intriguingly, when interpreted in this way, the literature suggests that BMI systems that follow physical laws are more usable than those that do not. Further, in on-line control it appears that BMI systems that reduce to equivalent physical forms tend to be equally-well controlled. These results have implications not only for BMI design, but may also shed light on the brain’s ability to conceptualize motor effectors of varying forms.

2 Linear physical systems with control signals

Before proposing our view of BMI design from a control system perspective, we first look at a general control system, as shown in Fig. 1. In this system, control signals generated by the controller (the brain) are used to drive the plant (the decoder and cursor). Feedback about the system is observed by the sensors (the eyes) and used by the controller to generate new control signals. If it is a physical system, the system dynamics, which dictate the plant’s evolution as driven by the control signal, should obey physical laws. Formally, the dynamics for a control system can be represented as

$$ d \boldsymbol{x} / d t = f(\boldsymbol{x}, \boldsymbol{u}, t), $$
(1)

where \(\boldsymbol {x} \in \mathbb {R}^{d}\) is the plant’s state in d dimensional space, \(\boldsymbol {u} \in \mathbb {R}^{n}\) is the control signal and t is the current time. To make the control system a legitimate physical system, f should obey physical laws, i.e., velocity should be the derivative of position, acceleration should be the derivative of velocity, etc.

For computational simplicity, a linear discretized approximation of Eq. (1) is typically used in practice. To do that, continuous time is discretized into a series of steps with the step size as a small time interval Δ (Δ = 1/60s for example), and the system dynamics at the t-th time step is expressed as

$$ \boldsymbol{x}_{t+1} = A_{t} \boldsymbol{x}_{t} + B_{t} \boldsymbol{u}_{t}. $$
(2)

Here, \(A_{t} \in \mathbb {R}^{d \times d}\) is a matrix that dictates the internal physics of the control system, and \(B_{t} \in \mathbb {R}^{d \times n}\) is a matrix that dictates how the control system responds to its control inputs.

If we are studying a physical system driving a prosthetic effector’s movement, then x is often the kinematics of the prosthetic effector, which could include the position, p, the velocity, v, and potentially higher-order terms as well. For example, the 1st order linear discretized approximation, where the state is the position x t = p t , has the form

$$ \boldsymbol{p}_{t+1} = A_{t} \boldsymbol{p}_{t} + B_{t} \boldsymbol{u}_{t}. $$
(3)

To ensure that position is the integral of velocity, we need to have the following relationship:

$$ \left( \begin{array}{cc} \boldsymbol{p}_{t+1} \\ \boldsymbol{v}_{t+1} \end{array}\right) = \left( \begin{array}{cc} I & \varDelta I \\ \alpha_{t} & \beta_{t} \end{array}\right) \left( \begin{array}{c} \boldsymbol{p}_{t} \\ \boldsymbol{v}_{t} \end{array}\right) + \left( \begin{array}{cc} \boldsymbol{0} \\ B_{t} \end{array}\right) \boldsymbol{u}_{t}. $$
(4)

Note that the control signal affects the next velocity, only.

To gain some insight into the meaning of these control system parameters, consider a simple 2nd order physical system (Fig. 2). This system consists of a point mass, m p , attached to a viscous damper and an elastic spring connected in parallel. The forces operating on the point mass come from external sources, F e x , as well as internal from the spring and damper. The force from the spring is F s = −k p (where k is the spring constant and p is the position of the point mass, measured relative to the equillibrium position of the spring), and the force from the damper is F d = −η v (where η is the damping coefficient and v is the velocity of the point mass). Because these elements act in parallel, these two forces sum at the point mass, to create a total internal force F i n , where

$$ F_{in} = F_{d} + F_{s} = - \eta v - k p. $$
(5)

Together with the external force the acceleration of the point mass may be derived from Newton’s laws of motion:

$$ a = F_{tot} / m_{p} = (F_{in} + F_{ex}) / m_{p} = (- k p - \eta v + F_{ex}) / m_{p}. $$
(6)

Discretized by a small time duration Δ, the updated velocity is given as

$$ v_{t+1} = v_{t} + \varDelta a_{t} = v_{t} + \varDelta (- k p_{t} -\eta v_{t} + F_{ex,t}) / m_{p}. $$
(7)

and the updated position is given as

$$ p_{t+1} = p_{t} + \varDelta v_{t}, $$
(8)

Putting this all together into matrix form, the discretized control system for the Kelvin-Voigt model in Fig. 2 is seen to be:

$$ \left( \begin{array}{c} p_{t+1} \\ v_{t+1} \end{array}\right) = \left( \begin{array}{cc} 1 & \varDelta \\ -\varDelta k/m_{p} & 1 - \varDelta \eta/m_{p} \end{array}\right) \left( \begin{array}{c} p_{t} \\ v_{t} \end{array}\right) + \left( \begin{array}{c} 0 \\ \varDelta/m_{p} \end{array}\right) F_{ex,t}. $$
(9)
Fig. 2
figure 2

A simple 2nd order physical control system. Here, a point mass m is hooked up to a parallel combination of a spring (with spring constant k) and a damper (with damping coefficient η). This configuration corresponds to the well-known Kelvin-Voigt model from material science

Full size image

Comparing Eq. (9) with the general 2nd order physical control system presented in Eq. (4) reveals that the external force F e x,t plays the role of the control signal u t , the spring-like elastic effects account for α t , and the viscous damping effects account for β t .

3 BMI systems from a physical control systems perspective

In this section we re-interpret some fairly common decoding algorithms, including the PVA, OLE and Kalman filter, in terms of the physical control systems to which they correspond. Figure 1 casts the general BMI system into the control system framework. Here the motor cortex drives the BMI system by generating a proper set of neural activation patterns, which are the end result of a sequence of brain computations that take both visual feedback and goal information into account. Visual feedback of the prosthetic effector’s movement is used to correct future movement. Therefore, BMI system design is essentially a control system design problem where one attempts to construct a BMI system that is as usable as possible. Formally, the dynamics of a BMI system can be represented as

$$ d \boldsymbol{x} / d t = g(\boldsymbol{x}, \boldsymbol{y}, t). $$
(10)

where \(\boldsymbol {y} \in \mathbb {R}^{n}\) is the recorded firing rates from a population of n neurons. The only difference between Eqs. (10) and (1) is that the control signal u takes the form of neuronal firing rates y.

When discussing the physical implementation of different BMI decoders, it will be necessary to distinguish the following three variables: motor intent, estimated motor intent, and implemented movement. Note that here and elsewhere in this document, we use the star notation to refer to intended movements (e.g., the intended velocity will be denoted as v ∗). Estimates of those intents will be indicated by an overhead carat (e.g., the estimated velocity will be denoted \(\hat {\boldsymbol {v}}\)). Finally, the movement that is implemented by the prosthesis will come without any notation (e.g., the implemented velocity will be denoted as v). Please refer to Table 1 for all the notations.

Table 1 Notations
Full size table

3.1 Linear estimators

We start with one of the earliest decoding algorithms, the population vector algorithm (PVA), first introduced by Georgopoulos and colleagues in 1986 (Georgopoulos et al. 1986). The central assumption of the PVA is that neurons are tuned to the direction of desired movement. Formally, with this assumption the firing rate of the neuron, y, can be written as

$$ y = b_{0} + b_{x} d_{x}^{\ast} + b_{y} d_{y}^{\ast} + \varepsilon $$
(11)

where b 0, b x , and b y are linear regression coefficients, \((d_{x}^{\ast }, d_{y}^{\ast })^{T}\) is a unit vector that points in the direction of the intended movement, and the noise ε is assumed to be Gaussian. For simplicity, we have written this encoding model in two dimensions, though it has been found to generalize to three dimensions as well (Georgopoulos et al. 1986). This linear regression reduces to the well known “cosine-tuning” function (Georgopoulos et al. 1982).

If we allow m to be the magnitude of the regression coefficient vector (b x ,b y )T (also known as the “modulation depth”), and we denote the “preferred direction” of the neuron as d = (b x ,b y )T/m, we can write down the decoding equation for the PVA as:

$$ \hat{\boldsymbol{v}} = \frac{k_{s}}{n} \sum\limits_{i = 1}^{n} \frac{y_{i} - b_{0,i}}{m_{i}} \boldsymbol{d}_{i}. $$
(12)

Here, i indexes each of the n recorded neurons and k s is a constant that scales the unitless decoded direction into a velocity (Chase et al. 2012). The general interpretation of this equation is that each neuron “pushes” the cursor in its preferred direction, with the amount of the push being proportional to its firing rate. The sum of all of these pushing vectors is decoded as the population vector, which gets implemented in the BMI system. In this document, for simplicity, we assume the firing rates are alway centralized so we replace (y i −b 0,i ) with y i in the following sections. We can rewrite Eq. (12) in matrix form by gathering the preferred directions into a single matrix D of size 2×n, where each column corresponds to the preferred direction of a single neuron, and gathering the modulation depths m i into a single diagonal matrix M:

$$ \hat{\boldsymbol{v}} = (k_{s}/n) D M^{-1} \boldsymbol{y}. $$
(13)

To implement the PVA estimate of velocity into an actual device, the velocity of the prosthesis is set equal to its estimate, i.e., \(\boldsymbol {v}_{t} = \hat {\boldsymbol {v}}_{t}\). The implemented position is then set equal to the integral of these velocity commands, p t+1 = p t +Δ v t . Therefore, to represent the decoding from the physical system perspective, we have

$$ \text{PVA physical system: } \boldsymbol{p}_{t+1} = \boldsymbol{p}_{t} + \varDelta (k_{s}/n) D M^{-1} \boldsymbol{y}_{t}, $$
(14)

which is a special case of a 1st order linear physical control model (3), where A t = I and B t = Δ(k s /n)D M −1.

The PVA is a biologically-inspired algorithm. In practice, however, it has been shown to return biased estimates of motor intent when the preferred directions of the recorded neurons are not uniformly distributed (Salinas and Abbott 1994; Kass et al. 2005; Chase et al. 2009). To compensate for this bias, the optimal linear estimator (OLE) has been proposed. As the name implies, the OLE computes the optimal linear estimate according to square errors. The encoding model is written as

$$ \boldsymbol{y} = B \boldsymbol{v}^{\ast} + \boldsymbol{\varepsilon}, $$
(15)

and the decoding model is

$$ \hat{\boldsymbol{v}} = (B^{T} B)^{-1} B^{T} \boldsymbol{y}, $$
(16)

We can see that if the preferred directions of recorded neurons are uniformly distributed, then B T B∝I and OLE is equivalent to PVA. Similar to PVA, the implemented velocity is also set equal to the estimated velocity in the decoding model of OLE and the position of the prosthesis is derived by integrating velocity. Thus, the physical system corresponding to OLE decoding is

$$ \text{OLE physical system: } \boldsymbol{p}_{t+1} = \boldsymbol{p}_{t} + \varDelta (B^{T} B)^{-1} B^{T} \boldsymbol{y}_{t}. $$
(17)

This is again another special case of a 1st order linear physical control system (Table 2 and 3).

Table 2 BMI decoders under physical system perspective
Full size table
Table 3 BMI decoders comparison
Full size table

Both the PVA and the OLE correspond to first-order physical control systems, albeit with slightly different mappings from neural firing to cursor movement. Experimental results from (Salinas and Abbott 1994; Kass et al. 2005; Chase et al. 2009) have shown that the OLE overcomes that bias. From an estimation standpoint, the OLE should be a better decoding algorithm than the PVA.

However, Chase and colleagues (Chase et al. 2009) also demonstrated that the PVA and the OLE perform equivalently on-line: subjects were just as adept at controlling the PVA as they were at controlling the OLE, despite the fact that neural activity was mapped to different cursor movements under the two algorithms. One possible interpretation of this is that subjects learn the mapping from neural activity to decoder, be it a biased decoder like the PVA or an unbiased decoder like the OLE. The difference between the mapping from neural activity to cursor movement produced under these two decoders is akin to a visuomotor distortion, and visuomotor distortions are learned very quickly (Krakauer et al. 2000; Paz et al. 2005; Wu and Smith 2013). From this standpoint, the neurons are rapidly changing their activity to provide appropriate control signals to the device. Once learning is accomplished, both the OLE and the PVA give the subject a 1st order linear physical control system to control, and there appears to be no significant difference between the usability of these systems. This learning process is somewhat intuitive. Experiments demonstrate subjects can even learn the shuffled decoder with no predictive power after several days of practice (Ganguly and Carmena 2009; 2010).

3.2 Linear state-space decoders

A problem with the linear estimators described in the previous section is smoothing: when implemented on small time bins Δ, the movement estimates can be quite noisy. To compensate for this it is common to smooth the firing rate estimates (Koyama et al. 2010a). Naturally, any smoothing also affects the physical system.

Linear state space models handle smoothing in a more elegant manner by applying a smooth prior to the evolution of the intended kinematics. This prior takes the form of a linear dynamics equation:

$$ \boldsymbol{{x}}^{\ast}_{t} = A \boldsymbol{{x}}^{\ast}_{t-1} + \boldsymbol{\omega}_{t}, $$
(18)

where \(\boldsymbol {\omega }_{t} \sim \mathcal {N}(\mathbf {0}, R)\) is assumed to be zero mean Gaussian noise. The observations (recorded spike counts) are assumed to be linearly related to the kinematics:

$$ \boldsymbol{{y}}_{t} = B \boldsymbol{{x}}^{\ast}_{t} + \boldsymbol{\varepsilon}_{t}, $$
(19)

where again the noise in the fit is assumed to be Gaussian, \(\boldsymbol {\varepsilon }_{t} \sim \mathcal {N}(\mathbf {0}, Q)\).

In this implementation, the decoding problem is a Bayesian inference problem where the goal is to estimate the a posteriori probability of the intended kinematics when given the recorded spike counts. The solution in this case is the well-known Kalman filter (Kalman 1960), first introduced for BMI use by Wu et al. (Wu et al. 2003). The Kalman filter provides an efficient recursive algorithm to compute the posterior probability \(p(\boldsymbol {{x}}^{\ast }_{t} | \boldsymbol {{y}}_{1,\ldots ,t})\). The optimal estimate \(\hat {\boldsymbol {{x}}}\) of the intended kinematics x ∗ is the mean of this distribution, which can be estimated through the following set of recursive equations:

$$ K_{t} \,=\, (A \hat{\Sigma}_{t-1}A^{T} \!+ R) B^{T} \!(B (A \hat{\Sigma}_{t-1}A^{T} \!+ R) B^{T} \!+ Q)^{-1}\!, {} $$
(20)
$$ \hat{\boldsymbol{{x}}}_{t} = (A - K_{t} B A) \hat{\boldsymbol{{x}}}_{t-1} + K_{t} \boldsymbol{{y}}_{t}, $$
(21)
$$ \hat{\Sigma}_{t} = (I - K_{t}B) (A \hat{\Sigma}_{t-1}A^{T} + R). $$
(22)

Here K t is the familiar Kalman gain which computes the optimal mixing between reliance on information from the prior and information from the observations according to the noise assumed to be in each. \(\hat {\varSigma }_{t}\) is the estimate of the covariance of \(p(\boldsymbol {{x}}^{\ast }_{t} | \boldsymbol {{y}}_{1,\ldots ,t})\).

To use the Kalman filter, an initial covariance matrix \(\hat {\varSigma }_{0}\) and kinematics \(\hat {\boldsymbol {{x}}}_{0}\) must be given. In practice, it is common to center the prosthesis so all terms in \(\hat {\boldsymbol {{x}}}_{0} = \boldsymbol {{0}}\) with certainty \(\hat {\varSigma }_{0} = 0\). The Kalman gain is time dependent. However, it tends to converge to a stable value within a few timesteps, independent of the observations y (Chui and Chen 2009; Malik et al. 2011). Another common practice is to initialize the system with \(\hat {\boldsymbol {{x}}}_{0}\) and \(\hat {\varSigma }_{0}\) as the matrix that ensures that K t is stable (Dethier et al. 2011; Sadtler et al. 2014).

The relationship between the kinematics x t that are actually implemented and the estimated kinematics \(\hat {\boldsymbol {{x}}}_{t}\) depends on exactly how the Kalman filter is implemented. Here we introduce two popular implementations, the velocity Kalman filter (VKF) (Kim et al. 2008; Hochberg et al. 2012) and the position-velocity Kalman filter (PVKF) (Wu et al. 2006; Homer et al. 2013; Gowda et al. 2014).

3.2.1 Velocity Kalman filter

The VKF assumes that neurons are tuned to the intended velocity. Thus, the state evolution equation is

$$ \boldsymbol{{v}}^{\ast}_{t} = A \boldsymbol{{v}}^{*}_{t-1} + \boldsymbol{\omega}_{t} $$
(23)

and the observation equation is

$$ \boldsymbol{{y}}_{t} = B \boldsymbol{{v}}^{\ast}_{t} + \boldsymbol{\varepsilon}_{t}. $$
(24)

From Eq. (21), we have the estimated intended velocity as

$$ \hat{\boldsymbol{{v}}}_{t} = (A - K_{t} B A) \hat{\boldsymbol{{v}}}_{t-1} + K_{t} \boldsymbol{{y}}_{t}. $$
(25)

The implemented velocity v t in this case is set equal to \(\hat {\boldsymbol {{v}}}_{t}\) and the position is the integral of the velocity (Kim et al. 2008; Hochberg et al. 2012). In matrix form, the implemented movement can be written as

$$\begin{array}{@{}rcl@{}} &&\text{VKF physical system: }\\ &&\left( \begin{array}{c} \boldsymbol{{p}}_{t} \\ \boldsymbol{{v}}_{t} \end{array}\right) = \left( \begin{array}{cc} I & \varDelta I \\ 0 & A-K_{t} B A \end{array}\right) \left( \begin{array}{cc} \boldsymbol{{p}}_{t-1} \\ \boldsymbol{{v}}_{t-1} \end{array}\right) + \left( \begin{array}{c} \mathbf{0} \\ K_{t} \end{array}\right) \boldsymbol{{y}}_{t}. \end{array} $$
(26)

Comparing with Eq. (4), we can see that the VKF is a 2nd order linear physical control system where the system state includes position and velocity, x t = (p t ,v t )T, with an elastic term, α t , that is equal to zero and a viscous term, β t = A−K t B A. It is interesting to note the parallel between force and velocity representations that emerge in this implementation of the Kalman filter. Even though the VKF makes the assumption that neurons are driven by intended velocities, they play the role in Eq. (26) of providing a force input to the system.

To our knowledge, nobody has directly compared the VKF to either the OLE or the PVA in on-line control. However, Koyama and colleagues demonstrated that a variant of the VKF called the Laplace-Gaussian filter (LGF) outperformed the PVA and OLE on-line (Koyama et al. 2010a). The LGF and the VKF are both state-space decoders, and differ in only two main respects: the LGF fits an observation model that assumes Poisson noise statistics and log-linear tuning to intended velocity, whereas the VKF assumes linear Gaussian tuning of neurons to intended velocity. Importantly, this implies that the physical control system representing the VKF and the LGF will be of the same form: second order with no elastic terms.

Koyama and colleagues performed simulations and off-line trajectory reconstructions to determine that the key factor that allowed the LGF to outperform the PVA and OLE was its state-space formulation: they found no significant differences in the performance of the LGF relative to the VKF. We therefore take this as indirect evidence that the VKF would outperform the OLE and the PVA on-line. Why should this be the case? From an estimation standpoint, the interpretation would be that intended velocities really do evolve smoothly over time as implied by Eq. (23), and so incorporating the fact enables better estimates of the velocity intent. However another interpretation is that the VKF and OLE are fundamentally different control systems, and 2nd order physical control systems may simply be easier to control than 1st order physical control systems.

3.2.2 Position velocity Kalman filter

Another widely used Kalman filter model is the PVKF. In contrast to the VKF, the PVKF assumes that neurons are tuned to both the intended position and intended velocity. Thus, the state is \(\boldsymbol {{x}}^{\ast }_{t} = (\boldsymbol {{p}}^{\ast }_{t}, \boldsymbol {{v}}^{\ast }_{t})^{T}\) and to encourage the state evolution model to obey physical laws, it is typically set to be

$$ \left( \begin{array}{cc} \boldsymbol{{p}}^{\ast}_{t} \\ \boldsymbol{{v}}^{\ast}_{t} \end{array}\right) = \left( \begin{array}{cc} I & \varDelta I \\ \mathbf{0} & A_{v} \end{array}\right) \left( \begin{array}{cc} \boldsymbol{{p}}^{\ast}_{t-1} \\ \boldsymbol{{v}}^{\ast}_{t-1} \end{array}\right) + \left( \begin{array}{cc} \mathbf{0} \\ \boldsymbol{\omega}_{v,t} \end{array}\right). $$
(27)

The PVKF observation model is

$$ \boldsymbol{{y}}_{t} = (B_{p}, B_{v}) \left( \begin{array}{cc} \boldsymbol{{p}}^{\ast}_{t} \\ \boldsymbol{{v}}^{\ast}_{t} \end{array}\right) + \boldsymbol{\varepsilon}_{t}. $$
(28)

From Eq. (20) we can compute K t . Dissociating K t into two parts corresponding to position and velocity as K t = (K p,t ,K v,t )T, we can write the estimated intended position and velocity as

$$ \left( \begin{array}{c} \hat{\boldsymbol{{p}}}_{t} \\ \hat{\boldsymbol{{v}}}_{t} \end{array}\right) = \left( \begin{array}{ccc} I - K_{p,t} B_{p} & \varDelta I - K_{p,t} H \\ -K_{v,t} B_{p} & A_{v} - K_{v,t} H \end{array}\right) \left( \begin{array}{cc} \hat{\boldsymbol{{p}}}_{t-1} \\ \hat{\boldsymbol{{v}}}_{t-1} \end{array}\right) + K_{t} \boldsymbol{{y}}_{t}, $$
(29)

where H = Δ B p +B v A v

In the PVKF, the estimated position is not, in general, equal to the integral of the estimated velocity. Even though the state evolution equation (27) biases estimates of position and velocity to obey this rule, it is not a hard constraint: a compromise between position and velocity will be estimated that best explains the observed spike rates. This then leaves one with a choice when trying to implement the PVKF, since both the position and velocity estimates cannot be simultaneously implemented. One common method is to use the estimated position as the implemented position, and to allow the implemented velocity to evolve as \(\boldsymbol {{v}}_{t} = (\hat {\boldsymbol {{p}}}_{t+1} - \hat {\boldsymbol {{p}}}_{t}) / \varDelta \) (Wu et al. 2006):PVKF, position implementation:(Not a simple physical control system)

$$ \left( \begin{array}{c} \boldsymbol{{p}}_{t} \\ \hat{\boldsymbol{{v}}}_{t} \end{array}\right) = \left( \begin{array}{cccc} I - K_{p,t} B_{p} & \varDelta I - K_{p,t} H \\ -K_{v,t} B_{p} & A_{v} - K_{v,t} H \end{array}\right) \left( \begin{array}{cc} \boldsymbol{{p}}_{t-1} \\ \hat{\boldsymbol{{v}}}_{t-1} \end{array}\right) + K_{t} \boldsymbol{{y}}_{t}. $$
(30)

With this implementation, the estimated velocity, \(\hat {\boldsymbol {{v}}}_{t}\), becomes a latent variable that keeps the system from having a simple physical control system correlate. Experimental results from (Kim et al. 2008) show that when the PVKF is implemented in this way, its performance is inferior to the performance of the velocity-only Kalman filter. This fact is hard to rationalize from the estimation viewpoint, since the encoding model assumed by the PVKF typically fits the firing rates better than the encoding model assumed by the VKF. From the control system viewpoint, however, a possible explanation of this result is that the PVKF with position implementation is not a simple physical control system, and is therefore difficult to use.

There are other ways to implement the PVKF. One way is to use the estimated velocity as the implemented velocity and treat the estimated position \(\hat {\boldsymbol {{p}}}_{t}\) as the hidden variable. Another is to use a linear combination of estimated velocity and estimated position as the implemented velocity. This latter method was shown in (Homer et al. 2013) to work better than the position-implementation of the PVKF. However, none of these versions give the subject a simple physical system without hidden states to learn to control.

There is one implementation of the PVKF, however, that does correspond to a simple physical control system with no hidden states. To do this, the implemented velocity is made equal to the estimated velocity, and the estimated position becomes the integral of the estimated velocity:PVKF, velocity implementation:(Physical control system)

$$ \left( \begin{array}{c} \boldsymbol{{p}}_{t} \\ \boldsymbol{{v}}_{t} \end{array}\right) = \left( \begin{array}{cc} I & {\Delta} I \\ -K_{v,t} B_{p} & A_{v} - K_{v,t} H \end{array}\right) \left( \begin{array}{c} \boldsymbol{{p}}_{t-1} \\ \boldsymbol{{v}}_{t-1} \end{array}\right) + \left( \begin{array}{c} \mathbf{0} \\ K_{v,t} \end{array}\right) \boldsymbol{{y}}_{t}. $$
(31)

As it turns out, this implementation of the PVKF has been used by Gilja and colleagues to achieve the best BMI control demonstrated to date (Gilja et al. 2012). Their implementation of this equation was one of two design innovations of their ReFIT-KF algorithm, the other being a new method for re-calibrating the device from on-line training data. Although they derived (31) in a different way, using a causal-intervention step that forces the variance of the position estimate to go to zero, the net effect is the same, and leads to a second order physical system with both elastic and damping terms.

The velocity-implementation of the PVKF has been shown to handily outperform the VKF (Gilja et al. 2012). One interpretation of why this is true is that the causal intervention step better captures the fact that subjects know the position of the cursor from sensory feedback, so there should be no uncertainty about it. While this clearly cannot be exactly true, due to sensory feedback delays and sensory noise (Golub et al. 2013), it may be true enough to allow for better control. Another interpretation might be that neurons are driven by a non-volitional, positional signal representing the real position of the cursor, and this equation allows that “nuisance variable” term to be removed. However, yet another interpretation might be that this implementation of the PVKF results in a simple, 2nd order physical control system that still allows for neurons to modulate as a function of position. It may further be the case that 2nd order physical control systems that include a certain elastic component (a non-zero α t in Eq. (4)) are more usable than those systems, like the VKF, that do not have that component.

4 Discussion

The optimal way to implement a BMI decoding algorithm is an important question relevant to clinical deployment of neural prostheses. Here we recast this problem from the perspective of control system design, and derive the physical control systems corresponding to various types of decoders commonly used in BMI cursor control. This process enables new insights into BMI design, and suggests novel explanations about why some decoders have been shown to perform better than others. In particular, the literature suggests that: 1) 2nd order physical systems tend to be more usable than 1st order physical systems, 2) decoders that cannot be expressed as simple physical control systems do not appear to work as well as those that can be expressed this way, and 3) a 2nd order control system with elastic terms seems to work better than one without.

Recent work has highlighted the utility of approaching BMI design as a separate problem from inferring natural behavior (Tillery and Taylor 2004; Chase and Schwartz 2010). Marathe and Taylor (Marathe and Taylor 2011) studied the effect of mapping one control parameter (e.g., position, velocity, or goal) to the control of another. They found that the optimal mapping was not necessarily one-to-one, but rather changed as a function of different types of decoding noise. Gowda and colleagues (Gowda et al. 2014) have presented a thorough investigation of the dynamical systems properties of the PVKF. They found that certain implementations of the decoder could create workspace attractor points that might be detrimental to BMI control. These studies point out the gains that may be realized when BMI control is not constrained to reflect the neural encoding of natural arm dynamics, and emphasize the importance of the physical control system perspective when interpreting BMI performance.

4.1 Suggestions for new BMI systems

Our analyses suggest a number of new approaches to BMI decoding algorithm design that might prove fruitful. Of all the decoding algorithms we reviewed, none went beyond 2nd order. Given that 2nd order systems appear more usable than 1st order systems, it is interesting to speculate as to whether a 3rd or 4th order system would be even easier to control. These higher-order systems may actually be a closer match to the human arm: in (Liu and Todorov 2007), Liu and colleagues model the arm as a 3rd order linear physical system and are able to capture many of the emergent features of natural reaching movements. There have been instances in the literature that have included acceleration and higher order terms in their decoding algorithms. For example, Wu and colleagues compared the performance of a Kalman filter decoder with up to 6th order terms, and found that the 3rd order model consisting of position, velocity, and acceleration terms provided the best performance in their off-line trajectory reconstruction (Wu et al. 2006). However, they used the position-implementation of their decoder, which we have already demonstrated does not correspond to a simple physical system beyond 1st order. It would be interesting to test how physical implementations of higher order control systems might perform on-line.

Another fruitful approach might be the design of state-dependent control systems. The PVA and the OLE are both static systems, i.e., the physical control parameters do not vary with time. The Kalman filter technically has time varying parameters, but in practice the Kalman gain converges to a constant within a few timesteps, and some researchers even initialize it at the converged-values to keep it time-invariant (Malik et al. 2011; Sadtler et al. 2014). On the other hand, state-dependent control is also studied to deal with the situations where static system is no longer appropriate. One such situation is the instability of tuning curves between sessions (Rokni et al. 2007; Chestek et al. 2009). To compensate for such instability, different adaptive decoding algorithms are proposed where decoding parameters are updated iteratively over time (Li et al. 2011; Shpigelman et al. 2009; Zhang Y and Chase SM 2013; Suminski et al. 2013; Kowalski et al. 2013; Orsborn et al. 2014). Another kind of state-dependent control is used to assist the subject during training with the BMI system, where the assistance is adjusted according to the subject’s performance (Velliste et al. 2008). Other state-dependent control algorithms include (Shanechi et al. 2013), where in order to simultaneously estimate movement trajectory or target intent, the decoding parameters are adjusted as targets are approached. Such state-dependent control was the goal behind the recently-developed speed dampening Kalman filter (Golub et al. 2014), which effectively manipulated the viscous damping forces as a function of trajectory curvature to allow for more stable transitions from moving to stopping. The physical control system viewpoint we espouse here may be a way of integrating these approaches into a simple, unified framework for robust prosthetic applications.

Finally, we should note that we have focused exclusively in this review on kinematic BMI decoders of the type that are commonly used to drive cursors on computer screens. Another class of decoders attempts to extract kinetic information (forces and torques) from recorded neural activity for direct control of force output (Fagg et al. 2009; Nazarpour et al. 2012; Nishimura et al. 2013; Chhatbar and Francis 2013; Oby et al. 2013). It is unclear at present how to best integrate these two approaches to BMI design. One intriguing idea is that the brain represents impedances (relationships between kinetics and kinematics) rather than desired forces or kinematics per se (Hogan 1985; Tin and Poon 2005; Hogan and Sternad 2012). An interesting, active direction of research is to design decoders that seemlessly transition between free movement and object interaction. It is possible that decoding impedance control signals directly from the brain would enable this transition.

4.2 Do details of the control signal mapping matter?

We have focused primarily on the overall form of the physical control system represented by different decoders: e.g., whether it is 1st order or 2nd order and whether it contains both viscous and elastic elements, etc. We have not focused on the particular values of those parameters as much, except to point out that the PVA and OLE, which have different mappings between firing rates and cursor movement, perform similarly on-line (Chase et al. 2009). Do the details of the control system parameters matter at all?

The answer is certainly yes. Sadtler and colleagues have recently demonstrated that some mappings between neural activity and cursor movement are easily learned, while others are very difficult to learn (Sadtler et al. 2014). In that study they discovered that subjects could easily learn to map an existing pattern of neural activity to an arbitrary movement; what was difficult was creating new correlation patterns in the neurons themselves. This suggests that, provided the mapping from neural activity to cursor movement preserves the natural correlations within the neural population, two physical systems of the same form may be equivalently usable.

Also, while Sadtler and colleagues noted substantial improvement in performance with short amounts of practice on novel ‘within-manifold’ (correlation preserving) mappings, the performance under those mappings was still substantially worse than performance with the initially-estimated natural tuning curves. Further, Gilja and colleagues (Gilja et al. 2012) have noted a significant improvement in control if tuning curves are re-estimated from on-line training data, which demonstrates different mappings between firing rates and cursor movement can lead to different on-line performance.

Under the estimation-framework, performance will be best with mappings from neural activity to cursor movement that best capture the natural neural tuning. Under the control framework, performance will be best with those mappings that most accurately capture the volitionally usable correlation structures within the neural population. Of course, all of this ignores long-term learning, which may allow subjects to become proficient even at non-intuitive mappings (Ganguly and Carmena 2009; 2010). Further work will need to be done to determine the precise details that determine which mappings from neural activity to movement perform better than others, and how learning may ultimately play a role in sculpting their performance (Orsborn et al. 2014).

4.3 Alternate views of motor cortical recruitment

Decades of motor control studies have established that neural activity in motor cortex correlates with various features of movement. In BMI design, it is common to interpret these correlations by thinking of neurons as being tuned to the desired outcome or intended movement of an effector. However, neural activity can be flexibly dissociated from the effector (Schieber 2011): using operant conditioning, individual neurons can be trained to correlate and decorrelate from particular muscles (Fetz 1969; Fetz and Finocchio 1971), even when spike triggered averaging of EMG traces provide evidence of monosynaptic connections between the neuron and the muscle (Davidson et al. 2007). Evidence of a flexible relationship between neural activity and behavior also comes from natural movement paradigms. Neural activity in motor cortex readily changes during motor learning (Jarosiewicz et al. 2008; Paz et al. 2003; Wise et al. 1998; Gandolfo et al. 2000; Ganguly and Carmena 2009; Sadtler et al. 2014). Neural activity has also been shown to change during associative learning: neurons in M1 will respond to the color of a target during an associative learning task, and will often maintain that tuning after color is no longer relevant (Zach et al. 2008). Even context can modulate neural tuning. Hepp-Reymond and colleagues (Hepp-Reymond et al. 1999) demonstrated that neurons in primary motor cortex are sensitive to the context of an isometric force task, and will change their firing to particular force levels according to how many force targets are presented in the task.

These studies suggest that motor cortical activity is extremely fungible, restructuring itself in the face of new task demands in a manner that is not at all well understood. It is hard to reconcile this view with a static view of tuning to intended movements, unless one assumes that these motor intent signals can themselves be dissociated from the motor outcome (Chase and Schwartz 2010). If this is the case, there may be only subtle differences among the viewpoints that neurons tune to flexible motor intent signals, that neurons act as control signals to drive an effective behavior, or that populations of neurons act as a flexible pattern generator on which movements can be built (Shenoy et al. 2013). In this review, we mainly wish to highlight the importance of these multiple viewpoints when attempting to interpret BMI performance and, ultimately, design the optimal decoding algorithm.

4.4 Relationships to embodiment, internal models, and natural motor control

Natural motor control is fraught with computational difficulties, not least of which is the necessity to compensate for noisy, delayed sensory feedback. To generate fast, dexterous movements, it is necessary to compensate for these sensory delays. It is widely believed that we do this with the aid of internal models that allow us to predict, in real time, the outcomes of our motor commands before sensory feedback becomes available (Crapse and Sommer 2008; Shadmehr et al. 2010). These internal (forward) models are thought to take as input efference copies of our motor commands, and use them to predict the sensory consequences, such as the new arm or eye position, that result from those commands (Sommer and Wurtz 2002; Wolpert et al. 1995). Essentially, these models embody our internal conception of the physics of our limbs and how they respond to our motor commands. These predicted locations can then be used as the basis for planning the next movement before real sensory feedback becomes available, allowing for faster motor sequence production.

Internal models may also be used in BMI control. Motor commands of subjects using a BMI to control a computer cursor in a center-out movement task are more appropriate to the real time position of the cursor than the last sensory feedback position, indicating that subjects compensate for sensory feedback delays while using a BMI (Golub et al. 2012). Further, differences between a subject’s internal model of the decoder and the actual cursor dynamics may explain errors in on-line control (Golub et al. 2013).

It is intriguing to speculate on the utility of internal models in BMI control. Internal models are thought to be a key component in motor adaptation (Shadmehr et al. 2010; Kawato 1999), and thus subjects who build internal models of BMIs may be able to take advantage of all of the computational motor adaptation machinery that natural motor control relies upon. A reliable internal model of the BMI device may also be the key to embodiment, when the device feels like a natural extension of one’s body (Schwartz et al. 2006; Giummarra et al. 2008). However, can subjects build internal models of any BMI device? While the answer to this question is well beyond the scope of this review, we posit that it might be easier to build internal models of physical systems than it is to build internal models of non-physical systems, simply because all of the internal models we build in the context of natural reaching are, by definition, models of physical systems. It could be that the reason decoding algorithms corresponding to physical systems appear to be more useable than those that do not is because they are more readily conceptualized with an internal model.