Intracortical recording stability in human brain–computer interface users

Two subjects were implanted with platinum-coated silicon intracortical microelectrode arrays (Blackrock Microsystems, Salt Lake City, UT) in left primary motor cortex (M1). Subject 1 was a 54 year old female diagnosed with spinocerebellar degeneration without cerebellar involvement [43]. She was implanted with two 4 mm  ×  4 mm, 96-channel arrays for 33 months. Subject 2 was a 28 year old male with a C5 motor/C6 sensory ASIA B spinal cord injury. He had two 4 mm  ×  4 mm  ×  1.5 mm, 88 channel arrays implanted in motor cortex. Data from the first six months of recording from Subject 2 are included here. Subject 2 also had two 4 mm  ×  2.4 mm  ×  1.5 mm, 32 channel arrays implanted in primary somatosensory cortex, however, the recordings from those arrays are not included in this analysis. The arrays in somatosensory cortex were implanted for the purpose of providing cutaneous sensations through intracortical microstimulation [44], and were not recorded from as consistently as the arrays in motor cortex.

Neural data were recorded during BCI test sessions, typically performed 2–3 d a week, using the NeuroPort data acquisition system (Blackrock Microsystems, Salt Lake City, UT). We used 3 min of data collected during closed-loop BCI decoder calibration in order to identify the units recorded for each test session. We recorded 287 sessions across 33 months from Subject 1, and 51 sessions across six months from Subject 2. The neural data recorded during calibration for each of these days were used to quantify between-day unit stability.

We also wanted to quantify stability within a single day, so on a subset of days, 5 min long recordings were taken during closed-loop BCI calibration approximately every hour for 7 h. Subject 1 completed five sessions of testing in order to quantify within-day unit stability. These sessions were completed 276, 297, 307, 466, and 467 d after implant. Subject 2 participated in four sessions of within-day stability testing at 101, 128, 163, and 193 d after implant.

For all recording sessions, voltages were sampled on each electrode at 30 kHz with a 250–4500 Hz bandpass filter. Each time the voltage signal crossed a pre-defined spike threshold, the time of the crossing and a 48-sample (1.6 ms) snippet of the signal, starting ten samples before the threshold crossing, was saved. The spike threshold was set to either  −5.25 (Subject 1 prior to day 565) or  −4.50 (all other test sessions) times the baseline root-mean-square voltage on each channel, and then rarely manually adjusted if necessary.

Threshold crossing events were later sorted offline in Matlab (MathWorks, Natick, MA) to identify 'units', which we defined here as a collection of threshold crossing events from a single electrode that have similar waveform shapes. These units could be single-neuron action potentials, multi-neuron activity without well-isolated action potential shapes, or a collection of indistinct low-amplitude threshold crossings. Our sorting method relied on principal component analysis (PCA) to separate units based on their waveform shape. For each recording session and channel, the researcher viewed two plots of neural data prior to sorting the units. The first plot showed up to 400 overlaid spike snippets from a single channel (figure 1(A)). The second plot showed all recorded snippets from that channel as points in the plane of the first two principal components based on snippet waveform shape (figure 1(B)). The researcher then specified the number of sortable units (including single-neuron, multi-neuron, and indistinct low-amplitude units) on a given channel and initialized a Gaussian mixture model expectation-maximization algorithm to create the specified number of clusters in the principal component space. Every recorded snippet was assigned to the cluster to which it had the highest probability of belonging according to the expectation-maximization algorithm. Once all snippets were assigned to a unit, any unit that had a firing rate less than 0.2 Hz was discarded prior to further analysis.

Zoom In Zoom Out Reset image size

Once the snippets were sorted into units, we used a method developed by Fraser and Schwartz to determine if the units were stably recorded from one session to the next [6] (posted to Matlab Central as 'Tracking neurons over multiple days', identification no. 30113). A summary of this method is presented here. The method used four classifying characteristics: waveform shape, auto-correlation, mean firing rate, and cross-correlation with other units. The differences in these characteristics between a unit and any unit recorded during a different session generates a 4D similarity score. These characteristics 'are individually weak classifiers, but, because they represent independent sources of information, they can be combined into a strong classifier' [6]. The method also makes use of the assumption that a neuron can only be recorded by a single electrode, and so neurons recorded on different electrodes must be unique. Based on this assumption, the similarity scores for each unit when compared to all units recorded on all the other electrodes is calculated to generate a 4D null distribution of similarity scores for units that are known to be different. Similarity scores were then calculated between sessions for units recorded on the same channel, creating a population of potentially-stable scores. Using expectation-maximization to fit a mixture of Gaussian models to the single-label null distribution of known-different units, and the similarity scores of the multi-label potentially-stable units allows a quadratic classifier to compute an optimal decision boundary with a set error rate. For this analysis, the error rate was set to 5%.

