Self-recalibrating classifiers for intracortical brain–computer interfaces

Animal protocols were approved by the Stanford University Institutional Animal Care and Use Committee. We trained two adult male monkeys (Macaca mulatta, monkeys I and L) to perform center-out reaches for liquid rewards. As illustrated in figure 1, visual targets were back-projected onto a vertical screen in front of the monkey. A trial began when the monkey touched a central yellow square. Following a randomized (300–500 ms) hold period, a visual reach target appeared. After a brief (35 ms) delay, the animal was permitted to reach. Reach targets could appear at one of 28 locations. Each of the 28 targets was typically located at a distance of 4.5, 7.0, 9.5, or 12 cm from the central target and at an angle of 40°, 85°, 130°, 175°, 220°, 310°, or 355° from horizontal so that the targets were arranged in 4 concentric rings of 7 targets each. For all analyses in this paper, we collapse across reach distance and only consider the direction of reach on each trial. After the reach target was held for 200 ms for monkey L and 300 ms for monkey I, a liquid reward was delivered.

Zoom In Zoom Out Reset image size

Monkeys sat in a custom chair (Crist Instruments, Hagerstown, MD) with the head braced and nonreaching arm comfortably restrained. A 96-channel silicon electrode array (Blackrock Microsystems, Salt Lake City, UT) was implanted in motor cortex contralateral to the reaching arm. When identifying neural spikes on each electrode, threshold crossings (Chestek et al 2009, Fraser et al 2009) were used in place of spike sorting. For monkey L, we also required voltage waveforms to pass a shape heuristic. Threshold levels for each electrode were set at the beginning of each recording day with an automated method and then, with possible rare exception, held fixed for the remainder of data collection. The automated method of setting threshold levels consisted of collecting 2 s of filtered voltage data sampled at 30 kHz from each electrode and calculating a root-mean-square like value from this data. Thresholds levels were then set at three times this value for each electrode.

Recordings were made on a total of 41 (monkey L) and 36 (monkey I) days spanning a period of 48 (monkey L) and 58 (monkey I) days in total. In all of the work that follows, we limit ourselves to only analyzing successful trials. In the data analyzed, there were (mean ± s.d.) 1737 ± 482 (monkey L) and 1688 ± 368 (monkey I) successful trials per day spanning 53 ± 16 (monkey L) and 88 ± 16 (monkey I) min per day.

We aim to develop BCI classifiers which maintain stable day-to-day performance without any daily supervised retraining. For clarity, we will speak of developing classifiers within the context of the datasets used in this study, but the techniques presented here should generalize to other intracortical BCI classification scenarios.

Throughout this work, we rely on probabilistic methods to infer the intended reach direction from neural data. For trial t on day d, we form a feature vector, $\boldsymbol{x}_{d,t} \in \mathbb {N}^{E \times 1}$, from our neural data by counting the number of threshold crossings on each of E electrodes of the recording array. The notation x_{d, e, t} will refer to the specific spike count for electrode e¹⁴. Similarly, y_{d, t} ∈ {1, ..., J} indicates which one of J directions a reach is made in for a specific trial and day. Throughout this work we follow the convention of using bold script exclusively to denote vectors. The symbol P will be used to indicate probability distributions.

On any given trial there is a probabilistic model relating observed spike counts to reach direction. We model this with the distribution $P(\boldsymbol{x}_{d,t}| y_{d,t}; \Theta _{d,t})$, where Θ_{d, t} is a set of parameters characterizing the distribution. In this work for electrode e, these parameters consist of the means μ_{e, j} and variances v_{e, j} of spike counts when reaches are made in each direction j. We will refer to this set of parameters as tuning parameters. While both means and variances can drift in principle, we focus on achieving stable classification through tracking changes in mean spike counts. We term the set of μ_{e, j} parameters for a fixed reach direction j as class means, as they are mean counts conditioned on a class of movement. We can index class means across days and trials as μ_{d, e, j, t}, where d indexes day and t indexes a trial within a day. In some of the work that follows, we will speak of class means across all trials of a day, which will we denote as μ_{d, e, j}.

Movements of the recording array relative to the brain (Santhanam et al 2007), scar tissue buildup (Polikov et al 2005), behavioral changes (Chestek et al 2007) and neural plasticity (Taylor et al 2002, Jarosiewicz et al 2008, Ganguly and Carmena 2009, Chase et al 2012) may all cause tuning parameters to drift. In this section, we track drift in class means across days and time by convolving spike counts for all reach directions in a particular direction with a Gaussian kernel. This is a simple way of smoothing spike counts to estimate tuning parameters when true reach direction is known on each trial. This is an appropriate assumption for this analysis which aims to establish basic properties of tuning parameter drift. However, the techniques introduced in this section should not be confused with the methods that are introduced later to track drifting tuning parameters by our self-recalibrating classifiers when knowledge of true reach direction is not assumed known. Mathematically, estimates of tuning parameters are formed through convolution with a kernel as

$$\begin{equation} \hat{\mu }_{d,e,j,t} = \frac{1}{C} \sum_{s: y_{d,s} = j} K(t-s; w)x_{d,e,s}, \end{equation} \tag{ 1 }$$

where K(t − s; w) is Gaussian density function with 0 mean and standard deviation w and $C = \sum _{s: y_{d,s} = j} K(t-s; w)$ is a normalizing constant. Mathematically, the function K indexes over trial number, and therefore, t − s in the argument is the distance in number of trials between two trials and not time. There remains the challenge of selecting the proper value of the width, w, of the Gaussian density function. For the purposes of this paper, we will present results for a single value of w and then show in Supplementary Results (available at stacks.iop.org/JNE/11/026001/mmedia) that the qualitative results we obtain are robust as the value of w is varied.

We are interested in comparing the amount of drift between and within days. If we smooth only within days, we enforce smoothing within but not across days and we may bias our results to show larger between than within-day drift. To account for this, we ignore day boundaries and smooth across all trials for a given electrode and reach direction as if they were on the same day. This may tend to 'over-smooth' between days, but to foreshadow what is to come, we would like to argue that within-day variation in tuning parameters is small relative to between-day variation. Therefore, this procedure is conservative for our purposes.

We desire to examine the impact of tuning parameter drift on a standard BCI classifier. We select a Gaussian, naïve Bayes decoder as our standard classifier, which has been used for BCI classification before (Santhanam et al 2006, Maynard et al 1999). Briefly, this classifier models spike counts on electrodes as independent Gaussian distributions conditioned on reach direction. The notation μ_{e, j} and v_{e, j} refer to the mean and variance of the distribution for electrode e and reach direction j. For notational convenience, we will group all the means and variances for a single electrode e into the vectors $\boldsymbol{\mu }_e \in \mathbb {R}^{J \times 1}$ and $\boldsymbol{v}_e \in \mathbb {R}^{J \times 1}$, respectively. The joint probability of observing spike counts for all electrodes can be computed as the product of probabilities for individual electrodes. Formally, we write this as

