Improved prediction of bimanual movements by a two-staged (effector-then-trajectory) decoder with epidural ECoG in nonhuman primates

Two healthy rhesus macaque monkeys (Macaca mulatta, denoted as M23 and M24) were trained for the experiment. For the electrode implant surgery, the monkeys were prepared with sterile, anesthetic surgical procedures. A licensed veterinarian was present throughout surgery to induce anesthesia and to monitor and record all measured physiological variables. 1 h before the surgery, the animal was intramuscularly (I.M.) injected with atropine sulfate (0.08 mg/kg) to prevent excessive salivation during the surgery. One-half hour later, it was sedated with tiletamine-zolazepam (Zoletil^®, 10 mg kg⁻¹, I.M.), intubated, and placed under isoflurane anesthesia. A saline drip was maintained through an intravenous catheter placed in a leg vein. Throughout the surgery, body temperature, heart rate, blood pressure, oxygen saturation and respiratory rate were continuously monitored. The monkey was placed in a stereotaxic frame, the scalp was incised, and a craniotomy with a 2.5 cm radius was performed at both hemispheres, but the dura was left intact.

Two 32-channel platinum ECoG electrode arrays (Neuronexus, USA) were implanted in each hemisphere of each monkey covering from premotor cortex (PMC) to primary somatosensory cortex (S1). The diameter of the ECoG electrodes was 300 μm and inter-electrode distance was 3 mm. The ground electrode was placed in the midline. All electrodes were implanted in the epidural space (figure 1). The ECoG electrodes were connected to connectors (Omnetics, USA) affixed to the skull with dental cement and titanium screws (Vet implants, USA). We also implanted customized titanium head holders to fix the monkey's head in place.

Zoom In Zoom Out Reset image size

After the implant surgery, the monkeys were allowed to recover in their cages, the size of which was consistent with the Guide for the Care and Use of Laboratory Animals [31]. The temperature was maintained at 24  ±  4 °C and humidity was maintained at 50  ±  10%. The light was controlled as 12 h for day and 12 h for night. All experimental procedures were approved by the Seoul National University Hospital Animal Care and Use Committee (IACUC No. 13-0314).

The monkey sat on a chair, and its head was fixed in place with a head holder. The monkey was trained to push the buttons when lit. The button, trigger and juice reward systems were operated by MATLAB software (Mathworks, USA) and NI-PCI-6221 (National Instruments, USA), and the training and task operating system was set as shown in figure 2. The ready buttons were set in the natural position when the monkey sat on the chair, and cue buttons were set inside of the monkey's line of sight. The distance from the ready button to the cue button was approximately 30 cm, which is almost the maximum of the monkey's arm movement range. The ready buttons flashed green, and the cue buttons could flash red and blue.

Zoom In Zoom Out Reset image size

The experiment trial started when the monkey pushed the ready button. The time interval between the trial start and the cue on was randomly distributed between 3 and 7 s. Each trial is defined as one of three tasks, which were distinguished by the color of the light and the movement type. The cue button's light color indicates the target arm. For example, a blue light means, 'move the left arm', and red light means, 'move the right arm'. The cue buttons were turned on randomly so that each trial was a random task. For unimanual tasks, the target arm should move and push the cue button, while the non-target arm should push the ready button until each task was finished. For bimanual tasks, the cue buttons were simultaneously lit blue and red so that both arms should move and push the cue buttons at the same time. If the monkey pushed the cue buttons on fixed time (1 s) and if the monkey pushed the correct buttons, the trial was finished the monkey received a reward of some water. We recorded the monkey's ECoG and motion data for four months, one to four sessions per day, and each session included 200 trials.