As a first estimate of stability, pairwise comparisons were made from recordings taken on the same electrode during adjacent testing sessions. Each test session was compared to the session that immediately followed it. We refer to this as the iterative-comparison stability method. Pairwise comparisons continued longitudinally across all recording sessions. If a unit was stable from the first session to the second, and then again from the second to the third, it was considered stable across all three sessions. This yielded a monotonically decreasing number of stable units across time. For the within-day stability analysis, iterative pairwise comparisons were made to the recording sessions taken each hour.

The iterative-comparison method provided a conservative estimate of the number of stable units across three or more sessions due to its iterative nature. A false negative result during any single session comparison marked a unit as unstable for the rest of the test period, even if it was labeled as stable in all session comparisons following the false negative. While false negatives could inflate estimates of instability, a complementary feature of the iterative approach was that false positives, which occur at a 5% rate (based on the decision boundary being set at 95% of the null distribution), were rapidly corrected after only a few comparisons. Therefore, the iterative-comparison method for determining stability provided a lower bound on the estimated number of stable units over the test period.

An alternative method to determine stability is to compare each session to all following sessions directly. This direct method of comparison allowed units that may not have been recorded during some sessions to be marked as stable if they were recorded again later, similar to the method used by Vaidya et al [8]. The direct-comparison method provided a more generous estimate of unit stability over time. With this method, each comparison has a 5% false positive rate. However, false negatives which incorrectly mark a unit as unstable, are limited to single pairwise session comparisons rather than every subsequent session. By using both the iterative-comparison method and the direct-comparison method, we estimated a range within which the true number of stable units fell.

To perform the direct-comparison, we compared each session directly to the subsequent 50 recording sessions, or all remaining sessions if there were fewer than 50, to identify stable units. We chose to only make direct-comparisons across 50 sessions at a time because that would always yield comparisons for at least 100 calendar days, and continuing the comparisons beyond that point provided minimal additional information because the shape of the stability decay curve was well defined by the first 100 d. Once those comparisons were done for each session, we had the stability rates for units from every recording session to all subsequent recording sessions. For the within-day comparison, the first recording in each day was directly compared to the subsequent seven recordings.

The following predictive characteristics were calculated using data from the first day that a given unit was recorded: peak-to-peak voltage, baseline noise, firing rate, firing rate consistency, and 'model-free tuning', a task-independent metric of a unit's firing rate correlations. These predictive characteristics were selected separately from the classifying characteristics that were used to determine stability between sessions, (see iterative-comparison stability determination section). The stability determination used individual unit characteristics as well as information about the entire recorded population of units to enable a robust characterization of unit identity where we assume that a unit can only be recorded on a single electrode tip. These same characteristics were used previously in non-human primate studies to track unit stability across days [6]. Here, we chose to compare units that we classified as having different levels of stability in terms of characteristics that could be computed for each unit individually rather than considering the relationship to properties of the recorded population. The predictive characteristics computed for each unit are defined as follows:

• 1.  
  Baseline noise was calculated as two times the standard deviation of the first five and last five samples of each 48-sample unit snippet. These samples were chosen because each snippet is aligned such that the first ten samples occur before the threshold is crossed, and most waveforms will have returned to baseline approximately 1 ms (30 samples) after crossing threshold. Units with high baseline noise might be expected to be less stable because the recording noise could overwhelm the signal.
• 2.  
  Peak-to-peak voltage was calculated as the difference in voltage from the trough to the peak of the average waveform for each unit.
• 3.  
  Firing rate was calculated as the average number of threshold crossings per second across the 3 or 5 min recording. Firing rate was square-root transformed to better approximate a normal distribution [45].