$$\begin{equation} P(\boldsymbol{x}_{d,t}| y_{d,t} ) = \prod _{e=1}^E P(x_{d,e,t}| y_{d,t}). \end{equation} \tag{ 2 }$$

A uniform prior probability of reaching in any direction, $P(y_{d,t}) = 1/J$, is employed and Bayes' rule can then be used to find the posterior probability of reaching in direction j given neural data, $P(y_{d,t} = j| \boldsymbol{x}_{d,t})$. This posterior is calculated for each reach direction, and the direction with the highest probability is selected as the decoded direction for each trial.

In this work we will examine the performance of two versions of this standard classifier, which differ only in the way they are trained. The first version, referred to as the standard retrained classifier, learns the $\boldsymbol{\mu }_e$ and $\boldsymbol{v}_e$ parameters fresh with trials from the beginning of each day on which classification will be performed. This emulates standard practice in the BCI field. The second version, referred to as the standard non-retrained classifier, is trained initially on all trials collected on a small number of days in which neural data and desired reach direction is known. After this training, the parameters of the standard non-retrained classifier are held fixed, and it is then used to predict reach direction for trials on subsequent days without any type of retraining or recalibration.

We now describe the intuition and probabilistic model behind our two self-retraining classifiers. The first version, which we refer to as the SR classifier, takes a principled probabilistic approach which provides insight into the simplified SR (SR_S) classifier we subsequently present. Readers interested only in the SR_S classifier may skip to section 2.6 after reading the part of this section outlining the probabilistic model common to both classifiers without loss of understanding.

The SR classifier can be considered an augmented standard classifier. It continues to assume tuning parameters are fixed within a day, but it no longer assumes class means are known at the beginning of each decoding day. Mathematically, it achieves this by introducing a distribution over class mean values that a decoding day begins with. As a decoding day progresses, the SR classifier produces refined distributions over class means while simultaneously estimating reach direction on a trial-by-trial basis.

A formal description of the probabilistic model for the SR classifier follows. When explicitly referring to the class mean for electrode e, reach direction j on day d, we will use the notation μ_{d, e, j}. In results we show that the μ_{d, e, j} parameters for the same electrode drift from day to day in a highly correlated manner. Intuitively, this means that the average spike counts observed for two different reach directions on the same electrode tend to move up and down together across days. We capture this intuition by referencing all class means for electrode e to the same base value b_{d, e} for day d, and we allow this base value to be drawn from a Gaussian distribution independently on each day. Mathematically,

$$\begin{equation} b_{d,e} \sim \mathcal {N} (m_e, s_e), \end{equation} \tag{ 3 }$$

where m_e is the mean and s_e is the variance of the distribution from which base values are drawn. For notational convenience we will group all the base values for all electrodes for day d into a vector $\boldsymbol{b}_d \in \mathbb {R}^{E \times 1}$. We model the base values for different electrodes as independent, and the prior probability for the set of base values for all electrodes is therefore given by

$$\begin{equation} P(\boldsymbol{b}_d) = \prod _{e=1}^E P(b_{d,e}). \end{equation} \tag{ 4 }$$

We relate the class mean for reach direction j for an electrode to its base value through an offset, o_{e, j}, so that

$$\begin{equation} \mu _{d,e,j} = o_{e,j} + b_{d,e}. \end{equation} \tag{ 5 }$$

The offset, o_{e, j}, is fixed across days. Class means for a single electrode can then vary from day to day but they must move up and down across days together. Having defined a probabilistic model for b_{d, e}, we specify an observation model for a single electrode by defining the probability of a spike count on electrode e for trial t on day d given b_{d, e} and reach direction j as

$$\begin{equation} x_{d,e,t} | b_{d,e}, y_{d,t} = j \sim \mathcal {N} (o_{e,j} + b_{d,e}, v_{e,j}). \end{equation} \tag{ 6 }$$

This observation model is almost identical to that of the standard classifier with one important exception: spike counts now depend on the base value, which is a random variable that can vary from day to day. To specify the complete observation model, we again model electrodes as conditionally independent given reaching direction so that

$$\begin{equation} P(\boldsymbol{x}_{d,t}| \boldsymbol{b}_d, y_{d,t}) = \prod _{e=1}^E P(x_{d,e,t} | b_{d,e}, y_{d,t}). \end{equation} \tag{ 7 }$$

To complete the specification of the probabilistic model for the SR classifier, we define the prior probability of reaching in any direction as uniform, exactly as specified for the standard model.

2.5.1. Decoding with the self-recalibrating classifier

We now show how to decode reach direction while simultaneously refining knowledge of class means with the SR classifier. For each trial, we seek to find two distributions. The first is $P(y_{d,t} | \lbrace \boldsymbol{x}_{d,s} \rbrace _{s=1}^t)$, the posterior distribution over reach directions. In this work the notation $\lbrace \boldsymbol{x}_{d,s} \rbrace _{s=1}^t$ will be used to indicate the set of trials from 1 through t on day d. The most probable direction from $P(y_{d,t} | \lbrace \boldsymbol{x}_{d,s} \rbrace _{s=1}^t)$ is selected as the decoded reach direction on each trial. The second distribution is $P(\boldsymbol{b}_d | \lbrace \boldsymbol{x}_{d,s} \rbrace _{s=1}^t)$, which embodies the decoder's knowledge of the base values for a day given the neural activity that it has seen.

As a decoding day progresses, the decoder's knowledge of the true base values for a day, as represented in the distribution $P(\boldsymbol{b}_d | \lbrace \boldsymbol{x}_{d,s} \rbrace _{s=1}^t)$, becomes more refined. This is possible even though the decoder is never provided with knowledge of the user's true reach directions. The key intuition is that spike counts will tend to cluster around class means for a day, and this structure can be exploited by our classifier, as illustrated for a one electrode and two reach direction decoding scenario in figure 2.

Zoom In Zoom Out Reset image size

We now describe the formal decoding procedure for arriving at the distributions $P(y_{d,t} | \lbrace \boldsymbol{x}_{d,s} \rbrace _{s=1}^t)$ and $P(\boldsymbol{b}_d | \lbrace \boldsymbol{x}_{d,s} \rbrace _{s=1}^t)$ for each trial. The algorithm is recursive and will use the results from the previous trial, the distributions $P(y_{d,t-1} | \lbrace \boldsymbol{x}_{d,s} \rbrace _{s=1}^{t-1})$ and $P(\boldsymbol{b}_d | \lbrace \boldsymbol{x}_{d,s} \rbrace _{s=1}^{t-1})$, as its starting point in its calculations for trial t. When decoding the first trial, the prior distribution on reach directions, P(y_{d, t}), is used in place of $P(y_{d,t} | \lbrace \boldsymbol{x}_{d,s} \rbrace _{s=1}^{t-1})$ and the prior distribution on base values, $P(\boldsymbol{b}_d)$, is used in place of $P(\boldsymbol{b}_d | \lbrace \boldsymbol{x}_{d,s} \rbrace _{s=1}^{t-1})$.