The recording started three weeks after the electrode implant surgery. After the surgery, some of the electrodes developed bad channels because the connectors were damaged, so neural signals were recorded from 63 electrodes for M23 and 51 electrodes for M24. The signals were recorded by an EEG 1200 (Nihon Kohden, Japan) at a sampling rate of 1 kHz per channel. To record the arm trajectory, six wireless IMU (inertial measurement unit, XSENS, Netherlands) were used for motion tracking. The trackers were attached to both wrists, the upper arms, and on the back of the monkey. Between the shoulders and recorded with a 20 Hz sampling rate.

For data preprocessing, we used the MATLAB software and the EEGLAB toolbox [32]. The ECoG data was band-pass filtered from 0.3 to 200 Hz, and notch filtered 60, 120, and 180 Hz to remove power noise. The time-frequency representation, or the scalogram, was then made for decoding. This high-dimensional predictor vector, ${\rm Scal}{{{\rm o}}_{{\rm ch,freq,lag}}}$ describes the spatio-spectro-temporal information of the signals during the previous 500 ms at each electrode ch, frequency bin freq, and time lag lag. For each electrode, a 500 ms duration ECoG data was transformed into a spectrogram using wavelet transform and converted into a scalogram of 10 frequency bins (10–110 Hz, 10 Hz interval) and 5 time lags (100 ms interval). Scalograms were calculated at 50 ms intervals (90% overlap), which are matched with motion tracker recoding time (20 Hz, 50 ms). Except for the bad channels, we used 63 channels for M23 and 51 channels for M24. Thus, the number of variables in the scalogram as a predictor for M23 was 3150 (63 electrodes, 10 frequency bins, and 5 time lags), and for M24 was 2550 (51 electrodes, 10 frequency bins, and 5 time lags).

The body-centered 3D arm trajectories were converted into xyz coordinates, calculated by Euler angles of each motion tracker's obtained data. The trajectories were high pass filtered at 0.05 Hz to remove baseline drift and DC offset.

Previous research suggested that the PLS (partial least square) regression method could be used for predicting arm trajectories [10, 11]. PLS regression is useful when the number of predictors is much greater than the number of responses. Furthermore, the ${\rm Scal}{{{\rm o}}_{{\rm ch},{\rm freq},{\rm lag}}}$ has a high dimensionality and has high correlations between scalogram components, so that the PLS regression could prevent over-fit problems [33]. Therefore, using this method we found a linear relationship between predictor variables (scalogram) and response variables (prediction of arm movements).

The decoding model was calculated from the training data, where a 10-fold cross validation was performed. r values, correlation coefficient between the predicted value and the observed values, were also calculated for comparing two-staged method and the conventional method, which decodes brain signals without distinguishing between effector and motion states.

To decode the brain signals, we proposed the use of a two-staged decoder which combines the classification of movement conditions and PLS regression (figure 3). The proposed method is comprised of two steps: The first step is an effector classifier, which discriminates the movement condition, and the second step is predicting hand trajectory using PLS regression.

Zoom In Zoom Out Reset image size

For preprocessing, we prepared a training set based on each movement condition. The scalogram was gathered by four movement conditions—left unimanual movement, right unimanual movement, bimanual movements, and the stationary state from the observed trajectory and stacked every 50 ms. It would be used to create the effector classification model and the movement prediction model. All movement prediction model, including a conventional method (single stage decoder, which does not include the use of an effector classifier) and the two-staged decoder, encoded the movement of both arms, regardless of movement condition.

When we obtained new ECoG data, we prepared a new scalogram. It was then used to detect the movement condition through the use of a classification model created from the training set. A linear discriminant analysis and 500 ms of scalogram (five columns of time lags in scalogram) were used for the classification with the MATLAB software.