• 4.  
  Firing rate consistency was calculated as the inverse of the variance of the inter-spike interval divided by the mean firing rate of the unit. A unit with a low level of firing rate consistency might be more likely to be unstable if this is an indication that its firing behavior is more transient.
• 5.  
  Model-free tuning, which we defined as the percentage of firing rate variance explained through factor analysis. Typically, neural tuning is calculated by determining how much of the observed firing rate variance can be explained by task-relevant variables, such as position or velocity. However, since different BCI tasks, including 2 degree-of-freedom cursor control and 3, 5, 7, and 10 degree-of-freedom hand and arm control [2, 28, 31, 36], were conducted across days, we could not evaluate one specific model to quantify neural tuning. Instead, a 10D factor analysis model was used to identify the amount of shared firing rate variance for each unit as described in Williamson et al [32] and Bittner et al [46]. Factor analysis approximates the covariance matrix of the neural population activity as the sum of two components: a shared component, that represents shared fluctuations across the population, and a private component that represents variation in activity that is specific to each neuron. Briefly, if Σ is the n-by-n covariance matrix measured for all n neurons in a single session, Σ  ≈  Λ  +  Φ where Λ is a low-rank covariance matrix that captures the ten dimensions explaining the most shared covariance of the population and Φ is a diagonal matrix that summarizes the remaining 'private' variability not accounted for by shared fluctuations in the population activity. Essentially, Φ represents the independent variance of each neuron that is not explained by other recorded neurons. As a proxy for neural tuning across various BCI tasks, model-free tuning was calculated as the 'percent shared variance', i.e. the ratio of a neuron's shared variance to its total variance:

  $$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{\boldsymbol{m}}_{i}}=~\frac{{{\boldsymbol{\Lambda }}_{ii}}}{({{\boldsymbol{\Lambda }}_{ii}}+{{\boldsymbol{\Phi }}_{ii}})}\nonumber \end{align} \tag{ 1 }$$

  where m_i is the model-free tuning value for unit i and Λ_ii and Φ_ii are the scalar diagonal entries of their respective matrices corresponding to unit i. A 10D model was chosen to match the largest number of degrees-of-freedom decoded during the study [36]. Higher levels of model-free tuning indicate that a neuron's activity is better predicted by the population. All of the tasks during which the data was collected involved the subject using the recordings from M1 to control the velocity of an end-effector. We therefore assume that the majority of shared variance in the population will be related to the task, and that units with a higher percentage of shared variance are therefore better tuned to the task.

A Wilcoxon rank-sum test was used to compare each predictive characteristic between units that were stable for at least 1 d and units that were not stable between days. We also compared these characteristics for units that were stable for 1–6 d, 7–21 d, or greater than 21 d using the Kruskal–Wallis test.

Given that there may be an interaction between the different predictive characteristics and their impact on unit stability, we sought to model the relative utility of the different predictive characteristics for predicting the length of stability using Cox regression modeling [47]. The Cox regression model uses a survival function of the form:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle S\left(t,i \right)={{S}_{0}}{{\left(t \right)}^{{\rm exp}\left(\underset{k}{\mathop \sum }\,{{b}_{k}}({{C}_{k}}\left(i \right)-{{{\bar{C}}}_{k}}) \right)}}\nonumber \end{align} \tag{ 2 }$$

where $S(t,i)$ is the survival rate for unit i at time t, days after a selected recording session. ${{S}_{0}}(t)$ is the baseline survival rate at time t, the average survival rate across the entire population. The exponent of the equation is called the relative risk. The relative risk is calculated using the b regression coefficients for each characteristic: noise, peak-to-peak voltage, firing rate, firing rate consistency, and model-free tuning. The subscript k denotes each different characteristic. ${{C}_{k}}(i)$ is the measured characteristic of each unit, and ${{\overline{C}}_{k}}$ is the population mean for characteristic k. The higher the survival rate for a unit at a given time after the initial recording session, the more likely it is to still be stable. The model was fit based on each observed unit's predictive characteristics on the first day it was observed, and the number of days it was subsequently observed before becoming unstable. The model was only fit to Subject 1's data because the time period of recordings for Subject 2 was insufficient to confidently fit a model.

