Artifact-free recordings in human bidirectional brain–computer interfaces

A 28 year old male participant with a C5 motor/C6 sensory ASIA B spinal cord injury sustained 10 years prior was enrolled in the clinical trial. Presurgical magnetoencephalography imaging, as described previously [12], was used to identify the hand areas of primary motor cortex (M1) and area 1 of somatosensory cortex (S1). Four Utah microelectrode arrays (Blackrock Microsystems, Salt Lake City, Utah) were implanted in the participant's left sensorimotor cortex, as shown in figure 1. Two of the arrays were designed exclusively for recording and were implanted in the hand and shoulder areas of M1. The other two arrays were designed primarily for stimulation and were implanted in the hand area of S1. Each recording array was 4 mm  ×  4 mm in size with 88 functional electrodes with platinum-coated tips. Each stimulation array was 2.4 mm  ×  4 mm in size with 32 functional electrodes coated with a sputtered iridium oxide film, which increased the charge injection capacity for stimulation. The arrays were wired to two percutaneous connectors attached to the skull. It was not possible to record from the S1 arrays during ICMS, and the goals of our study do not necessitate simultaneously stimulating and recording from the same array.

Zoom In Zoom Out Reset image size

This study was performed under an Investigational Device Exemption granted by the US Food and Drug Administration. The study received approval from the Institutional Review Boards at the University of Pittsburgh and the Space and Naval Warfare Systems Center Pacific and is registered at ClinicalTrials.gov (NCT01894802). Informed consent was obtained prior to performing any experimental procedures.

Neural recording was accomplished using two Neuroport clinical electrophysiology systems, each consisting of a neural signal processor (NSP) and Front-End Amplifier (Blackrock Microsystems, Salt Lake City, UT). All system components are shown in figure 2. The amplifier firmware was modified to add a digital sample-and-hold feature, which allowed the digital outputs on each channel of the amplifier to be fixed at their last sampled value when the voltage on a TTL input on the amplifier was set high. By setting this TTL input high during stimulus pulses, we could effectively 'blank' out the stimulus artifact prior to subsequent signal processing onboard the NSP, which was unmodified and could only perform standard signal processing operations for spike extraction. Spike extraction was performed on the NSP by high-pass filtering and thresholding the digitized and blanked signal to detect neural spikes.

Zoom In Zoom Out Reset image size

ICMS was delivered using the Cerestim R96 microstimulator (Blackrock Microsystems). The stimulator has a 'sync' output, which produces a TTL pulse beginning 60 µs before each stimulus pulse. A monostable multivibrator (Texas Instruments CD74HC123, Dallas, TX) 'blanking circuit' was triggered by the rising edge of the sync signal and generated a TTL pulse of an arbitrary duration set by a potentiometer. This TTL blanking signal was then connected to the amplifier to blank the recordings during stimulation pulses. A timing diagram of this process is shown in figure 3. The Cerestim R96 and NSPs interfaced with a custom BCI software suite, developed in Matlab and C++, and which was configured to ignore registered spikes that occurred in the first sample after the offset of the sample-and-hold blanking period. These 'spikes' were the result of discontinuous sampling of the analog signal.

Zoom In Zoom Out Reset image size

ICMS consisted of trains of charge-balanced, cathodic-leading asynchronous pulses delivered at 100 Hz. Each pulse had a 200 µs cathodal phase duration, 100 µs interphase duration, and 400 µs anodal phase duration. Cathodal amplitudes ranged from 2 to 100 µA; anodal amplitudes were always half of the cathodal amplitude. Up to twelve electrodes could be simultaneously stimulated, with a maximum total charge per phase of 144 nC.

The system presented above has two primary components that can be tuned to reject stimulation artifacts while preserving neural data. First, the digital filter used for spike extraction can be tuned to minimize the effect of stimulation artifacts. Second, the raw signals can be blanked during stimulation for an arbitrary duration. This blanking duration can be tuned to balance removing the artifact from the raw signal, while also preserving neural data. These features were tuned using raw voltage recordings collected from electrodes in M1 when ICMS was applied to S1.