Calculating $P(y_{d,t}| \lbrace \boldsymbol{x}_{d,s}\rbrace _{s=1}^t)$.

For trial t, we calculate $P(y_{d,t}| \lbrace \boldsymbol{x}_{d,s}\rbrace _{s=1}^t)$ using Bayes' rule as

$$\begin{eqnarray} \fl P\big(y_{d,t} &=& j| \lbrace \boldsymbol{x}_{d,s}\rbrace _{s=1}^t\big) \nonumber\\ &=& \frac{ P\big(\boldsymbol{x}_{d,t} | \lbrace \boldsymbol{x}_{d,s} \rbrace _{s=1}^{t-1}, y_{d,t} = j\big) P\big(y_{d,t} = j| \lbrace \boldsymbol{x}_{d,s} \rbrace _{s=1}^{t-1}\big) }{P\big(\boldsymbol{x}_{d,t}| \lbrace \boldsymbol{x}_{d,s} \rbrace _{s=1}^{t-1} \big)}\nonumber\\ &=& \frac{c_{t|j}P(y_{d,t})}{c_t}, \end{eqnarray} \tag{ 8 }$$

where in moving to the second line we have defined $c_{t|j} = P(\boldsymbol{x}_{d,t} | \lbrace \boldsymbol{x}_{d,s} \rbrace _{s=1}^{t-1}, y_{d,t} = j)$, made use of the fact that reach direction is independent of neural activity so that $P(y_{d,t} = j| \lbrace \boldsymbol{x}_{d,s} \rbrace _{s=1}^{t-1}) = P(y_{d,t})$, and recognized that the denominator of the first line is a normalizing constant calculated as $c_t = \sum _{j=1}^J c_{t|j} P(y_{d,t})$.

The scalar c_t|j is the probability of the neural activity for trial t if the subject reaches in direction j given all past neural activity seen during the day. We will soon describe a means of calculating its value, but first note that it appears in the calculations for another distribution, $P(\boldsymbol{b}_d | \lbrace \boldsymbol{x}_{d,s}\rbrace _{s=1}^{t}, y_{d,t} = j)$, which we will need shortly. This distribution gives the probability of base values conditioned on all observed neural activity if the subject were to reach in direction j for the current trial. To calculate it, we again use Bayes' rule to write

$$\begin{eqnarray} &&\fl P\big(\boldsymbol{b}_d | \lbrace \boldsymbol{x}_{d,s}\rbrace _{s=1}^{t}, y_{d,t} = j\big) \nonumber \\ &&= \frac{P(\boldsymbol{x}_{d,t}| \boldsymbol{b}_d, y_{d,t} = j)P\big(\boldsymbol{b}_d | \lbrace \boldsymbol{x}_{d,s}\rbrace _{s=1}^{t-1}\big)}{c_{t|j}}, \end{eqnarray} \tag{ 9 }$$

where $P(\boldsymbol{x}_{d,t}| \boldsymbol{b}_d, y_{d,t} = j)$ is given by equation (7) and $P( \boldsymbol{b}_d | \lbrace \boldsymbol{x}_{s}\rbrace _{s=1}^{t-1})$ is the posterior over base values obtained when decoding the previous trial. In forming the numerator of equation (9), we have used the fact that $\boldsymbol{x}_{d,t}$ conditioned on $\boldsymbol{b}_d$ is independent of all previous neural observations and that base value is independent of y_{d, t} conditioned on the previous neural activity. We see that the numerator on the right-hand side of equation (9) is in the form of the product of two Gaussians. $P(\boldsymbol{b}_d | \lbrace \boldsymbol{x}_{d,s}\rbrace _{s=1}^{t-1})$ is defined to be a Gaussian distribution with mean $\boldsymbol{m}_{t-1}$ and covariance matrix Σ_{t − 1}. Viewed as a function of $\boldsymbol{b}_d$, $P(\boldsymbol{x}_{d,t}| \boldsymbol{b}_d, y_{d,t} = j)$ is Gaussian with mean $\boldsymbol{x}_{d,t} - \boldsymbol{o}_j$, where $\boldsymbol{o}_j \in \mathbb {R}^{E \times 1}$ is a vector of offsets for each electrode for reach direction j, and covariance matrix $V_j \in \mathbb {R}^{E \times E}$, where V_j is a diagonal matrix with the values of v_{e, j} along its diagonal. We can then calculate the mean, $\boldsymbol{m}_{t|j}$, and covariance matrix, Σ_t|j, of $P(\boldsymbol{b}_d | \lbrace \boldsymbol{x}_{d,s}\rbrace _{s=1}^{t}, y_{d,t} = j)$ of the left-hand side of equation (7) as

$$\begin{equation} \Sigma _{t|j} = \big( V_j^{-1} + \Sigma _{t-1}^{-1}\big)^{-1} \end{equation} \tag{ 10 }$$

$$\begin{equation} \boldsymbol{m}_{t|j} = \Sigma _{t|j}\big(V_j^{-1}(\boldsymbol{x}_{d,t} - \boldsymbol{o}_j) + \Sigma _{t-1}^{-1}\boldsymbol{m}_{t-1} \big), \end{equation} \tag{ 11 }$$

where the notation · ⁻¹ indicates matrix inversion. The scalar c_t|j is then the integral over the numerator of equation (9) and can be calculated as

$$\begin{eqnarray} \fl c_{t|j} &=& \int P(\boldsymbol{x}_{d,t}| \boldsymbol{b}_d, y_{d,t} = j)P\big(\boldsymbol{b}_d | \lbrace \boldsymbol{x}_{d,s}\rbrace _{s=1}^{t-1}\big)\ {\rm d}\boldsymbol{x}_{d,t} \nonumber\\ &=& \sqrt{(2\pi )^E {\rm det}( \Sigma _{t|j}) } P(\boldsymbol{x}_{d,t}| \boldsymbol{b}_d = \boldsymbol{m}_{t|j}, y_{t,d} = j) \nonumber \\ && \times P(\boldsymbol{b}_d = \boldsymbol{m}_{t|j}| \lbrace \boldsymbol{x}_{d,s}\rbrace _{s=1}^{t-1}), \end{eqnarray} \tag{ 12 }$$

where $P(\boldsymbol{x}_{d,t}| \boldsymbol{b}_d = \boldsymbol{m}_{t|j}, y_{t,d} = j)$ is obtained by evaluating equation (7) for the current observation at $\boldsymbol{b}_d = \boldsymbol{m}_{t|j}$ and $P(\boldsymbol{b}_d = \boldsymbol{m}_{t|j}| \lbrace \boldsymbol{x}_{d,s}\rbrace _{s=1}^{t-1})$ is obtained by evaluating the posterior over base values from the last trial at the same point, $\boldsymbol{b}_d = \boldsymbol{m}_{t|j}$.