The distributions of all of the predictive characteristics were heavily skewed. Data were log transformed prior to performing the Cox regression analysis in order normalize the distribution of observed values.

We also sought to examine whether prior stability was predictive of future stability. For example, was a unit that had been stable for a week more likely to be stable for another week than a unit that had only been recorded for one day? Units were labeled by the number of days that they had been stable, ranging from 1 to 60 d. For each duration of stability (n  =  1–60 d), we determined how many days longer each unit would remain stable. In this way, we were able to compare the expected length of stability for units based on how long they had already been recorded. This analysis was performed only for Subject 1, as Subject 2's data was collected over a shorter period of time. Additionally, we chose to investigate prior stability separately, rather than as part of the Cox regression modeling because prior stability was highly correlated with the predictive characteristics.

The largest number of recorded units (defined here as a mix of well-isolated single neurons, multi-neuron waveforms, and low-amplitude indistinct threshold crossings) was observed in the first two months after implant for both subjects. Subject 1 had a maximum of 301 recorded units on Day 17 and Subject 2 had a maximum of 407 units on Day 52. Since the participants had a different number of wired electrodes, we computed a ratio of the number of recorded units per channel, as shown in figure 2(A). This ratio peaked at 1.6 units per channel for Subject 1 and 2.3 units per channel for Subject 2. The number of sorted units declined over time, consistent with previous non-human primate and human subject studies. While the yield varies widely across subjects, studies tend to report many fewer units in implants that are more than a year old, compared to implants that are only a few months old [30, 48–50].

Zoom In Zoom Out Reset image size

The peak-to-peak voltage of the recorded units increased over the first month of the implant before slowly decreasing (figure 2(B)), again consistent with previous literature [30, 48–50]. The median unit size reached a maximum value of 89 and 102 µV, for Subjects 1 and 2 respectively. Note that the recording quality metrics in figure 2 exhibit a discontinuity at day 565 for Subject 1 due to a change in the spike threshold from  −5.25 to  −4.5 RMS. By moving the threshold closer to 0, more units with lower peak-to-peak voltages were recorded. This change was motivated by reports that small amplitude units can carry significant movement-related information and contribute to decoding [34, 35]. Subject 2 had more units per recording channel compared to Subject 1. This could be attributed to using the more sensitive threshold (−4.5 RMS) for the duration of the Subject 2's implant. Subject 2 tended to have more large amplitude units, as seen by the higher bound on the interquartile range for peak-to-peak voltage in figure 2(B).

These recording quality results demonstrate that the arrays performed as expected based on previous results [30, 48–50]. We do not believe this slow decline in recording quality is directly related to the recording instability observed over the course of hours and days described below. If units were lost from the recorded population because they became permanently unrecordable, then the number of available units would decline much faster than we observed.

Data were compiled across all recording sessions to obtain the median trajectory for unit stability over the subsequent days and weeks of testing. Figure 3 shows the median and interquartile ranges for estimated unit stability rates using both the iterative (figure 3(A)) and direct (figure 3(B)) comparison methods. In both cases, we observed two primary responses: (1) a sharp initial decline in the number of stable units and (2) a small group of long-lasting stable units. Therefore we fit the time course of unit loss with a two-term exponential function. This implies that units do not decay at a constant rate. Instead, there is a large proportion of units (~75%) that were moderately stable over multiple days but became unstable within one week. The remaining units were highly stable and were lost at a slower rate. The exponential function used to describe unit stability for both the iterative- and direct-comparison methods was:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {\rm Proportion}\,{\rm of}\,{\rm Units}\,{\rm Still}\,{\rm Stable}=a{{{\rm e}}^{-bt}}+c{{{\rm e}}^{-dt}}\nonumber \end{align} \tag{ 3 }$$

where t is the number of days since the initial recording session. We define this equation such that a is the proportion of the recorded population that were moderately stable, c is the proportion of the recorded population that were highly stable, $(1-{{e}^{-b}})$ was the proportion of moderately stable units that become unstable each day, and $(1-{{e}^{-d}})$. was the proportion of the highly stable units that become unstable each day. The proportion of units that were unstable is the remainder of the population not described by equation (3) when t  =  0, i.e. 1  −  (a  +  c).

Zoom In Zoom Out Reset image size