After the new scalogram had determined which effector as ${\rm Scal}{{{\rm o}}_{{\rm effector}\,{\rm ch},{\rm freq},{\rm lag}}}(t)$, the predicted motion ${{M}_{{\rm effector}}}(t)$ could be modeled as each of their linear combinations:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \begin{array}{@{}ccccccc@{}} & {{M}_{{\rm lm}}}(t)={{a}_{{\rm lm}0}}+\underset{{\rm ch}}{\mathop \sum }\,\underset{{\rm freq}}{\mathop \sum }\,\underset{{\rm lag}}{\mathop \sum }\,{{a}_{{\rm lm}\,{\rm ch},{\rm freq},{\rm lag}}}\cdot {\rm Scal}{{{\rm o}}_{{\rm lm}\,{\rm ch},{\rm freq},{\rm lag}}}(t)+{{\varepsilon }_{{\rm lm}}}(t); \nonumber \\ & {{M}_{{\rm rm}}}(t)={{a}_{{\rm rm}0}}+\underset{{\rm ch}}{\mathop \sum }\,\underset{{\rm freq}}{\mathop \sum }\,\underset{{\rm lag}}{\mathop \sum }\,{{a}_{{\rm rm}\,{\rm ch},{\rm freq},{\rm lag}}}\cdot {\rm Scal}{{{\rm o}}_{{\rm rm}\,{\rm ch},{\rm freq},{\rm lag}}}(t)+{{\varepsilon }_{{\rm rm}}}(t); \nonumber \\ & {{M}_{{\rm bm}}}(t)={{a}_{{\rm bm}0}}+\underset{{\rm ch}}{\mathop \sum }\,\underset{{\rm freq}}{\mathop \sum }\,\underset{{\rm lag}}{\mathop \sum }\,{{a}_{{\rm bm}\,{\rm ch},{\rm freq},{\rm lag}}}\cdot {\rm Scal}{{{\rm o}}_{{\rm bm}\,{\rm ch},{\rm freq},{\rm lag}}}(t)+{{\varepsilon }_{{\rm bm}}}(t); \nonumber \\ & {{M}_{{\rm ss}}}(t)={{a}_{{\rm ss}0}}+\underset{{\rm ch}}{\mathop \sum }\,\underset{{\rm freq}}{\mathop \sum }\,\underset{{\rm lag}}{\mathop \sum }\,{{a}_{{\rm ss}\,\,{\rm ch},{\rm freq},{\rm lag}}}\cdot {\rm Scal}{{{\rm o}}_{{\rm ss}\,\,{\rm ch},{\rm freq},{\rm lag}}}(t)+{{\varepsilon }_{{\rm ss}}}(t), \nonumber \\ \end{array}\nonumber \end{align}$$

where ${{a}_{{\rm effector}\,0}}$ is the intercept, ${{a}_{{\rm effector}\,{\rm ch},{\rm freq},{\rm lag}}}$ is the coefficient for the scalogram component at electrode $ch$, frequency bin $freq$, and time lag $lag$, and ${{\varepsilon }_{{\rm effector}}}(t)$ is the residual error. The subscripts $lm,~rm,~bm$ and $ss$ are abbreviated for left unimanual movement, right unimanual movement, bimanual movements, and stationary state, effector as movement type.

Using the coefficients of effector model $\left\{{{a}_{{\rm effector}~0}},\right.$$\left.{{a}_{{\rm effector}\,{\rm ch},{\rm freq},{\rm lag}}} \right\}$ we calculated the predicted x, y, and z arm coordinates and by repeating this process, we were able to predict both arm motion trajectories. In a test across the number of factors, performance reliably saturated at around 300 utilized factors. Thus, we empirically selected 300 PLS components as PLS parameters for all the models in our study.