Calculating $P(\boldsymbol{b}_d | \lbrace \boldsymbol{x}_{d,s}\rbrace _{s=1}^{t})$.

Finally, we compute $P(\boldsymbol{b}_d | \lbrace \boldsymbol{x}_{d,s}\rbrace _{s=1}^{t})$, the posterior distribution over the base values for all electrodes. An analytical expression for this can be written as

$$\begin{eqnarray} &&\fl P\big(\boldsymbol{b}_d | \lbrace \boldsymbol{x}_{d,s}\rbrace _{s=1}^{t}\big) = \sum _{j=1}^J P\big(\boldsymbol{b}_d | \lbrace \boldsymbol{x}_{d,s}\rbrace _{s=1}^{t}, y_{d,t} = j\big) \nonumber \\ &&\times P\big(y_{d,t} = j| \lbrace \boldsymbol{x}_{d,s}\rbrace _{s=1}^t\big), \end{eqnarray} \tag{ 13 }$$

where the two terms on the right-hand side of equation (13) are obtained from equations (8) and (9). The right-hand side of equation (13) is a mixture of multivariate Gaussians. This will be problematic as in general there will be an exponentially increasing number of components in the representation of $P(\boldsymbol{b}_d | \lbrace \boldsymbol{x}_{d,s}\rbrace _{s=1}^{t})$ as t increases. Therefore, we approximate the right-hand side of equation (13) with mean $\boldsymbol{m}_t$ and covariance matrix Σ_t using moment matching as has been done in previous hypothesis merging algorithms (Blom and Bar-Shalom 1988). We use the standard formulas for moment matching with two multivariate Gaussians to write

$$\begin{equation} \boldsymbol{m}_t = \sum _{j=1}^J \boldsymbol{m}_{t|j}P\big(y_{d,t} = j | \lbrace \boldsymbol{x}_{d,s}\rbrace _{s=1}^t\big) \end{equation} \tag{ 14 }$$

$$\begin{eqnarray} &&\fl \Sigma _t = \sum _{j=1}^J [\Sigma _{t|j} + (\boldsymbol{m}_{t|j} - \boldsymbol{m}_t)(\boldsymbol{m}_{t|j} - \boldsymbol{m}_t)^T ] \nonumber \\ &&\times P\big(y_{d,t} = j | \lbrace \boldsymbol{x}_{d,s}\rbrace _{s=1}^t\big), \end{eqnarray} \tag{ 15 }$$

where $P(y_{d,t} = j | \lbrace \boldsymbol{x}_{d,s}\rbrace _{s=1}^t)$ is calculated from equation (8). This completes the decoding procedure for trial t. For each new trial, this set of calculations is repeated in an iterative manner.

2.5.2. Learning the parameters for the self-recalibrating classifier

While the SR classifier is designed to produce stable day-to-day classification without daily retraining, it must be initially trained once. To perform this initial training there must exist a small number of days for which neural activity, $\boldsymbol{x}_{d,t}$, has been recorded along with associated reach directions, y_{d, t}. We will use the notation $\mathcal {D}_{{\rm Train}}$ to refer to this set of days. To avoid any possible confusion, we stress that after the classifier has been trained on this initial set of days, it requires no further supervised retraining and will be equipped to autonomously calibrate itself on new data that it encounters during normal decoding use.

The purpose of the initial training is to learn the mean, m_e, and variance, s_e, characterizing the distribution of base values for each electrode from equation (3) and the offset, o_{e, j}, and variance, v_{e, j}, characterizing the observation model in equation (6) for each electrode-class pair. Fitting the SR classifier is more difficult than fitting the standard classifier because the daily base values for each electrode, b_{d, e}, cannot be directly observed. Instead, these values are randomly drawn each day from the distribution specified by equation (3), and after the base values are drawn on each day, they specify the mean of the Gaussian distributions in equation (6) from which neural spike counts are drawn. Thus, the neural data we directly observe is essentially one step removed from the distribution that is characterized by the parameters we desire to estimate. Figure 3 illustrates the conceptual difference between the standard and SR classifier models.

Zoom In Zoom Out Reset image size

The EM algorithm (Dempster et al 1977) is a general procedure for finding the maximum likelihood model parameters in the presence of unobserved variables, in this case b_{d, e}. The EM algorithm is an iterative procedure which is initialized with an initial guess for the model parameters. In this case, these are the m_e, s_e, o_{e, j} and v_{e, j} parameters. Each iteration of the algorithm produces new estimates for these parameters through a two step procedure. In the first step (referred to as the E-step) of each iteration, the latest parameter estimates and observed data are used to form a posterior distribution over the unobserved variables conditioned on the observed data. For the SR classifier, the unobserved variables are the daily base values for each electrode, b_{d, e}, and the observed data are the neural observations, $\boldsymbol{x}_{d,t}$, and reach directions, y_{d, t}, for all trials on all days in the training set. In the second step of each iteration (referred to as the M-step), parameter values are updated so that the expectation of the log likelihood of the observed and unobserved data, computed using the posterior distribution found in the first step of each iteration, is maximized. This algorithm is known to converge to a local optimum. See Supplementary Methods (available at stacks.iop.org/JNE/11/026001/mmedia) for specific details of the EM algorithm used in this work.

2.5.3. Robust classification through outlier detection

The classification procedure described above can be modified to improve robustness on real world data. Key to the probabilistic model underlying the SR classifier is the specification that firing rates can change between but not within days. While minor changes in firing rate throughout a day are tolerable, large changes, which may be observed in real data, can be problematic.

In particular, the SR classifier represents knowledge of the base values for a day obtained from classifying trials 1 through t with the distribution $P(\boldsymbol{b}_d| \lbrace \boldsymbol{x}_{d,s}\rbrace _{s=1}^t)$. This distribution is refined for each trial that is classified, and over time the uncertainty for the base values for a day, represented in the covariance matrix Σ_t, is reduced. This presents no problem in normal operation. However, if the firing rate of an electrode suddenly changes after the uncertainty in base values for a day is reduced, spike counts recorded on that electrode will strongly and erroneously affect classification.

The classification procedure can be modified by adding an additional step to be performed before any of the steps already described when classifying reach direction for trial t. In this step, electrodes with large firing rate changes are flagged, and the uncertainty in the estimate for the base value for those electrodes is increased to a predetermined value. This prevents those electrodes from dominating when determining reach direction probability and maintains classification accuracy. As classification continues, uncertainty around the base value estimates for the flagged electrodes is reduced starting from the predetermined value as the SR classifier refines the distribution $P(\boldsymbol{b}_d| \lbrace \boldsymbol{x}_{d,s}\rbrace _{s=1}^t)$.