Electrode recordings from both M1 arrays were collected during trials in which ICMS was delivered to the S1 electrodes for the purpose of tuning and validating the artifact rejection system. Analog voltage recordings were band-pass filtered (0.3–7500 Hz) and digitized at 30 000 samples per second. Signal blanking and digital filters were not applied online. This allowed for the artifact rejection and spike detection features to be implemented and tuned offline on the raw, artifact-contaminated signals.

The final online artifact rejection implementation was validated by analyzing the distribution of spikes detected in the inter-stimulus interval. Spike times were detected and recorded online from both recording arrays while the participant performed a motor imagery task consisting of reaching and grasping movements that modulated neural activity in M1 [2]. While the participant performed this task, we collected three 60 s long trials under three conditions: no stimulation, stimulation without artifact rejection, and stimulation with artifact rejection. The participant was informed that the stimulation was irrelevant to his motor task. During trials with stimulation, three electrodes on the medial stimulation array delivered 60 µA pulses at 100 Hz in 500 ms bursts. Each burst was separated by 500 ms. The spikes during each stimulation train were binned into 0.5 ms increments within the 10 ms interval between stimulus pulses. A histogram of the mean binned spike counts with 95% confidence intervals was constructed for each stimulation condition. Stimulation times from the trials with artifact rejection were then applied to the trials without stimulation to generate equivalent histograms for the no stimulation condition.

Digital filters affect both the artifacts caused by ICMS as well as the neural signals of interest. In a typical neural recording setup, a digital fourth-order Butterworth filter with a high-pass cutoff frequency of 250–300 Hz is used to remove low frequency components of the recorded signal [21, 22]. To study the effects of cut-off frequency on the ability to record single-units, a filter analysis based on Lempka et al [23] was performed on data recorded from M1 in the absence of stimulation. Raw voltage signals were recorded from both M1 arrays during trials lasting 45 s in five sessions across a 9 month period, with at least two months in between each session. We identified high-amplitude single units by selecting channels with more than 30 samples that crossed a  −150 µV threshold after filtering with a fourth-order high-pass Butterworth filter with a 250 Hz cutoff frequency [1]. A total of 58 channels across the five sessions were identified using this method and all sorted single units and multiunit signals from these channels were included in order to have a mixture of high- and low-SNR signals. A snippet detection threshold was set at  −5.25 times the root-mean-square voltage (V_RMS) of each filtered signal. The  −5.25 V_RMS multiplier has been used in previous studies (e.g. Collinger et al [1]) to extract units of varying SNR without being overly sensitive to low-amplitude noise. A biased estimate of V_RMS was calculated using the Blackrock Microsystems algorithm described in Christie et al [24] to reduce the effect of outliers. The extracted spike snippets were then sorted using principal component analysis and a Gaussian mixture model expectation-maximization algorithm [25]. Spike times were saved for each sorted unit. Units with less than five spikes were excluded from further analysis. A total of 167 sorted units across the 58 channels were analyzed further, which included both high-SNR single units and low-SNR multi-unit activity. The raw voltage signals were then re-filtered with a first-order high-pass Butterworth filter with varying cutoff frequencies. Sorted snippets were extracted for each filter using the previously determined spike times. The mean peak-to-peak voltage, noise estimate, and signal-to-noise ratio (SNR) were calculated for each unit and filter setting. The noise estimate was calculated for each channel as two times the standard deviation of the continuous filtered signal after removing each spike [23]. The SNR was calculated as the mean peak-to-peak voltage divided by the noise estimate. Unsorted threshold crossings were extracted for each cutoff frequency, using a threshold of  −4.5 V_RMS, to determine the effect of changing the cutoff frequency on the number of threshold crossings, as used for online BCI control.