Figure 3(C) shows the two-term exponential models of unit stability for Subjects 1 and 2 using both the iterative- and direct-comparison methods. The participants exhibited very similar trends in the rate of unit stability.

Using the more conservative iterative-comparison method 47.1% of recorded units from Subject 1 were unstable. However, 44.7% of units recorded from Subject 1 showed moderate stability. Of these moderately stable units, 16.3% became unstable each day. Finally, 8.2% of units showed high stability, with only 2.6% of these units becoming unstable each day. Subject 2 showed similar trends although a larger percentage of units (61.4%) were unstable. Another 32% of units showed moderate stability, with 15.3% of those becoming unstable each day. Another 6.6% of units showed high stability, with 3.7% of those becoming unstable each day. In figure 3(C), the y-intercepts of the solid lines are equal to the sum of the percentage of units that are moderately stable and the percentage of units that are highly stable for each subject.

As expected using the direct-comparison approach, a larger percentage of units were considered to be highly stable (19.7% for Subject 1 and 12.5% for Subject 2). Nearly 20% of units were moderately stable for both subjects. Using the direct-comparison approach, we estimated that approximately two-thirds of units were unstable. The estimated decay rates were lower using the direct estimation of stability as compared to the iterative approach, which is much more conservative across repeated sessions (table 1). All data in table 1 is derived from the two-term exponential models of the stability estimates.

The between-day stability analysis revealed that approximately 60% of units were only recorded on a single day (table 1). As such, we sought to examine unit stability over a shorter timescale by collecting eight recordings, one every hour for 7 h, to identify units that were stable on the order of hours, rather than days.

For both subjects, the iterative- and direct-comparison methods estimated that 72%–80% of recorded units were stable across two consecutive recordings, separated by 1 h (table 2). Figure 4 shows the unit stability profiles for the iterative (4(A)) and direct (4(B)) comparison methods. Overall, Subject 1 appeared to have a slightly higher percentage of stable units, which is consistent with the between-day results. Both subjects showed the expected trend, in that fewer units were recorded stably as the time since the initial recording increased.

Zoom In Zoom Out Reset image size

We fit the within-day stability estimates with a single-term exponential function as illustrated in figure 4(C) because the collected data did not extend far enough in time to judge whether highly stable units had a significantly different decay rate as was observed between days. The exponential regression models were used to estimate the rate at which stable units could no longer be recorded. Using the more conservative iterative-comparison approach, we estimated that stable units were lost at a rate of 8.7% or 12.5% per hour for Subjects 1 and 2 respectively. When we instead used the direct-comparison method, we estimated that the number of stable units decreased at a rate of 3.6% or 5.8% per hour for Subjects 1 and 2 (table 2). These decay rates are markedly faster than those from the between-day stability analysis. The shorter length of data collection makes it impossible to properly identify units that will be highly stable, leading to the assumption that all units will become unstable as quickly as the ones that are observed to become unstable during the 8 h of recording. This is why it is important to consider the between-day and within-day analyses together.

Data from the between- and within-day experiments can be combined to provide a rough estimate of overall unit stability. Using the iterative-comparison method, approximately 25% of units were stable on the order of hours, 47.5% of units were stable on the order of days, and 7.5% were stable on the order of weeks. Using the direct-comparison method, approximately 45% of units were stable for hours, 17.5% for days, and 12.5% for weeks.

When each predictive characteristic was considered individually, unstable units had lower noise levels (12.5 versus 14.0 µV and 23.9 versus 25.8 µV for Subjects 1 and 2 respectively), lower peak-to-peak voltages (37.6 versus 43.9 µV and 62.5 versus 75.4 µV), lower firing rates (0.36 versus 1.59 Hz and 1.84 versus 3.52 Hz), less firing rate consistency (0.47 versus 0.67 and 0.74 versus 0.99), and lower model-free tuning (0.010 versus 0.022 and 0.016 versus 0.022) than units that were stable for at least one day (p  <  10⁻¹⁸ for all predictive characteristics for both subjects, Wilcoxon rank-sum). When comparing across different durations of stability (less than 1 week, 1–3 weeks, or greater than 3 weeks), we observed differences in all 5 predictive characteristics in Subject 1 (p  <  10⁻¹⁶, Kruskall-Wallace, figure 5). Subject 2 showed a significant difference (p  <  0.01) for all predictive characteristics other than firing rate consistency (p  =  0.32) for units that were stable for the same three time periods.