To flag erratically behaving electrodes, a distribution which predicts the probability of the spike count on trial t for electrode e given all previous neural data, $P(x_{d,e,t}| \lbrace \boldsymbol{x}_{d,s} \rbrace _{s=1}^{t-1})$, is formed for each electrode as shown in Supplementary Methods (available at stacks.iop.org/JNE/11/026001/mmedia). This can be used as a gauge for determining if the observed spike count, x_{d, e, t}, for a trial and an electrode is anomalous by setting upper and lower thresholds on the value of x_{d, e, t} so that a normally behaving electrode will only exceed these values approximately 1% of the time. Whenever an electrode does exceed these values, it is flagged as behaving erratically. As with any anomaly detection procedure, thresholds have to be set to balance true and false positive rates, and a false positive rate of 1% should not significantly affect classification performance.

After one or more electrodes are flagged as behaving erratically for trial t, uncertainty around the estimates for their associated base values in the distribution $P(\boldsymbol{b}_d| \lbrace \boldsymbol{x}_{d,s}\rbrace _{s=1}^{t-1})$ are increased. For ease of notation, we will momentarily switch from representing the set of base values in vector notation to set notation, so that each entry in the vector $\boldsymbol{b_d}$ will be represented by a corresponding member in the set $\lbrace b_{d,e}\rbrace _{e=1}^E$. We will refer to the set of flagged electrodes for trial t as E_{t, flagged} and the set of non-flagged electrodes as E_{t, non-flagged}. The uncertainty for the base values of erratically behaving electrodes represented in the distribution for $P(\lbrace b_{d,e} \rbrace _{e=1}^E| \lbrace \boldsymbol{x}_{d,s}\rbrace _{s=1}^{t-1})$ is reset in a two step procedure. In the first step, the distribution over only non-flagged electrodes, $P(\lbrace b_{d,e} \rbrace _{e\in E_{{\rm t,non{\scriptsize-}flagged}}}| \lbrace \boldsymbol{x}_{d,s}\rbrace _{s=1}^{t-1})$, is formed. Formally, this is accomplished by forming the marginal distribution over the electrodes in E_{t, non-flagged} from $P(\lbrace b_{d,e} \rbrace _{e=1}^E| \lbrace \boldsymbol{x}_{d,s}\rbrace _{s=1}^{t-1})$ and is easily accomplished since $P(\lbrace b_{d,e} \rbrace _{e=1}^E| \lbrace \boldsymbol{x}_{d,s}\rbrace _{s=1}^{t-1})$ is a Gaussian distribution. In the second step, a new distribution over the base values for all electrodes, which we refer to as $P^*(\lbrace b_{d,e} \rbrace _{e=1}^E| \lbrace \boldsymbol{x}_{d,s}\rbrace _{s=1}^{t-1})$, is formed as

$$\begin{eqnarray} &&\fl P^*\big(\lbrace b_{d,e} \rbrace _{e=1}^E| \lbrace \boldsymbol{x}_{d,s}\rbrace _{s=1}^{t-1}\big) \nonumber \\ &&= P\big(\lbrace b_{d,e} \rbrace _{e\in E_{{\rm t,non{\scriptsize-}flagged}}}| \lbrace \boldsymbol{x}_{d,s}\rbrace _{s=1}^{t-1}\big) \prod _{e \in E_{{\rm t,flagged}}} P^{\prime }(b_{d,e}), \end{eqnarray} \tag{ 16 }$$

where P'(b_{d, e}) are Gaussian distributions. The mean of the distribution P'(b_{d, e}) is set to, m_{e, t − 1}, which is the element for electrode e of the mean vector for the distribution $P(\lbrace b_{d,e} \rbrace _{e=1}^E| \lbrace \boldsymbol{x}_{d,s}\rbrace _{s=1}^{t-1})$. This ensures that the mean of the final distribution, $P^*(\lbrace b_{d,e} \rbrace _{e=1}^E| \lbrace \boldsymbol{x}_{d,s}\rbrace _{s=1}^{t-1})$, is the same as that of the original distribution $P(\lbrace b_{d,e} \rbrace _{e=1}^E| \lbrace \boldsymbol{x}_{d,s}\rbrace _{s=1}^{t-1})$. The variance of the distribution P'(b_{d, e}) is set to s_e, the variance of the distribution that base values are modeled as being drawn from at the start of each decoding day. Computationally, these two steps can be carried out through a series of simple manipulations of the covariance matrix for $P(\lbrace b_{d,e} \rbrace _{e=1}^E| \lbrace \boldsymbol{x}_{d,s}\rbrace _{s=1}^{t-1})$ and increase the uncertainty around base values for suspect electrodes in $P^*(\lbrace b_{d,e} \rbrace _{e=1}^E| \lbrace \boldsymbol{x}_{d,s}\rbrace _{s=1}^{t-1})$. After the distribution for $P^*(\lbrace b_{d,e} \rbrace _{e=1}^E| \lbrace \boldsymbol{x}_{d,s}\rbrace _{s=1}^{t-1})$ is formed, classification, as described by the algorithm in section 2.5.1, is performed where $P(\lbrace b_{d,e} \rbrace _{e=1}^E| \lbrace \boldsymbol{x}_{d,s}\rbrace _{s=1}^{t-1})$ is replaced with $P^*(\lbrace b_{d,e} \rbrace _{e=1}^E| \lbrace \boldsymbol{x}_{d,s}\rbrace _{s=1}^{t-1})$ as a starting point when decoding trial t.

The SR classifier addresses the problem of base values, $\boldsymbol{b}_d$, which can drift from day to day by starting with a prior distribution, $P(\boldsymbol{b}_{d})$, over these values and then forming refined distributions, $P(\boldsymbol{b}_{d} | \lbrace \boldsymbol{x}_{d,s} \rbrace _{s=1}^{t})$, over base values based on neural activity seen throughout the decoding day. When predicting reach direction for a trial, the SR classifier then considers the probability of reach directions under all possibly values for $\boldsymbol{b}_d$ and weights its estimate according to how likely each value of $\boldsymbol{b}_d$ is, given by $P(\boldsymbol{b}_{d} | \lbrace \boldsymbol{x}_{d,s} \rbrace _{s=1}^{t})$.

Thus, the SR classifier never assumes any single value of $\boldsymbol{b}_d$ is correct, but decodes by considering an infinite number of $\boldsymbol{b}_d$ values, each weighted by some probability. This represents a fully probabilistic means of decoding. However, it is not the only possible means of accounting for drifting tuning parameters. In particular, we can forego modeling uncertainty in $\boldsymbol{b}_d$ and instead decode with a single point estimate for $\boldsymbol{b}_d$ formed for each trial. For the simplified classifier, $\boldsymbol{b}_d$ is defined so that b_{d, e} is the average number of spike counts on electrode e for any trial on day d. We can express this in mathematical notation as $b_{d,e} = \mathbb {E}[ x_{d,e,t}]$, where the expectation is taken across trials.