The recording and stimulation electrophysiology system described above was used to develop a bidirectional BCI system that was based upon the high degree of freedom BCI system that our group has described previously [1, 2]. Multichannel spiking activity was recorded using the digital filters and thresholding techniques described above, but implemented in real-time on the NSPs. No spike sorting was applied; threshold crossings on each recording channel were treated as a single-unit. Firing rates were estimated using 20 ms bins of spike counts, which were then smoothed using an exponential function with a 440 ms sliding window. The two-phase calibration procedure previously described by Collinger et al [1] was used to train an optimal linear estimator (OLE) decoder to control the Modular Prosthetic Limb (Johns Hopkins University Applied Physics Lab). The resulting decoder predicts five degree of freedom (DOF) limb velocity by relating the square-root transformed firing rate, f, of each unit to 5D velocity, v, using equation (1).

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle f={{b}_{0}}+{{b}_{x}}{{v}_{x}}+{{b}_{y}}{{v}_{y}}+{{b}_{z}}{{v}_{z}}+{{b}_{r}}{{v}_{r}}+{{b}_{g}}{{v}_{g}}.\nonumber \end{align} \tag{ 1 }$$

The 5D velocity includes the Cartesian endpoint components x, y, and z, an endpoint roll orientation component, r, and a grasp component, g. Optimal linear estimation with ridge regression was used to determine decoding weights from each b coefficient to predict arm kinematics from the recorded firing rates. During testing, firing rates were scaled using equation (2) prior to decoding to account for global changes in firing rate [26]. Individual unit firing rates, f, were scaled by the mean population firing rate during training, F_training, and divided by the mean population firing rate over the previous 300 ms, F_{300 ms}, to produce a scaled firing rate, f_scaled.

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{f}_{{\rm scaled}}}=f\cdot \frac{{\rm mean}\left( {{F}_{{\rm training}}} \right)}{{\rm mean}({{F}_{300\,{\rm ms}}})}.\nonumber \end{align} \tag{ 2 }$$

To generate stimulation commands, torque sensors, located in the motors at the base of the digits of the Modular Prosthetic Limb, were mapped to electrodes in S1 [12]. Torque sensor values, T, were linearly converted to stimulation amplitude, A (see equation (3)), such that electrodes began stimulating at a minimum amplitude, A_min, when the corresponding sensor value crossed a defined threshold, T_min. The stimulation amplitude then increased linearly with increasing torque until a maximum value, A_max, was reached at a defined maximum torque value, T_max. Stimulus frequency was fixed at 100 Hz. Pulse phase durations were the same as those described previously. A 50% duty cycle was applied to limit stimulation durations to 1500 pulses within 30 s, i.e. after 15 s of continuous 100 Hz stimulation electrodes were temporarily disabled for 15 s.

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle A=\left( \frac{T-{{T}_{{\rm min}}}}{{{T}_{{\rm max}}}-{{T}_{{\rm min}}}} \right)\left( {{A}_{{\rm max}}}-{{A}_{{\rm min}}} \right)+{{A}_{{\rm min}}};{{T}_{{\rm min}}}\leqslant T\leqslant {{T}_{{\rm max}}}.\nonumber \end{align} \tag{ 3 }$$