Zoom In Zoom Out Reset image size

We created a Cox regression model to predict the survival rate of a given unit over time based on its noise level, peak-to-peak voltage, firing rate, firing rate consistency, and model-free tuning. For both subjects, higher peak-to-peak voltage, firing rate and model-free tuning all correlated with a higher survival rate. Firing rate had the strongest correlation, followed by peak-to-peak voltage, and then model-free tuning. Firing rate consistency did not have a significant correlation with survival rate. Baseline noise had a strong correlation with survival rate in the data from Subject 1 but not Subject 2. This could be because all of Subject 2's data were collected early in the implant when signal quality is generally higher. The Cox regression model allows us to assess the relative risk that a unit will become unstable. The relative risk is the exponent from the survival function (equation (1)), which is calculated from observed predictive characteristic values and the regression coefficients. An increase in relative risk decreases the survival rate of a unit. To illustrate the relative strength of each characteristic we calculated the relative risk for hypothetical units. For the characteristic of interest, we evaluated the 25th or 75th percentile value measured across all units, while the other characteristics were set at the median value. The change in relative risk between the 25th and 75th percentile hypothetical unit for each characteristic is shown in table 3.

The Cox regression model also allows us to estimate the survival rates for the 'best' and 'worst' units. Here we defined the best-case unit as one in the 95th percentile for firing rate, peak-to-peak voltage, and model-free tuning, and in the 5th percentile for noise and firing rate consistency. Conversely, the worst-case unit was one in the 5th percentile for firing rate, peak-to-peak voltage, and model-free tuning, and in the 95th percentile for noise and firing rate consistency. By plotting the predicted stability for these two cases and for a unit with the median value in all predictive characteristics, we can visualize the range of possible survival rates (figure 6). The best-case unit had a relative risk of becoming unstable that was 3.7×  lower than the worst-case unit and 2.1×  lower than a unit with the median value of each characteristic. This means a best-case unit has a 36.8% chance of remaining stable for one week, compared to 12.8% for a median unit and only 2.4% for a worst-case unit. The best-case unit is also much more likely to be stable for one month (9.5%) compared to the worst-case unit (0.015%).

Zoom In Zoom Out Reset image size

For Subject 1, we could also examine whether the amount of time that a unit had been stable predicted future stability. As shown in figure 7, the longer a unit had been stable, up to a point, the longer it was likely to remain stable in the future. We see a steady increase in both the amount of time a unit will remain stable and the variance about that estimate, up to approximately 40 d of stable recording. After 40 d there appears to be a plateau in the amount of time a unit is expected to remain stable.

Zoom In Zoom Out Reset image size

We found that approximately 60% of units are unstable after a single day. Even within a single day, we found that 45% of units became unstable after 7 h. The rate at which units became unstable was well-fit by a one-term exponential for recordings made within a day, and by a two-term exponential for recordings made on multiple days. Within a single day, units became unstable at a rate of approximately 11% per hour. For units that were stable for at least one day, the majority were lost at a rate of 16% per day, while a smaller percentage of units were highly stable for multiple weeks and were lost at a rate of 3% per day. Based on the differences in decay rates estimated by the two-term exponential models, we classified units as unstable, moderately stable, and highly stable. This is a coarse classification of the population that likely has a more continuous distribution of stability rates, however we are limited by the frequency (e.g. three times per week) of our recordings. In the future, continuous wireless [52, 53] recordings will likely make it possible to describe more detailed temporal profiles of unit stability.

Our results confirm prior reports [6–9] of signal instabilities that are commonly observed with intracortical recordings of rhesus macaques. When using the iterative-comparison method, we found slightly fewer stable units than Fraser et al who used the same method to quantify stability in two macaques [6]. However, our direct-comparison estimates of stability were comparable to what Vaidya et al observed in their macaque subject [8]. An important difference is that both studies showed substantially higher rates of stability over the first two weeks than we observed in our study, but this may be due to the fact that the previous studies only included well-isolated units. Given the similarities in long-term stability estimates between previous studies and this one, a macaque model appears to be appropriate for investigating issues surrounding chronic recording stability. However, further investigation is needed to determine if short term stability, particularly within a single day, is comparable between macaques and humans.