An estimate for b_{d, e} on trial t, which we refer to as $\hat{b}_{d,e,t}$, can be formed by taking the empirical average of all of the observed spike counts for electrode e on day d up to and including trial t. This estimate can be recursively formed as

$$\begin{equation} n_{d,t} = n_{d,t-1} + 1, \end{equation} \tag{ 17 }$$

$$\begin{equation} \hat{b}_{d,e,t} = \frac{n_{d,t-1}\hat{b}_{d,e,t-1} + x_{d,e,t}}{n_{d,t}} . \end{equation} \tag{ 18 }$$

We must select initial values for n_{d, 0} and b_{d, e, 0}, where the hat has been dropped on b_{d, e, 0} to indicate this is a parameter of the model specifying a value with which to initialize equation (18). If we select n_{d, 0} = 0, $\hat{b}_{d,e,t}$ will be the true empirical average for trials 1 through t. However, this has one disadvantage. Because of the very fact that equation (18) forms estimates as the empirical average of spike counts, at the start of the day there will be no previous trials with which to average spike counts for initial trials with and hence values for $\hat{b}_{d,e,t}$ will be strongly influenced by the spike counts of each trial. This will result in large fluctuations in $\hat{b}_{d,e,t}$ for the early trials of a decoding session. An alternative is to select b_{d, e, 0} so that it represents the expected average spike count on electrode e for any given day and assign this some initial weight by setting n_{d, 0} > 0. Intuitively, we can think of starting each day with an initial estimate of the average spike count for the electrode that was obtained as the empirical average of n_{d, 0} 'virtual samples'. This initial estimate is analogous to the prior distribution for the fully probabilistic classifier specified in equation (3). Importantly, as shown in Supplementary Methods (available at stacks.iop.org/JNE/11/026001/mmedia), as long as n_{d, 0} and b_{d, e, 0} are chosen to be finite, $\hat{b}_{d,e,t}$ approaches the true mean firing spike count for electrode e as the number of trials observed during a day grows. While this is reassuring, there are still values of these parameters that will result in optimal decoding, and below we discuss how values for n_{d, 0} and b_{d, e, 0} can be learned from the training data.

After an estimate for $\hat{b}_{d,e,t}$ is formed for each trial, an estimate for the class means, μ_{d, e, j}, can be obtained with equation (5), and a classification of reach direction can be made using the standard classifier described in section 2.4 with the updated class means. The SR_S classifier can be conceptualized as consisting of two algorithms acting together. The first algorithm tracks base values, and hence tuning parameters, using an accumulated average of spike counts. This simple first algorithm is independent of and provides input in the form of updated tuning parameters to the second algorithm, which serves as the decoding core and is nothing more than the standard classifier. Figure 4 presents a schematic depiction of how the SR_S classifier decodes a single trial and depicts the separation of the self-recalibrating and classification portions of the SR_S decoding algorithm. See Supplementary Methods (available at stacks.iop.org/JNE/11/026001/mmedia) for a more formal description of the relationship between the SR and SR_S classifiers. We note this method of tracking baseline firing rates is similar to the work of Li et al (2011), who in their work on adaptive decoders for continuous control signals experimented with tracking changes in baseline firing rates by measuring the empirical average of spike counts. However, since their algorithm runs in a batch mode and does not update parameters on a single trial basis as our does, they do not introduce the concept of selecting initial values and weights to begin the empirical average for a day with.

Zoom In Zoom Out Reset image size

2.6.1. Learning model parameters for the simplified self-recalibrating classifier

As with the SR classifier, the SR_S classifier must also be trained once on an initial set of training days. The goal of this training is to learn the electrode-class variances, v_{e, j}, and the electrode-class offsets, o_{e, j}, for all electrodes and reach directions. In addition, values for b_{d, e, 0}, must be learned for each electrode, and the initial weight, n₀, which will be assigned to these means must be determined.

Given a collection of days, $\mathcal {D}_{{\rm Train}}$, we estimate values for b_{d, e, 0}, o_{e, j} and v_{e, j} as follows. First, for each day in $\mathcal {D}_{{\rm Train}}$ we estimate the average spike counts on each electrode over all reach directions, $\hat{\mu }_{d,e}$, and the average spike counts on each electrode when the subject made reaches in direction j, $\hat{\mu }_{d,e,j}$. We calculate

$$\begin{equation} \hat{\mu }_{d,e} = \frac{ \sum _{t = 1}^{T_d} x_{d,e,t} }{T_d}, \end{equation} \tag{ 19 }$$

and

$$\begin{equation} \hat{\mu }_{d,e,j} = \frac{\sum _{t: y_{d,t}= j} x_{d,e,t}}{T_{d,j}}, \end{equation} \tag{ 20 }$$

where T_d is the number of trials on day d, and T_{d, j} is the number of trials with reaches in direction j on day d. Having calculated these intermediate values, an estimate of b_{d, e, 0} is then calculated as the average of the daily average spike counts as

$$\begin{equation} \hat{b}_{d,e,0} = \frac{\sum _{d \in \mathcal {D}_{{\rm Train}}}\hat{\mu }_{d,e}}{|\mathcal {D}_{{\rm Train}}|}, \end{equation} \tag{ 21 }$$

where $|\mathcal {D}_{{\rm Train}}|$ is the number of days in $\mathcal {D}_{{\rm Train}}$. Estimates for the class-electrode offsets, o_{e, j}, are calculated as the average of the difference between daily class-electrode mean spike counts and the daily average spike counts. We write this as

$$\begin{equation} \hat{o}_{e,j} = \frac{\sum _{d \in \mathcal {D}_{{\rm Train}}} \hat{\mu }_{d,e,j} }{|\mathcal {D}_{{\rm Train}}|} - \hat{\mu }_{d,e}. \end{equation} \tag{ 22 }$$

Finally, the electrode-class variances, v_{e, j}, are estimated as

$$\begin{equation} \hat{v}_{e,j} = \frac{\sum _{d\in \mathcal {D}_{{\rm Train}}} \sum _{t: y_{d,t}{= j}} (x_{d,e,t} - \hat{\mu }_{d,e,j})^2}{( \sum _{d \in \mathcal {D}_{{\rm Train}}} T_{d,j}) - 1}, \end{equation} \tag{ 23 }$$

which can be understood as a generalization of the standard formula for the variance of a sample generalized to account for values of μ_{d, e, j} which may change between days.

Learning the value of n₀.

The equations above are used to estimate values for the b_{d, e, 0}, o_{e, j} and v_{e, j} parameters of the model. We turn to cross-validation to learn a value for n₀, which is the amount of weight or number of virtual samples that should be assigned to the initial estimate of average electrode spike counts for a day, b_{d, e, 0}.