The performance of the bidirectional system was evaluated using an object transfer task. The goal of this experiment was to demonstrate that the artifact rejection system enabled bidirectional BCI performance that was functionally equivalent to BCI performance without stimulation. In each session, a five DOF velocity decoder was trained and the decoder's functional performance was then evaluated in a task where the participant was required to use the Modular Prosthetic Limb to grasp a cylinder on the left side of a table, pick it up, and transfer it to the right side of the table. A 20 cm wide region in the center of the table was indicated by two lines of tape. The participant was required to hold the object above this region while transporting it, rather than simply drag the object across the surface of the table. If the object contacted the table at any point in this region, the participant was required to bring the object back to the left to reset the transfer. The task was scored based on the number of successful transfers in 2 min. In each of two sessions, the system was configured with either a fourth-order high-pass Butterworth filter with a cutoff frequency of 250 Hz or a first-order high-pass Butterworth filter with a cutoff frequency of 750 Hz. Filter settings, spike thresholds, and decoders were fixed throughout each session. Within each session, three task conditions were tested: no stimulation, stimulation with signal blanking, and stimulation without signal blanking. Here, signal blanking refers to the combination of blanking before digital filtering and thresholding, and software rejection to discard transient threshold crossings immediately after the blanking period. Each task condition was tested with both filter settings, for a total of six conditions. Each condition was tested in a block format with three repetitions per block. In all cases, coincidence detection, a standard feature of the NSP, was enabled to reject threshold crossings that occurred simultaneously on at least 45 channels. During task conditions with ICMS, eight electrodes spanning both S1 arrays were mapped to the index finger torque sensor and simultaneously amplitude-modulated using the sensor-stimulus transform in equation (3). For this task, T_min  =  .05 Nm, T_max  =  .25 Nm, A_min  =  18 µA, and A_max  =  46 µA. These parameters were chosen with the intention of generating large ICMS artifacts on both M1 arrays during the task and were not necessarily intended to provide beneficial feedback for the task. The 50% stimulation duty cycle, a required safety feature of our system, allowed for the possibility of intervals in which index finger torque did not generate stimulus pulses during ICMS trials. For example, if stimulation was delivered for 15 continuous seconds, stimulation would automatically shut off for 15 s. This could result in periods of grasp without associated ICMS pulses.

The stimulation artifacts recorded in M1 during microstimulation of S1 were typically on the scale of millivolts, and often exceeded eight millivolts, saturating the amplifier. These artifacts were one to two orders of magnitude larger than most extracellular potentials, such as those shown in figure 4. However, even in the case of amplifier saturation, spikes could be observed in the raw voltage signal less than one millisecond after the offset of a microstimulation pulse. Therefore, we determined that it would be possible to discard artifacts and recover spikes a short time after each stimulus pulse. Our artifact rejection scheme is the result of applying a minimal number of changes to our existing clinical BCI system. Signal-blanking duration and digital filter parameters were tuned to allow us to reliably record artifact-free spikes between each stimulation pulse.

Zoom In Zoom Out Reset image size

A typical neural recording scheme involves recording a high bandwidth voltage signal (i.e. 0.3–7500 Hz bandpass-filtered signal sampled at 30 000 samples per second), applying a digital 'spike filter', and then applying a voltage threshold to the filtered signal to detect spikes. A good spike filter is typically designed to remove low-frequency local field potentials (LFPs) while passing high-frequency action potentials. For example, in previous studies we used a fourth-order Butterworth high-pass filter with a 250 Hz cutoff frequency [1, 2]. However, ICMS introduces another source of high-frequency content that can cause these standard filters to produce undesirable output after the stimulus pulse is completed. We wanted to explore whether modified filters alone could minimize these effects.

Microstimulation generates high-amplitude stimulus artifacts that are impulse-like in nature. As a result, these high-amplitude artifacts generate filter output resembling the filter's impulse response, distorting the artifact. For infinite impulse response (IIR) filters, this can cause ringing in the filter output, often resulting in spurious threshold crossings long after the initial stimulus artifact.

In many signal-processing applications, filter parameters are chosen to best approximate an ideal filter in the frequency domain. For this reason, higher-order Butterworth filters are often preferred for their short transition bands and flat passbands and stopbands. Here, preserving high frequency components and attenuating low frequency components of the signal are necessary to enable spike detection. However, the time domain properties of the filter are also critical because we wish to detect spikes as quickly as possible after each stimulation pulse. The post-stimulus filter ringing can be reduced or eliminated by choosing a filter with an impulse response exhibiting a fast settling time and few or no oscillations. The impulse responses to several high-pass Butterworth filter designs are shown in figure 5. Ultimately, we chose a first-order Butterworth filter design to eliminate oscillations in the impulse response. We further found that setting the high-pass cutoff frequency at higher values resulted in impulse responses with faster settling times.