To calculate the hand trajectory that could be extracted from each cortical area, each sub-band and each time lag in the ECoG signal, we quantified the spatio-spectro-temporal contribution of brain activity for predicting each arm x, y, z coordinate. Three different weight values were calculated from the coefficients of each effector model $\left\{{{a}_{{\rm effector}\,{\rm ch},{\rm freq},{\rm lag}}} \right\}.$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \begin{array}{@{}ccccccc@{}} & {{W}_{{\rm s}}}\left( {\rm ch} \right)=\frac{\mathop{\sum }_{{\rm freq}}\mathop{\sum }_{{\rm lag}}|{{a}_{{\rm effector}\,{\rm ch},{\rm freq},{\rm lag}}}|}{\mathop{\sum }_{{\rm ch}}\mathop{\sum }_{{\rm freq}}\mathop{\sum }_{{\rm lag}}|{{a}_{{\rm effector}\,{\rm ch},{\rm freq},{\rm lag}}}|}; \nonumber \\ & {{W}_{{\rm f}}}\left( {\rm freq} \right)=\frac{\mathop{\sum }_{{\rm ch}}\mathop{\sum }_{{\rm lag}}|{{a}_{{\rm effector}\,{\rm ch},{\rm freq},{\rm lag}}}|}{\mathop{\sum }_{{\rm ch}}\mathop{\sum }_{{\rm freq}}\mathop{\sum }_{{\rm lag}}|{{a}_{{\rm effector}\,{\rm ch},{\rm freq},{\rm lag}}}|}; \nonumber \\ & {{W}_{{\rm t}}}\left( {\rm lag} \right)=\frac{\mathop{\sum }_{{\rm ch}}\mathop{\sum }_{{\rm freq}}|{{a}_{{\rm effector}\,{\rm ch},{\rm freq},{\rm lag}}}|}{\mathop{\sum }_{{\rm ch}}\mathop{\sum }_{{\rm freq}}\mathop{\sum }_{{\rm lag}}|{{a}_{{\rm effector}\,{\rm ch},{\rm freq},{\rm lag}}}|} \nonumber \\ \end{array}\nonumber \end{align}$$

The spatial contribution ${{W}_{{\rm s}}}\left( {\rm ch} \right)$ of each recording electrode ch was calculated based on the ratio of the weight from the frequency bins and the time lags in the recording electrode to the total weight from all of the electrodes, frequency bins and time lags. The spectral contribution ${{W}_{{\rm f}}}\left( {\rm freq} \right)$ of each frequency band freq and the temporal contribution ${{W}_{{\rm t}}}\left( {\rm lag} \right)$ of each time bin lag were also calculated in the same manner as the spatial contribution. Thus, the spatial contribution ${{W}_{{\rm s}}}\left( {\rm ch} \right)$ was used to quantify the contribution of each channel for making predictions across all of the frequency and time bins. Spectral contribution ${{W}_{{\rm f}}}\left( {\rm freq} \right)$ and temporal contribution ${{W}_{{\rm t}}}\left( {\rm lag} \right)$ were used to quantify the contributions of each frequency band and time lag, respectively. These contributions can be interpreted as how each scalogram component contributes to the decoding performance.

To evaluate the performance of the two-staged decoder, we calculated the correlation coefficient for 32 sessions in M23, and 39 sessions in M24. The two-staged decoder showed a meaningful performance in decoding the brain signals for arm movement. It was more accurate than the conventional method (single stage decoder), which does not include the use of an effector classifier. The conventional method decodes the brain signal without distinguishing between the effector and motion states. Thus, the decoder used in the conventional method did not include different weight (or regression coefficients) for each effector, and it was evaluated with a fixed set of prediction weights for any hand or any motion state. In contrast, the two-staged method proposed herein considers the effector. As shown in figure 4, the correlation coefficient between the predicted arm's coordinates and the observed trajectory was significantly improved for all calculated axis, except for the right arm z-coordinate for M23 and the left arm x-coordinate for M24 (*: Bonferroni-corrected p  =  0.01, two tailed paired t-test). The correlation coefficients (r) are shown in table 1 for the left arm x-, y-, z-coordinate and the right arm x-, y-, z-coordinate for M23 and M24 respectively (mean $\pm$ SD, n  =  32 sessions in M23 and 39 sessions in M24). Due to the narrow range of movement along the x-axis in M24, the correlation coefficients for the y and z coordinates were higher than that of the x coordinate. The other axes were sufficiently good to decode the arm movement. Thus, the two-staged decoder can be used to predict arm movements in both unimanual and bimanual tasks.