In the cross-validation procedure, we seek the value of n₀ which will produce the highest classification accuracy on novel data. To do this, multiple rounds of training and validation are performed on the data that is present in $\mathcal {D}_{{\rm Train}}$. In each round, one day in $\mathcal {D}_{{\rm Train}}$ is set aside as a validation day. On the remaining days, the b_{d, e, 0}, o_{e, j} and v_{e, j} parameters of the model are learned. Importantly, these values are learned on all days in $\mathcal {D}_{{\rm Train}}$ except the validation day. This is critical because cross-validation uses the validation day to simulate classification on novel data, and it is therefore important that none of the trials in the validation day are used for learning any part of the model for a round of cross-validation. After learning these parameters, a value of n₀ is selected and the classifier is used to predict reach direction for each trial in the validation day. This is repeated for a range of n₀ values. For a single round of cross-validation, one or more values of n₀ will produce the highest classification accuracy on the validation day. To select the final value of n₀ for use in the trained classifier, a round of cross validation is performed for each day in $\mathcal {D}_{{\rm Train}}$, and the value of n₀ which produces the highest average accuracy across all rounds is chosen.

Figures 5 and 6 show the results of the analysis of drifting tuning parameters for monkeys L and I, respectively. Throughout this work threshold crossings were collected in a single window per trial. In this analysis and in all of the work that follows we count spikes falling in a single 250 ms long window for each trial starting 150 ms after the go cue. As described in section 2.3, spike counts from each electrode were then convolved with a Gaussian kernel to produce smoothed estimates of tuning parameters across trials. The results presented here were obtained with a kernel width of 100 trials. Results obtained with kernel widths of 25, 50 and 200 trials were qualitatively similar (see Supplementary Results, available at stacks.iop.org/JNE/11/026001/mmedia).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Panel A in figure 5 shows the traces of two tuning parameters for one electrode for monkey L. Traces for monkey I were similar but are not shown. Three important properties of tuning parameter drift are apparent in this plot. First, tuning parameters change much more between than within days. Second, tuning parameters for the same electrode tend to move up and down together across days. Finally, tuning parameters do not tend to drift progressively further from their values on day 1 over time.

Figures 5(B) and 6(A) provide a quantitative measurement of the observation that tuning parameters appear to drift more between than within days. The drift of a tuning parameter within a single day is quantified as the standard deviation of its estimated values over all of the trials within a day, σ_{d, e, j}. This can be measured for all days, electrodes and reach directions. The distribution of all σ_{d, e, j} values is plotted as the dashed line in each figure. Similarly, drift between days is quantified as the standard deviation of mean tuning parameter values for each day, σ_{e, j}. The distribution of σ_{e, j} parameters is plotted as the solid line in both figures. The mean of the between-day standard deviation values is significantly larger than that of the within-day standard deviation values (permutation test, p < 0.01). This finding motivates the decision to simplify the self-recalibrating classifiers so that they ignore within-day drift but compensate for between-day drift in tuning parameters.

Figures 5(D) and 6(B) provide evidence that tuning parameters on the same electrode drift in a correlated manner. Correlation coefficients for all pairs of daily mean tuning parameters (the black dots for the example pair of tuning parameters in figure 5(C)) on the same electrode were calculated for each electrode. We also examined correlation across electrodes by measuring the correlation of daily mean values for pairs of tuning parameters for the same reach direction across all pairs of electrodes. The distributions of these two sets of correlation coefficients are shown in figures 5(D) and 6(B). For both subjects, correlation coefficients for tuning parameters on the same electrode concentrate near 1. This insight has important practical implications for the design of self-recalibrating classifiers which will be detailed in section 4.

While tuning parameters drift in a correlated manner on the same electrode, it is important to point out that the distributions of correlation coefficients in figures 5(D) and 6(B) leave room for tuning parameters to move in a limited manner relative to one another while still moving up and down together. Indeed, this is visible in figure 5(A) where traces of the two tuning parameters move up and down together between days but are closer to one another on days 1–4 than they are for days 7–10, for example. We will show below that our decoders, which use a simplified design that models parameters on the same electrode as moving in perfect lockstep, achieve good performance, but future decoders could be developed to model additional aspects of tuning parameter drift. Interested readers may refer to Supplementary Results (available at stacks.iop.org/JNE/11/026001/mmedia), where we quantify changes in modulation depth across time.

Figures 5(E) and 6(C) measure the tendency of tuning parameters to progressively drift from their value on day 1. The percent change of the mean tuning parameter values from day 1, Δ_{d, e, j}, for days 2 and on for all tuning parameters was calculated. The average of the absolute value of the percent change in tuning parameter values across all electrodes and reach directions for each day was then plotted in these figures. When performing this analysis, for numerical reasons we removed any electrodes which had a value of less than .001 spikes/window for any reach direction. The results are shown in figures 5(E) and 6(C). If tuning parameters drift in an unconstrained manner, the average of the unsigned Δ_{d, e, j} values should monotonically increase with time. However, average values do not appear to systematically increase with time, supporting the conclusion that tuning parameters can drift from day to day but do so in a constrained manner.

We compare the performance of the SR, SR_S, standard retrained and standard non-retrained classifiers on the two prerecorded datasets described in section 2.1. The large number of days in each makes these ideal datasets for evaluating classifier stability.

The standard retrained classifier is trained fresh each day using the first 400 trials of the day. After training, the classifier is used to predict reach direction on each of the remaining trials of the day. Importantly, while true reach direction is known for each trial of the datasets and is used for training, this information is not made available to the classifier when predicting reach direction (testing) and is only used after a prediction is made to evaluate classifier accuracy.

The SR, SR_S and standard non-retrained classifiers must be initially trained on a small number of days where neural activity is recorded and the desired reach direction for trials is known. After this they require no further training when classifying reach direction on days when only neural activity is available. For this analysis, all trials on the first ten days of each dataset are used for training these classifiers. After this training, the classifiers are used to classify reach direction for trials on the remaining days of each dataset. To allow a fair head-to-head comparison with the standard retrained classifier, classification on each test day for these classifiers begins at trial 401 of each day. In this way, the set of trials that each classifier predicts reach direction for is the same on any given day. The breakup of data for training and testing each of the classifiers is shown in figure 7(B).

Zoom In Zoom Out Reset image size

When training all classifiers, any electrode with an average of less than two spike counts per trial on the training data (the first 400 trials of each day for the standard retrained classifier and all trials on days 1–10 for the standard non-retrained, SR and SR_S classifiers) is not used when subsequently decoding reach direction on test data. This step prevents these electrodes from causing problems if their average spike counts were to rise at a later point in time when performing classification. Critically, this can be done in clinical practice as the criterion for exclusion is based only on training data.

We now examine the implications of tuning parameter drift on the standard non-retrained classifier. The performance of this classifier indicates what might happen if we simply used the standard classifier without retraining and did nothing to account for drifting tuning parameters. Of course, this is an offline analysis and subjects may be able to adapt their neural signals through cognitive strategies in a closed-loop setting to improve BCI performance (Chase et al 2009), but the results presented here may still be considered as a type of lower bound on performance without retraining.