As noted above, some previous studies of stability in macaques [6, 8] only included well-isolated units. However, we chose to include poorly-isolated units, in addition to well-isolated units, because they contain significant information about movement parameters [34, 35] and are often used for decoding [3, 36]. Poorly-isolated unit recordings are often a summation of activity from multiple units and are therefore referred to as multi-unit activity. We recognize that differences between single- and multi-unit activity may be important for interpretation of our stability results. For example, the voltage measured from multi-unit activity, is more variable than that measured from well-isolated single units [35]. Additionally, well-isolated units tend to provide information about multi-dimensional movement velocities, while multi-unit activity may be more representative of movement speed [34]. The Cox regression model showed that units with lower peak-to-peak voltage, higher noise levels, and less distinct tuning were less stable overall (see worst-case curve in figure 6). These predictive characteristics are more typically observed for multi-unit recordings suggesting that multi-unit activity is more unstable than single-unit activity. This result will need to be reconciled with the relatively high levels of stability in LFP and unsorted threshold crossings found by Flint et al [17]. One potential explanation is that by averaging activity over a much larger population for LFPs, or the combination of well-isolated single-units and multi-unit activity that contribute to threshold crossings, individual neuron stability becomes less important.

To be sure that our results were not driven primarily by poorly-isolated units, we conducted a secondary analysis using the within-day stability data. With this analysis, we defined well-isolated units as those with the most consistent waveform shapes. Consistency was calculated by taking the dot product of each recorded snippet with the subsequent one after normalizing them to remove the effect of peak-to-peak voltage. This allows us to identify units with very consistently shaped waveforms (high average dot product). The population of units was separated based on a mixture of Gaussians model, with 47.7% and 54.1% of units being labeled well-isolated for Subjects 1 and 2 respectively, and the rest being labeled poorly-isolated. While waveform similarity is a criteria used to determine stability, stability determination is based on the comparison of mean waveform shape between sessions, while the isolation determination was made based on the variability of the waveform from one spike to the next within a single recording session. After 1 h, 81% of well-isolated units were stable as compared to 63% of poorly-isolated units. Approximately 8.5% of well-isolated units became unstable each hour as compared to 16.4% of the poorly-isolated units. This result agrees with the between-day analysis of the Cox regression model and suggests that multi-unit activity is less stable, although it is informative for BCI. Developers of long-term BCI systems will need to balance the need for the extra information available from multi-unit activity [34, 35] with the need for recorded units to remain stable for long periods of time.

In addition to measuring the rate of unit instability, we also sought to determine whether there were particular characteristics that distinguished stable units from those units that were more likely to be unstable. We investigated unit recording noise, peak-to-peak voltage, firing rate, firing rate consistency, and model-free tuning because they provide information about a unit's activity, consistency, clarity, and modulation and are independent of the task being performed. In general, stable units tended to have a higher firing rate, peak-to-peak voltage, and model-free tuning value. Units that were not stable had lower values of these characteristics, likely because many of them were low-amplitude activity without clearly defined action potentials. We used a regression model to predict unit stability and found higher firing rates were most predictive of prolonged stability with higher peak-to-peak voltage and model-free tuning values also providing significant improvements in expected length of stability. These characteristics appear to be the most robust predictive characteristics that could be used to identify units that would provide more stable long-term recordings, creating more reliable decoders.

The within- and between-day stability time courses provide some insight into the causes of recording instability. The typical explanation of recording instability has been that the electrode tip moves relative to the recorded neurons so that the distance becomes too large for clear recording [11, 12]. Our within-day analysis shows that the loss of units is gradual over this time period suggesting that electrode movement does not occur suddenly, for instance from sudden head movement. It remains a possibility that micromotion occurs at the electrode-tissue interface due to ongoing biological processes [54]. Another possible explanation for instability is that neurons go through active and inactive states, not unlike the inactivation that occurs during slow wave sleep [55]. A change in unit state could drastically alter the firing rate causing the unit to be recognized as unstable. Finally, neuronal death and glial scarring have also been proposed as potential causes of unit instability [54, 56]. With glial scarring, neurons are further away from the electrode tip and at some point could become too far away to be recorded. However, we observed that the total number of recorded neurons decreases much more slowly, by orders of magnitude, than the number of stable units so we do not expect that neuronal death or an increase in glial scar thickness are the primary cause of instability (figures 2 and 3).