Zoom In Zoom Out Reset image size

Our final spike filter implementation consisted of a first-order Butterworth high-pass filter with a 750 Hz cutoff frequency as it minimized filter ringing and exhibited a fast settling time after stimulation. The effects of this filter, and the more common fourth-order, 250 Hz Butterworth filter, on an example stimulus artifact are demonstrated in figure 6. The effects can particularly be seen in panels 6(c) and (d). Ultimately, this modified filter reduces, but does not eliminate, the distortion of the stimulus artifact. Further increases in cutoff frequency, while resulting in faster settling times, would attenuate much of the signal power in the spike band, commonly cited as approximately 300–6000 Hz [27, 28], with much of the information at about 1000 Hz [29]. While the 750 Hz cutoff frequency is within this band, most of the spike signal is preserved because of the gradual roll-off of the first-order filter, as shown in figure 7.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Stimulation waveforms are distorted by IIR spike filters, resulting in long duration artifacts. The first-order 750 Hz high-pass filter reduces, but does not eliminate, this distortion after the end of the stimulus pulse. We therefore implemented a method of blanking the recorded signal during pulse delivery to eliminate the bulk of the stimulus artifact itself, thus avoiding contaminating the filter input with a large impulse-like response. This was achieved using a digital sample-and-hold function that was implemented in the amplifier prior to digital filtering on the NSP.

While signal blanking was effective at removing stimulus artifacts, additional artifacts were caused by discontinuities at the offset of blanking, as highlighted in figure 8(a). These discontinuities, while much smaller in amplitude than the stimulus artifact, often resemble a step input. The resulting step response drives the filtered signal near or past the spike threshold, resulting in spurious artifacts after blanking. The ringing resulting from this step response can be seen in figure 8(a) between approximately 1.5 and 4 ms. Because filter ringing persists in the signal after blanking, the criteria used to tune the digital filter in response to stimulus artifacts still apply after blanking the signal during stimulation. In figure 8(b), the use of a first-order 750 Hz high-pass Butterworth filter effectively eliminates oscillations in response to the discontinuity at the offset of the blanking period and does not introduce additional artifacts after blanking. The combination of signal blanking and a digital filter with a fast settling time after impulse or step inputs can eliminate stimulus artifacts and resume spike detection shortly after each stimulus pulse.

Zoom In Zoom Out Reset image size

In our final implementation, the recorded voltage signals were blanked during each stimulation pulse using a sample-and-hold period of 1467 µs (44 samples at 30 kHz). This includes 60 µs before each stimulation pulse, 700 µs while the stimulator is delivering current, and an additional 707 µs as the voltage begins to recover from the anodal phase, which in some cases saturates the amplifier. One additional sample was blanked in software, after the application of the digital filter and threshold, to discard occasional false spikes caused by the discontinuity at the offset of the blanking period. This complete implementation allowed for reliable spike detection as soon as 740 µs after the offset of each stimulus pulse and a total signal-blanking duration of 1500 µs.

To quantify the ability of the system to detect spikes over the inter-stimulus interval, we measured the occurrences of spikes in 0.5 ms time bins between successive stimulus pulses. Spikes were detected and recorded online from both recording arrays while the participant performed an imagined reaching and grasping task. Task-irrelevant stimulation was delivered in bursts on a group of electrodes with and without artifact rejection. The histogram in figure 9(a) demonstrates that artifacts dominate the neural recordings when a fourth-order Butterworth 250 Hz high-pass filter is used without signal blanking. Most detected 'spikes' occur between 0–0.5 ms and 3.5–4.5 ms, which can be explained by the initial stimulus pulse and filter ringing, as seen in figure 6(c). When the first-order 750 Hz high-pass filter and signal blanking are implemented, spikes can be recorded throughout the interval between the blanking period and subsequent stimulus pulse (figure 9(b)). While a slightly elevated spike rate occurs from 1.5 to 5 ms after each stimulus pulse, the flat distribution from 5 to 10 ms is not significantly different from the distribution seen in the absence of stimulation (two-sample Kolmogorov–Smirnov test, p  =  0.058).