Examination of figure 8, which displays the average daily classification accuracy of the standard non-retrained classifier for both subjects, reveals two important findings. First, performance of the standard non-retrained classifier is substantially lower than that of the daily retrained classifier on most days. This motivates the development of classifiers which account for drifting tuning parameters. However, classification accuracy of the standard non-retrained classifier does not significantly decline with time, which constitutes the second finding. To quantify this claim, the slope of a linear fit to daily classification accuracy is not significantly different than zero for both subjects. Specifically, the 95% confidence interval for the slope, measured in units of percent accuracy per day, is ( − 0.006, 0.006) for monkey L and ( − 0.012, 0.001) for monkey I, which contains zero for both monkeys. While this finding might at first seem surprising, it is consistent with the finding above that tuning parameters do not show a systematic tendency to progressively wander from their values on day 1 but seem to vary within a relatively limited range.

Zoom In Zoom Out Reset image size

Figures 9 and 10 show the daily classification accuracies of the SR and SR_S classifiers for each subject. We emphasize that after these classifiers were trained on days 1–10 for a dataset, they were not provided with any information about the true reach direction of trials with which to retrain themselves when predicting reach direction on the following days. Instead, using the techniques described in section 2, the SR and SR_S classifiers simultaneously predicted reach direction and the daily value of tuning parameters from neural data alone. Therefore, the classification accuracies reported here are indicative of what might be achieved in the clinic without any further daily supervised retraining or recalibration after the system is initially trained.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 1 compares the overall accuracies, calculated as the mean of the average daily accuracies, for all four classifiers considered in this work. The overall accuracies of the SR and SR_S classifiers were both 14% higher than that of the standard non-retrained classifier for monkey L and 13% and 16%, respectively, higher for monkey I. These performance improvements were statistically significant for both self-recalibrating classifiers and both subjects (One Sided Wilcoxon Signed-Rank Test, p < 10⁻⁶ for both classifiers for monkey L, p < 10⁻⁵ for both classifiers for monkey I). Further, both self-recalibrating classifiers nearly recovered the accuracy achieved by the standard retrained classifier. The overall accuracy of the SR and the SR_S classifiers were both 1% higher than that of the retrained classifier for monkey L and only 5% and 3% lower for monkey I. These differences were not statistically significant for either self-recalibrating classifier and monkey L (One Sided Wilcoxon Signed-Rank Test, p > 0.1 for SR classifier and p > 0.05 for the SR_S classifier), while the small decrease in accuracies observed on the data for monkey I were (One Sided Wilcoxon Signed-Rank Test, p < 10⁻⁴ for SR classifier and p < 10⁻² for the SR_S classifier). To be clear, the self-recalibrating classifiers decoded the exact same trials, trials 401 and on for each test day, as the retrained classifier and achieved this performance without the benefit that the retrained classifier had of being retrained on the fully labeled data of the first 400 trials of each day.

An important consideration when evaluating the performance of a self-recalibrating classifier is the number of trials required to learn tuning parameter values. A self-recalibrating classifier is of little use if its performance starts near chance each day and then takes hundreds of trials to reach acceptable prediction accuracy. Figure 11 shows the average classification accuracy achieved by each of the self-recalibrating classifiers throughout the course of a decoding day for both subjects. These accuracies were calculated by dividing all of the trials decoded during a test day into sequential bins of 20 trials each and then calculating the average prediction accuracy in each of these bins over all days of data for a subject. The number of trials in each day varied, and the number of bins that accuracies could be calculated over was limited by the day with the fewest number of trials for each subject. Monkey L had at least 386 trials for classification on each test day and monkey I had at least 752. Examination of the plots in figure 11 reveal that both classifiers start a decoding day well above chance levels of performance. Further, the SR classifier shows little evidence of a 'burn-in' period, where accuracy might be lower as tuning parameter values are being learned. It is likely that if there were more days of data present and smaller bin sizes could be used in the analysis, a burn-in period of less than 20 trials might become apparent, but in any case the SR classifier rapidly achieves steady state performance. The SR_S classifier shows evidence of slightly decreased performance in the first bin of 20 trials for both subjects, and performance appears to be at asymptotic levels by trials 21–40.

Zoom In Zoom Out Reset image size

Figure 12 displays tuning parameter estimates on two electrodes produced by both self-recalibrating classifiers on the first five days of test data for monkey L. Results for monkey I were similar. For this analysis, estimates for true tuning parameters were obtained by convolved spikes with a Gaussian kernel with a standard deviation of 100 trials as described in section 2.3 with one exception: days were processed separately and not concatenated together. The SR classifier produces a distribution over tuning parameter values on each trial (shown as the distributions over μ_{d, e, j = 1} and μ_{d, e, j = 2} in figure 2), and the mean of these distributions, $\mathbb {E}[\mu _{d,e,j}| \lbrace \boldsymbol{x}_{d,s} \rbrace _{s=1}^t]$, was used to plot the traces in figure 12. Even though there are seven tuning parameters per electrode, traces for only two tuning parameters per electrode are shown for visual clarity. Characteristics of the way both classifiers predict tuning parameters become apparent with inspection of the plots in figure 12. First, on all days estimates for tuning parameters on the same electrode move up and down together. This immediately follows from the modeling assumption that tuning parameters on a given electrode are yoked to a single base value, b_{d, e}, in equation (5). Second, tuning parameter estimates vary more strongly at the beginning than at the end of a decoding session. This can be observed in the steep change in tuning parameter values at the beginning of certain days, such as day 14 for electrode 17 in figure 12. Intuitively, the reason for the reduced variability in tuning parameter estimates for both decoders is that confidence in estimates of base values increases as more neural activity is observed. Mathematically, this corresponds to less uncertainty in the distribution of base values, $P(\boldsymbol{b}_d | \lbrace \boldsymbol{x}_{d,s} \rbrace _{s=1}^t)$, for the SR classifier and an increasing value of n_{d, t} for the SR_S classifier.

Zoom In Zoom Out Reset image size

The 'spikes' in the parameter estimates for an electrode produced by the SR classifier visible in both figures occur when that electrode is flagged as behaving erratically by the mechanism for detecting anomalous spike counts. When this occurs, the uncertainty around the base value for that electrode is increased. This allows spike counts collected immediately after an electrode is flagged to strongly influence parameter estimates, and just as at the beginning of a day, multiple trials have to pass before the uncertainty for the base value of a flagged electrode is reduced and parameter estimates vary less. The reader may notice that the plotted electrodes are flagged as behaving erratically multiple times within the same day. As explained in section 2.5.3, the false positive rate of the flagging mechanism was set to 1%; therefore, the frequency of parameter resets seen in figure 12 should not be understood to indicate the true frequency of anomalous events. While we do not show results for this, in practice the SR classifier appears to be reasonably robust without this mechanism to detect and compensate for electrodes with anomalous spike counts. In fact, classification accuracy of the SR classifier was markedly improved on only three days for monkey L and on no days for monkey I by its inclusion.