Table 1. The correlation coefficient for predicted bimanual movements' coordinates using the two-staged decoder.

Subject	Left arm (correlation coefficients r)			Right arm (correlation coefficients r)
	x-axis	y-axis	z-axis	x-axis	y-axis	z-axis
M23	0.48  ±  0.10	0.64 $\pm$ 0.11	0.54  ±  0.11	0.54  ±  0.16	0.66  ±  0.09	0.52  ±  0.12
M24	0.24  ±  0.09	0.59 $\pm$ 0.08	0.64  ±  0.08	0.37  ±  0.08	0.52  ±  0.08	0.60  ±  0.08

Zoom In Zoom Out Reset image size

In the evaluation of the performance of the effector classifier in the two-staged decoder, a mean accuracy of over 80% was found in the classification between the movement and the stationary state. In addition, the accuracy for the four classes (left unimanual movement, right unimanual movement, bimanual movements and stationary state) was about 70%, as shown in table 2. Furthermore, we also calculated F1 score, which is the evaluation metrics for classification models. It takes both precision (positive predictive value) and recall (true positive rate) into consideration. F1 score reaches its best value at 1 and worst at 0 in the classification performance. Therefore, based on the classification accuracy and F1 score, the origin of some of the errors in the two-staged decoder appeared to be in the performance of the classifier.

To assess which of the electrodes, frequency bands, and time points (lags) play an important role in predicting arm movement in the two-staged decoder, a contribution analysis was conducted according to the status of the four effectors (left unimanual, right unimanual, bimanual movements and stationary state). In the electrode contribution analysis, the channel contribution of the movement state in M23 is generally involved in various areas of the brain's electrodes, compared to the stationary state's contribution which is involved in just the local areas such as SMA or M1, as shown in figure 5(a) (p  <  0.01, Wilcoxon signed-rank test). It is noteworthy that the right supplementary motor cortex showed a dominant role in arm movement states. In addition, the left and right posterior parietal cortices in the bimanual movements condition appeared to be more involved in the prediction compared to the unimanual movement conditions. Unfortunately, there were too many bad channels in M24, so we could not make a meaningful interpretation. In the frequency band contribution analysis, the alpha band and the beta band showed higher contributions in all conditions compared to the gamma band frequency (p  <  0.01, Wilcoxon signed-rank test, asterisks in figure 5(b)). Interestingly, in the time lags contribution analysis, all time points generally showed even contributions (p  <  0.01, Wilcoxon signed-rank test, asterisks in figure 5(c)).

Zoom In Zoom Out Reset image size

To validate the long-term durability of the signal quality, we calculated correlation coefficients r for each of the four months' sessions. (Same-session prediction) We acquired a high r as predictive accuracies that showed no significant monotonic decrease over the four month period. (Figure 6(a), p  >  0.01, $\left| \rho \right|$  <  0.45, Spearman correlation test). We measured the power spectrum density and in vivo impedance of whole channels, and, except for the bad channels, no differences were found over time (<40 kΩ at 1 kHz). The findings indicated that the signal quality and the epidural ECoG recording system showed long-term durability.

Zoom In Zoom Out Reset image size

Next, to demonstrate the long-term stability of the two-staged model, we verified that the decoding model would have applied other data recorded later (cross-session prediction) and used a different validation method in the cross-session. The decoding model constructed from training data from the first session was used to predict all validation data in the rest of the sessions for four months. We calculated the correlation coefficient difference ($\Delta r$) between the first session and the remaining sessions. Figure 6(b) shows the mean correlation difference for six coordinates for both arms (left and right x, y, z positions) and no significant monotonic decrease over time was found (p  >  0.01, $\left| \rho \right|$  <  0.3, Spearman correlation test). This demonstrates the long-term stability of the new method for BCIs.