Zoom In Zoom Out Reset image size

As discussed previously, the spike filter was tuned to quickly return to steady state after perturbations related to stimulation. However, the filter parameters also have an effect on overall signal quality. To investigate the effects of cutoff frequency on signal quality, we performed an analysis based on the filter analysis in Lempka et al [23] using neural recordings in the absence of ICMS. The relationships of cutoff frequency on peak-to-peak voltage, noise level, SNR, and unsorted threshold crossings are shown in figure 10. We validated that increasing high-pass cutoff frequency decreased in peak-to-peak voltage, while increasing SNR. These results agree with the trends presented in Lempka et al. Increasing the first-order high-pass cutoff frequency from 250 Hz to 750 Hz resulted in a 20.0 µV decrease in median peak-to-peak voltage (−13.8% change, figure 10(a)), a 3.6 µV decrease in the median standard deviation of the noise (−25.6% change, figure 10(b)), and an increase in SNR of 0.67 (13.0% change, figure 10(c)). The change in cutoff frequency had a minimal effect on the median number of detected unsorted threshold crossings (−1.8% change, figure 10(d)).

Zoom In Zoom Out Reset image size

The complete artifact rejection scheme was validated online in a bidirectional BCI control task where the participant was required to transfer an object across a table as many times as possible in 2 min using a BCI-controlled robot arm. Microstimulation delivered through a group of eight electrodes was triggered by the reaction torque measured at the index finger motor as the robot hand contacted the object. Stimulation electrodes and amplitudes were chosen to generate large artifacts and were not necessarily selected to provide feedback that would lead to improvements in the task. Task performance for each condition (no stimulation, stimulation without blanking, and stimulation with blanking for both the fourth-order 250 Hz high-pass filter and first-order 750 Hz high-pass filter) is summarized in figure 11, and representative trials with and without artifact rejection can be seen in supplementary video 1. A one-way analysis of variance revealed a significant difference in number of transfers between task conditions (F(5,12)  =  75.64, p  <  0.001). Post-hoc Tukey tests revealed that the type of filter caused no significant difference in performance without stimulation (p  =  0.462). Further, when the full artifact rejection scheme was implemented, including the first-order filter, task performance with stimulation was not significantly different than both the no stimulation fourth-order filter condition (p  =  0.761) and the no stimulation first-order filter condition (p  =  0.994). When artifact rejection was not included, stimulation significantly impaired performance compared to all other conditions (p  <  0.001) to the point that the task was essentially impossible. A mean (± standard deviation) of 3766  ±  322 pulses (i.e. 31% of the trial) were delivered during each trial with artifact rejection. 3411  ±  1619 pulses (i.e. 28% of the trial) were delivered during all other stimulation trials.

Zoom In Zoom Out Reset image size

During trials without artifact rejection, the participant typically had no volitional control of the arm during ICMS. In some of these trials, the hand continued to grasp the object, resulting in stimulation pulses that continued for the maximum 15 s. After this, the electrodes were temporarily, but automatically disabled for the subsequent 15 s due to a 50% duty cycle enforced to ensure a safe ICMS protocol. During these 15 s periods, the task appeared to revert to the no stimulation condition, allowing the participant to complete a small number of transfers, as seen in the two ICMS conditions without signal blanking in figure 11. During the fourth-order filter condition without signal blanking, the artifact-corrupted signal frequently resulted in a decoded velocity towards the top-right of the workspace. While the participant reported that he had little-to-no volitional control during stimulation, the object would sometimes drop out of the hand in the correct side of the workspace, or would be held over the correct side of the workspace until the 15 s duty cycle elapsed and he could drop the object in the correct location without ICMS. This strategy permitted an average of five transfers per trial during this condition.