We found the rate of stability was significantly increased when estimated with the direct-comparison method suggesting units frequently return after being marked as unstable. This is consistent with the theories of micromotion and state changes. Neurons that move away from electrodes could move back towards the tip if movement resembles a random walk, so they would be recognized as stable with the direct-comparison method, but not the iterative-comparison method. Similarly, neurons would be recognized as stable with the direct-comparison method if they returned to an active state after entering an inactive state. Further evidence to disentangle these hypotheses could be gathered using forthcoming wireless recording systems [52, 53] to conduct continuous long-term recordings. If most units become unstable after small changes in peak-to-peak voltage, this would suggest that micromotion may be the primary driver of instability. In contrast, sudden shifts in firing rates would suggest a change in neural state.

If micromotion of the electrode tips relative to the recorded neurons were the primary driver of instability, then the units with large peak-to-peak voltage would be the most stable because the recorded neuron is nearer to the electrode tip and could tolerate a small amount of motion while still being recorded. We did observe that a large peak-to-peak voltage was predictive of stability, however we also found that firing rate and model-free tuning were predictive of stability. There is not an obvious reason for firing rate and model-free tuning to correlate with a neuron's proximity to the electrode. We think that this suggests neuron-type may explain differences in stability rate, potentially because certain neuron-types undergo state changes differently [55] or because certain types of neurons are more or less stable for reasons that are yet to be determined. A thorough waveform analysis, for instance using the fact that interneurons have narrower spike waveforms than pyramidal neurons [57], may be able to determine whether specific types of neurons are more consistently recorded.

While the stability of recordings may improve with changes to the size, shape, materials, and coatings of recording electrodes [56], understanding the rate at which units are changing and which predictive characteristics are more likely to be stable using existing clinical electrodes provides useful information for the development of intracortical BCI devices. One important design criterion for a clinical BCI is to minimize the number of times that the device needs to be recalibrated. One possible approach is to identify very stable units and use only those units for the BCI. For example, Ganguly and Carmena [33] showed that a monkey could complete a cursor task using the same decoder for up to 19 d using just 15 recorded units. The results from the iterative-comparison analysis provide a conservative estimate that a 96-channel array (assuming ~1 unit per channel, figure 1(A)) would provide at least 15 stable units for 9 d. The direct-comparison analysis, predicts at least 15 stable units for 29 d. We would therefore expect that if we could identify these stable units, BCI performance on a cursor control task could be maintained for 9–29 d.

A significant limitation of selecting only highly stable units for a BCI is that it will not scale well to applications that require higher degrees of freedom because a limited amount of information can be extracted from such a small sample of units. However, the estimates of unit stability presented here should aid in the development of unsupervised recalibrating decoders [37–42]. These decoders adjust to the loss of units from the recorded population and the addition of the new units that replace them by recalibrating the decoder during use. The stability decay rates shown here provide recalibration algorithms with an a priori estimate of how quickly units will be lost and trigger a recalibration before performance decays too far. The results from the within-day comparisons are most pertinent here, showing that around half of the recorded units turn over within an 8 h period. This indicates that updates will need to occur at least every few hours, which will require a system that is used consistently through the day. The units that are identified as stable could also provide a baseline framework to which the activity of newer units can be compared to begin estimating their tuning. If unsupervised recalibration could run in the background of a BCI system, then a BCI user may only have to participate in supervised decoder training once while maintaining performance for at least the rest of the day without operator intervention.

An important additional consideration is that while a unit may remain stable, it is possible that its relationship to motor output may change over time. For example, some studies have shown that preferred directions and modulation depth can change with learning or task context as the subject's control of the BCI improves [15, 16, 58]. We did not test this explicitly in our study, however, mean firing rate and model-free tuning are likely related to a neuron's functional properties. Future studies will need to characterize stability in terms of BCI performance across different types of tasks. Here we provide information regarding the stability of the ability to record a given unit chronically in a human intracortical implant.