A Sliced Inverse Regression (SIR) Decoding the Forelimb Movement from Neuronal Spikes in the Rat Motor Cortex

Several neural decoding algorithms have successfully converted brain signals into commands to control a computer cursor and prosthetic devices. A majority of decoding methods, such as population vector algorithms (PVA), optimal linear estimators (OLE), and neural networks (NN), are effective in predicting movement kinematics, including movement direction, speed and trajectory but usually require a large number of neurons to achieve desirable performance. This study proposed a novel decoding algorithm even with signals obtained from a smaller numbers of neurons. We adopted sliced inverse regression (SIR) to predict forelimb movement from single-unit activities recorded in the rat primary motor (M1) cortex in a water-reward lever-pressing task. SIR performed weighted principal component analysis (PCA) to achieve effective dimension reduction for nonlinear regression. To demonstrate the decoding performance, SIR was compared to PVA, OLE, and NN. Furthermore, PCA and sequential feature selection (SFS) which are popular feature selection techniques were implemented for comparison of feature selection effectiveness. Among SIR, PVA, OLE, PCA, SFS, and NN decoding methods, the trajectories predicted by SIR (with a root mean square error, RMSE, of 8.47 ± 1.32 mm) was closer to the actual trajectories compared with those predicted by PVA (30.41 ± 11.73 mm), OLE (20.17 ± 6.43 mm), PCA (19.13 ± 0.75 mm), SFS (22.75 ± 2.01 mm), and NN (16.75 ± 2.02 mm). The superiority of SIR was most obvious when the sample size of neurons was small. We concluded that SIR sorted the input data to obtain the effective transform matrices for movement prediction, making it a robust decoding method for conditions with sparse neuronal information.

Introduction

In order to improve daily life activities for paralyzed patients, the establishment of a non-muscular communication interface between brain neurons and machines has rapidly developed over the last two decades (Schwartz, 1993, 1994; Donoghue, 2002; Schwartz et al., 2006; Velliste et al., 2008). With assistance from stable generated brain-derived control signals incorporated with prosthetic devices and motor functions, paralyzed patients now possibly regain their ability to move a computer cursor (Kennedy et al., 2000; Hochberg et al., 2006; Gilja et al., 2015), control an anthropomorphic prosthetic arm (Wodlinger et al., 2015), and drive a prosthetic device (Hochberg et al., 2012; Collinger et al., 2013) through a brain-machine interface (BMI). One important challenge to BMI is how to design an appropriate neural decoder (Pohlmeyer et al., 2014). To address the challenge, previous studies have carefully utilized training paradigms that have been designed for a BMI decoder and controller. For brain-derived control signals, neural decoding is an indispensable technique that translates neuronal activities to physical states, such as the position of a foraging rat (Brown et al., 1998), arm movement (Ashe and Georgopoulos, 1994), movement speed (Moran and Schwartz, 1999), hand position (Paninski et al., 2004), and joint angular velocity (Reina et al., 2001).

A population vector algorithm (PVA), one method for decoding motor cortical activity, assumed that a neuron's firing rate is related to the velocity vector of movement. PVA categorizes each neuron's contribution into directional and distance information of the movement by a directional tuning function under uniform variance conditions (Georgopoulos et al., 1988). A previous study showed that PVA decoding could expose the visuomotor coordinate transformations between visual and motor information by processing masses of neuronal activities recorded from relative brain regions (Takeda and Funahashi, 2004; Watanabe et al., 2009). PVA presented superior performance in predicting hand path throughout reaching tasks (Schwartz, 1994). However, a uniformity constraint is usually not the case for real experiments, and the equality of the tuning function is variable because of the small amount of unit recordings in realistic applications (Schwartz et al., 2001). To compensate for the non-uniform preferred directions in the population of recorded neurons, an optimal linear estimator (OLE) was proposed to define the preferred direction of each neuron using the center of mass of the tuning function (Salinas and Abbott, 1994). Requiring large numbers of neurons with a temporal solution of 10–100 ms, PVA and OLE studies successfully predicted the kinematic parameters of a primate arm movement (Schwartz et al., 2001; Takeda and Funahashi, 2004; Watanabe et al., 2009).

A Bayesian decoder, a probabilistic decoding technique, could achieve accurate offline trajectory reconstructions by combing simple trajectory models (Yu et al., 2007). However, off-line reconstruction may not be suitable for online prosthesis control because the essential features of a real prosthesis are not acquired, and the system dynamics may vary because the user is in a closed loop. Furthermore, offline and online approaches resulted in different parameter choices for decoding algorithms (Cunningham et al., 2011). Therefore, the neural decoders and the motor prosthesis must be tested online even though online control experiments are more expensive both in terms of physical resources and time (Gilja et al., 2011). A recursive Bayesian decoder, i.e., a Kalman filter, was developed to decode the neural data recorded in the monkey motor and premotor cortex in response to goal-directed reaching movements (Shenoy and Carmena, 2014). It yielded high decoding performance and accurate trajectory prediction when the probability modeling assumptions were satisfied. For online purposes, a modified Kalman filter that transforms the acquired neural signals into a controller input was further designed for online cursor-control tasks and resulted in high performance in rhesus monkeys (Gilja et al., 2012). To adapt decoders to the dynamics of a prosthetic device and its environment, a likelihood gradient ascent and a self-recalibrating classifier were proposed to update decoder parameters during closed-loop BMI operation and normal use (Dangi et al., 2013; Bishop et al., 2014). Additionally, neural networks developed from probabilistic aspects were designed in an online setting (Sussillo et al., 2012) and in a real-time setting (Dethier et al., 2013).

A selection of cortical neurons, instead of all available neurons, used in the encoding process could improve the control performance of the neuroprosthetic system, such as robotic arms (Wahnoun et al., 2006). However, neural coding mechanisms evolve with time, individual experience, and the learning process (Nicolelis, 2001), i.e., the contribution from individual neurons may vary considerably from day to day (Carmena et al., 2003). It has been observed that neuronal activity for monkey was not as stable from day to day (Sadtler et al., 2014). The decoding algorithms which require previous day's observation of neural activity, such as PVA and OLE, may be affected by neuron's stability. For this reason, the selection of cortical neurons becomes an essential issue for the decoding processes, especially after continuous practice and learning. Furthermore, long-term inflammatory responses lead to a gradual decrease of recording quality and the eventual breakdown of the electrode's recording ability (Polikov et al., 2005; Schwartz et al., 2006). Losing neural signals over time will result in chronic coding failure. Thus, a decoder that has the capability to process recorded signals from a small number of neurons will become more important.

Sliced inverse regression (SIR) is a data-analytic tool that can effectively perform nonlinear regression based on a small number of inputs (Li, 1991). It divides the range of output variable into several intervals and partitions the input data into several slices according to the corresponding output value. Each slice consists of data with a similar contribution to output estimation. SIR then applies a weighted principal component analysis (PCA) to theses slice means of data to find effective dimension reduction directions for general and flexible setups. Each slice gains weight according to its contribution to output estimation. With a simple inverse regression model, SIR requires a low computational cost and retains reliably stable subspaces to extract primary information from noisy data with effective dimension reduction directions. Because of the good performance in dimension reduction and data de-noising, SIR has been widely applied to data-intensive marketing environments (Naik et al., 2000), data classification (Dai et al., 2006; Wu, 2008), and medical images (Wu and Lu, 2007; Lu, 2008; Tu et al., 2015).

It has been known that the number of recorded neurons in rodents is less than primates in BMI applications. Owing to the abilities of processing small number of input variables and slicing data for inverse regression model, SIR was considered as a neural decoding algorithm to provide realistic behavior interpretation for brain-derived control tasks in this study. A surrogate driving task with lever-pressing was designed for evaluating the efficacy of the proposed decoding algorithm in predicting motor functions of rats. The signals of lever-pressing and the spike trains related to the intended forelimb movement trajectories were simultaneously acquired during the task. This study adopts SIR to predict intended forelimb movement trajectories according to the recorded neurons from primary motor (M1) cortex. We presented experimental validation of the proposed decoding system using the recorded neurons to predict forelimb movement. We demonstrated that the proposed SIR decoding algorithm can not only extract primary features from the small number of neurons but perform more accurate prediction of intended forelimb movement trajectories than PVA and OLE which usually require hundreds of neurons in primates.

Materials and Methods

Animal Preparation

Four male Wistar rats weighing 250–300 g (BioLASCO Taiwan Co., Ltd., Taiwan) were used in this study. All rats were individually housed in a 12 h light-dark cycle (light from 7.00 to 19.00 h) at room temperature (22 ± 1°C) with access to food and water ad libitum in the experimental animal center of National Yang Ming University. All experiments were conducted according to standards established in the Guide for the Care and Use of Laboratory Animals, which has been approved by the Institutional Animal Care and Use Committee at the National Yang Ming University.

Animal Training and Behavioral Tasks

The rats were trained to use their right forelimb to press a lever to obtain the water reward for a week before electrode implantation, as shown in Figure 1. The rat was placed in the lab-designed Plexiglas testing box with a 15-cm tall lever above the floor on the left side (Figure 1A) and a water-feeder with a flow rate of 1 ml/time on the right side as shown in Figure 1B (Lin et al., 2016). Figure 1C shows the experimental setup where a rat was pressing a lever for the water reward. Before achieving the successful lever-pressing training, all rats underwent water deprivation for 8 h per day. In this study, we have defined the criterion for the successful training was to consecutively repeat the lever-pressing and water-drinking for five times during daily 5-h sessions (9:00–14:00), for 4 days at the most. Once the rats learned the association, they always kept the skilled concept (Lin et al., 2016).

FIGURE 1

Figure 1. Experimental setup and protocol. A perspective drawing (A) and vertical view (B) of the Plexiglas testing box. A lever is on the left side of the barrier and a water supply is on the right side. (C) A rat is using his right forelimb to press the lever to obtain a water reward. Simultaneously, his forelimb movement trajectory is videotaped by a camcorder approximately 25 cm away from the box, and the neuronal activities are recorded by the implanted electrode.

Surgery and Electrophysiological Mapping

Animals were anesthetized with pentobarbital (50 mg/Kg, i.p.) and were placed on a stereotaxic apparatus (Model 900, Kopf Instruments, Tujunga, CA, USA). A 16-channel stainless microwire electrode array (diameter of 0.002 ft., California Fine Wire Co., Ltd, Grover Beach, CA, USA) was inserted vertically and was implanted into layer V of the M1 cortex (2–4 mm anterior and 2–4 mm left-lateral to bregma; 1.7 mm ventral to the cortical surface) by referring to a previous work . Here, a standard intracortical microstimulation (ICMS) technique was conducted to deduce maps of rat forelimb movement representations in the M1 cortex, which could help assess the functional integrity of M1 cortex and activate pyramidal cell fibers. ICMS delivered a 40 ms stimulus train with 0.2 ms square-wave monophasic cathodal pulses at 350 Hz to the electrodes (impedance: 200–400 kΩ at 1 kHz) by an isolated pulse stimulator at a rate of 1/s (Model 2100, A-M Systems Inc., Sequim, WA, USA). Because the intensity of the stimulating current depends on the distance between the neuron and the stimulating electrode, the threshold current intensity can be estimated by a strength-distance relationship as follows:

$\begin{array}{cc}\begin{array}{l}I=k{r}^2+I_m\hfill \end{array}& \left(1\right)\end{array}$

where k = 1292μA/mm² for direct activation, r is the distance, and I_m = 1μA. An implantation location of the electrode site was defined as valid when the rat forelimb was activated by ICMS with a current intensity less than 60μA. Then, a stainless steel screw was secured to the skull over the cerebellum as a reference electrode. Finally, the microwire electrode array was secured in the skull using dental acrylic (Type 1 Class 1, Hygenic Corp., Akron, OH, USA) and was covered with a small amount of 2% agar. For a better recovery, all rats were given an analgesia (Buprenorphine/Buprenex, 0.05 mg/kg s.c.; Reckitt Benckiset Healthcare Ltd, Hull, UK) every 8–12 h for 3 days and antibiotic treatment (Ampicillin, 100 mg/kg s.c. twice daily; Bristol Myers Squibb, New York, NY, USA) for 7 days after surgery. Following a 1-week post-surgery recovery period, all implanted rats received the behavioral task to use their forelimb to press a lever for water reward. The forelimb movements in the rat were captured by a charge coupled device (CCD) camera (DFK 21F04, Imaging Source, Bremen, Germany) and the neuronal signals were recorded by a Multi-channel Acquisition Processor (MAP, Plexon Inc., Dallas, TX, USA) through a 16-channel stainless microwire electrode array implanted in the rat M1 cortex. The detailed data recordings for forelimb movement and neuronal signals are described in the Supplementary Note 1.

Trajectory Prediction Model

This study assumed that lever-pressing forelimb movement, which was considered to be a stereotyped movement, was performed at a nearly constant distance from the CCD in each trial, i.e., the distance did not vary dramatically. The recorded trajectory might consist of major forelimb movement and minor whole body movement which led to the coupling mechanism of two-dimensional forelimb movement (see Supplementary Note 1). The coupling mechanism resulted in a nonlinear relationship between neural activity and forelimb movement, and thus caused general linear regression to fail at forelimb movement prediction. SIR performs as a nonlinear regression since it can recover the most severe nonlinearity of the data by estimating effective dimensional reduction (e.d.r.) space (Li, 1991). Therefore, SIR was adopted to predict the forelimb movement according to the neural activity in this study.

The two-dimensional trajectory movement vectors v_x and v_y in Cartesian coordinates are transformed into polar coordinates as follows:

$\begin{array}{cc}\begin{array}{c}{v}_r=\sqrt{v_x^2+v_y^2}\end{array}& \left(2\right)\end{array}$

$\begin{array}{cc}\begin{array}{c}{v}_{\theta }={\tan }^{-1}\frac{v_y}{v_x}\end{array}& \left(3\right)\end{array}$

where v_r and v_θ are the magnitude and direction, respectively.

The movement response g is assumed to be given by a deterministic function f with additive noise ε, so that

$\begin{array}{cc}\begin{array}{l}g=f\left({\beta }_{1}z,\dots ,{\beta }_{K}z,\epsilon \right)\hfill \end{array}& \left(4\right)\end{array}$

where g is v_r or v_θ, and z is firing rate data in R^p. Here, β's are unknown linearly independent projection row vectors, and K is the sufficient number of β's. The p-dimensional variable z is projected onto the K-dimensional space by functional relationship f, where p ≥ K. The combinations of β's are called the e.d.r. direction and the linear space produced by the β's are called e.d.r. space. The present study assumes the movement response g was predictable from K projected variables. To train the functional relationship f, a set of training data consisting of N training samples was prepared. According to the model assumptions, the centered inverse regression curve E(z|g) − E(z) is included in the linear subspace, which is spanned by β_kΣ_zz(k = 1, …, K), where Σ_zz represents the covariance matrix of z. SIR sorts and divides the whole data z into H intervals (slices) according to the g value. Each slice has almost equally number of observations. SIR then performs an eigenvalue decomposition of the weighted covariance matrix Σ_E(z|g) with respect to Σ_zz. The weighted covariance matrix Σ_E(z|g) is constructed as:

$\begin{array}{cc}{\Sigma }_{E\left(z|g\text{}\right)}=\frac{{\displaystyle \sum _{h=1}^H{m}_h\text{}\left({\overline{z}}_h-\overline{z}\right){\left({\overline{z}}_h-\overline{z}\right)}^{\prime }}}{\left(N-1\right)}& \left(5\right)\end{array}$

where m_h denotes the size of each slice, $\overline{z}$ is the sample mean of z, and ${\overline{z}}_h$ is the sample mean of the hth slice. The e.d.r. directions could be estimated by solving the generalized eigen-problem:

$\begin{array}{cc}{\Sigma }_{E\left(z|g\text{}\right)}{\beta }_j={\lambda }_j{\Sigma }_{zz}{\beta }_j& \left(6\right)\end{array}$

where j = 1, …, p and λ₁ ≥ λ₂ ≥ ⋯ ≥ λ_p. Then, z was further projected onto the e.d.r. space by the first K e.d.r. directions as follows:

$\begin{array}{cc}w=\left[{\beta }_{1}z,\dots ,{\beta }_{K}z\right]& \left(7\right)\end{array}$

Then, a linear combination of w was performed to predict the forelimb movement. Although a linear combination approach was adopted, SIR was considered as a nonlinear regression since there is no linearity constraint on the prediction rules. Note that the user-specified parameters of SIR are only the number of slices H and the number of components K. It has been known that SIR can provide root n consistent estimates regardless of the choice of H. A previous study has demonstrated that the performance of SIR is less sensitive to the selection of H when H was set 5, 10, and 20 (Li, 1991). Furthermore, it has been found that the first component (K = 1) is close to the e.d.r. space. Therefore, H and K were set to 10 and one, respectively, in this study.

Time-Lags and Temporal Orders

In fact, the physical relationship between the neuronal signal and the forelimb movement may imply time-lags in the neuronal signal. Previous work indicates that a model that assumes that all cells exhibit the same time-lags is computationally simple (Wu et al., 2004). Then, the optimal time-lags could be found with an empirical setting for further improvement of the decoding task. A number of time-lags (0–5 time bins, at levels corresponding from 33 to 165 ms) were evaluated for trajectory prediction by SIR.

In addition to the time-lags, the temporal order of the input is another interesting factor for the decoding issue. The information at the n^th time bin may have a relationship with that at the (n − 1)^th time bin. Hence, both current and previous neuronal activities are important and are considered as the inputs for prediction. Therefore, a tapped delay line model of neuronal activities is adopted in this study where a third-order model would consider the n^th, (n − 1)^th, and (n − 2)^th time bins as the input.

Performance Evaluation and Statistical Analysis

This study computed the root mean square error (RMSE) between true and predicted forelimb movements from movement start to endpoint in order to examine the performance of proposed decoding algorithm (Srinivasan and da Silva, 2011). The experimental trails were randomly split 70/30% into training and testing sets for each rat. Therefore, the performance of the proposed decoding algorithm on the testing set could be evaluated on each rat individually. A 10-fold cross validation was applied to avoid capitalization on chance (Efron and Tibshirani, 1994).

For statistical analysis, the predicted performance (RMSE) from 4 testing sets (145 trails; rat 74: 24 trials, rat 102: 18 trails, rat 106: 78 trails and rat 129: 25 trails) were represented as the mean ± standard error of mean (SEM). Two-way ANOVA was calculated using effects of time-lag (bin number 0, 1, 2, 3, 4, and 5) and temporal order (1-oder, 2-order, and 3-oder of time bins) as the independent variables in order to determine if there were any differences in the decoding ability based on which parameters was employed. Post-hoc comparisons were conducted using a Tukey HSD post-hoc test and the significance level was corrected to *P < 0.002 using a Bonferroni correction for the comparison of six time-lags and three temporal order. MATLAB (MathWorks, Natick, MA, USA) was used for all statistical analyses.

Results

To evaluate the decoding performance of the proposed algorithm for trajectory prediction, the majority of decoding methods, PVA, OLE, and NN, were implemented for comparison. Furthermore, two popular feature selection techniques, PCA (Wold et al., 1987) and sequential feature selection (SFS) (Aha and Bankert, 1996), were implemented for comparison of feature selection effectiveness. Although PCA was developed as dimension reduction of feature space, it could be considered as a way to select features from principle components. Then a linear regression approach was adopted to perform regression whose inputs were the features provided by PCA and outputs were forelimb movements. The number of principal components was selected according to the variance of the reconstruction error (Valle et al., 1999). The same regression procedure was applied to SFS. A set of time-lags (0–5 time-lags) was carried out to observe the effect of different delays between the neuronal activity and the forelimb movement in each method. Furthermore, a set of experiments was conducted to evaluate the effect of various temporal orders (1–3 temporal orders) used in each decoding method. The experimental data were recorded from four rats where the number of trials and number of neurons per trial for each rat are shown in Table 1. The average number of successful trials was 45.9 ± 4.3, and the average number of recorded neurons was 21.2 ± 2.6 units. The period of each trial was 0.4–0.7 s before lever-pressing and 0.2–0.4 s after the lever-pressing.

TABLE 1

Table 1. Experimental data characteristics.

Neuronal Signal Pattern during Behavior Task

The task timeline in Figure 2 presents the sequential images of the lever-pressing (Figure 2A) and the corresponding spike trains (Figure 2B). The spike trains were acquired from five neurons related to right forelimb movement and were represented as the neuronal activity histogram with 33 ms time bins. The neuronal activities distinctly increased approximately 0.4 s before lever-pressing (Figure 2B), corresponding to the second image at 01:56.332 s (Figure 2A). The maximum value in the histogram appears approximately 0.1 s before the lever-pressing. The neuronal activity has a substantial reduction after the rat completes the lever-pressing and then it re-strengthens gradually because of the redundant movement off of the lever.

FIGURE 2

Figure 2. One example of forelimb movement over time. (A) The movement video of forelimb movement. (B) The neuronal activities were recorded from five neurons during one movement displayed as spike trains and the neuronal activity histogram (a bin size of 33 ms). The red line indicated the moment when the rat presses down the lever with the right-forelimb. Moreover, our results showed that neuronal firing rates highly correlated with forelimb movement; >71% (41/57) neurons exhibited specific firing changes during movement used to discriminate directional pairs.

Effects of Different Time-Lags and Temporal Orders

Results of a particular test trial (Rat 106, 29 units, and 52th trial) were shown in Figure 3 where the trajectories were reconstructed by six decoding methods with delay activity of four time-lags. The actual trajectory (blue solid line) was compared with the decoded trajectories by SIR (black dashed line), OLE (red dashed line), PVA (green dashed line), PCA (magenta dashed line), SFS (cyan dashed line), and NN (yellow dashed line) in Figure 3. Furthermore, the results of one and three temporal orders were conducted to demonstrate the advantage of SIR with the requisite amount of input data as shown in Figures 3A,B, respectively. In the one temporal order experiment, the trajectories reconstructed by PVA and OLE obviously deviated from the actual trajectory more than that reconstructed by SIR. PCA and SFS, which perform feature selection, could achieve accurate prediction in the previous time steps but did not predict the latter time steps well. Similarly, a NN, which has learning ability for nonlinear regression, could predict the trajectory well in the first few time steps, but it did not have robust prediction performance because of the random initialization of weights that leads to prediction error. As the temporal order increased to three, all methods had more accurate prediction compared to those of one temporal order. Overall, SIR shows the best performance among the other methods, especially when the decoding methods used less neuronal activity information.

FIGURE 3

Figure 3. Reconstructed trajectories of the test trial with the use of delayed activity with four time-lags (132 ms). The actual trajectory (blue solid line) and the trajectories predicted by SIR (black dashed line), OLE (red dashed line), PVA (green dashed line), PCA (magenta dashed line), SFS (cyan dashed line), and NN (yellow dashed line) are shown for an example trial using (A) one- and (B) three- temporal orders. The trajectory reconstructed by SIR is more accurate than the other methods.

Figure 4 presents the effects of various time-lags and temporal orders in each method. Figure 4A shows the results of SIR where the smallest RMSE (8.47 ± 1.32 mm) was obtained by using four time-lags (132 ms) and one temporal order. It can be seen that SIR with four time-lags could achieve a significantly smaller RMSE than those with various time-lags [F_{(5, 54)} = 4.22, *P < 0.002 with Two-way ANOVA with Bonferroni correction, N = 145]. Furthermore, the RMSE of SIR had no conspicuous variations among the three different temporal orders. As shown in Figure 4B, OLE achieved the smallest RMSE (17.22 ± 3.80 mm) by using four time-lags and three temporal orders. However, there was no significant enhancement of the prediction performance using OLE decoding with different time-lags and temporal orders. Figure 4C shows the results of PVA where the RMSE decreased as the number of temporal orders increases. PVA resulted in an average RMSE of 21.76 ± 8.11 mm when using one time-lag and three temporal orders. Figures 4D,E showed the results of PCA and SFS, respectively, where the features were selected via these two algorithms. PCA achieved a decreasing RMSE as the number of time-lags increased and obtained the smallest RMSE (19.13 ± 0.75 mm) when using four time-lags and three temporal orders. The results of SFS did not present a decreasing RMSE as the number of time-lags increased. SFS achieved the smallest RMSE (22.75 ± 2.01 mm) when using five time-lags and three temporal orders. Figure 4F shows the results of NN where the smallest RMSEs (16.75 ± 2.02 mm) were achieved by using three time-lags and two temporal orders. The RMSEs of NN did not present a regular trend as the number of time-lags increased. The forelimb movement predictions using OLE, PVA, PCA, SFS, and NN were not affected by either time-lag or temporal order. No significant interaction of time-lag and temporal order was found in the decoding methods of OLE, PVA, PCA, SFS, and NN in comparison to SIR. These results indicated that SIR outperformed other methods for trajectory prediction.

FIGURE 4

Figure 4. RMSEs of (A) SIR, (B) OLE, (C) PVA, (D) PCA, (E) SFS, and (F) NN decoding methods plotted with various time-lags (33 ms/lag) and temporal orders. The error bars denote standard error of the mean (Mean ± SEM). The results showed that SIR is superior to the other methods for trajectory reconstruction. SIR is unaffected by temporal orders, and the best performance was achieved with four time-lags (132 ms). The symbol * indicates significant different means with P < 0.002 and analyzed by Bonferroni correction for multiple comparisons, N = 145. Mean ± SEM%.

Discussion

The main finding of this study is that a rat's forelimb movement could be successfully predicted and reconstructed using relatively few motor cortical neurons. In comparison with competing neural decoding algorithms including PVA, OLE, PCA, SFS, and NN, SIR presented an extremely superior RMSE in distance deviation between the reconstructed and real forelimb movement trajectories.

Previous studies indicated that neuronal activity discharged before the onset of the desired movement, such as the motor preparation period, and encoded behaviors (Chapin et al., 1999; Churchland et al., 2006). The kinematic parameters therefore were decoded and reconstructed with high accuracy using the neuronal activity before the onset of the movement. Hence, in this study, SIR and the competing algorithms decoded the neuronal activities during the motor preparation period for reconstruction of introduced upcoming lever-pressing. The results showed that SIR, OLE, and PCA achieved optimal efficiency when using the neuronal activities that led to the onset of forelimb movement for 132 ms. PVA, SFS and NN each achieved optimal efficiency by using the neuronal activities that led to the onset of forelimb movement for 33, 165, and 99 ms, respectively. In Figure 2B, the peak of the spike train occurred at four time-lags, prior to the onset of lever-pressing. The neuronal firing rate then declined for 0.2 s as a result of the completion of motor command transmission. Because the rat performed an unexpected forelimb swing, the neuronal firing rate increased again approximately 0.2 s after to the onset of lever-pressing.

The performance of cortical neural decoding hinged on the exploited information in chronically-recorded neuronal activities. Previous studies using PVA (Schwartz, 1994) and OLE (Salinas and Abbott, 1994) show that cortical neurons with known motor associations were chronically sampled and as many as possible were recorded. Because of the lack of a precise technique to target the modulated neurons that largely contributed to goal-directed behavior, PVA and OLE summed the weighted vectors across all neurons, performing a neuronal vote, to predict the kinematic parameters (Salinas and Abbott, 1994; Schwartz, 1994). A large number of electrodes and sample neurons (usually up to hundreds) was required for reconstruction of kinematic parameters with a high degree of accuracy (Chapin et al., 1999; Wessberg et al., 2000; Serruya et al., 2002; Taylor et al., 2002). However, the neuronal activity was not as stable from day to day (Sadtler et al., 2014). PVA and OLE may be affected by neuron's stability since they extract movement information from the selected cortical population.

In this study, we recorded only tens of neurons from rat M1 where the amount of recorded neurons was insufficient for PVA and OLE, which usually require hundreds of neurons to provide a robust neural decoding process (Takeda and Funahashi, 2004; Wahnoun et al., 2006). Compared to PVA and OLE, SIR can effectively achieve nonlinear regression from a small number of inputs (Li, 1991). SIR adopted a sliced regression framework with a sorting procedure to divide the neuronal dataset into several slices according to the sorted output variable value. Each slice contained neurons with a similar contribution to the introduced lever-pressing and was then modified by a weight. Slices containing neurons with tiny or even a null contribution to the lever-pressing may gain zero weight and can be removed from the decoding model. Multiplied by a proper weight according to the weighted PCA, a slice containing a few neurons with a high contribution presented a comparable influence on the prediction and reconstruction of the introduced lever-pressing to the neuronal vote from hundreds of neurons. Hence, SIR is able to perform forelimb prediction through a small number of neurons. On the other hand, dimensionality reduction technique factor analysis is usually adopted to describe population activity using low-dimensional set of factors and highlight feature of interest in data from a large number of recorded neurons (Sadtler et al., 2014). Although SIR could perform dimensionality reduction through weighted PCA, it preserved all recorded neurons and assigned weights to the slices according to their contribution. It learned forelimb movement prediction from whole neuronal activities regardless of neuron's stability across days. Thus, SIR was robust to uncertain variation of movements and neuronal activities across days due to the success of inverse regression and effective dimension reduction. Furthermore, using this SIR, the size of the neural decoding model topology was significantly reduced, burdened with data storage and reduced computational loading, indicating that efficiency in neural decoding in comparison to PVA and OLE was attainable. Compared to PCA and SFS, which perform feature selection and dimensional reduction, PCA could further project the data onto another space, which could lead to a better reconstruction than SFS. However, SIR outperformed PCA and SFS because SIR clusters data into each slice according to the output values. NN performed better prediction than PCA and SFS because of its nonlinearity and learning ability. Nevertheless, NN did not result in a robust reconstruction because of the mechanism of random initialization and the existence of many local optima. These comparisons indicate that the neural decoding based on SIR with one temporal order presents a smaller RMSE in reconstructing limb movement than those based on PVA or OLE with three temporal orders and those based on feature selection and learning ability. It indicates that SIR can be a more suitable solution than the commonly used linear progression methods using the neuronal ensemble inputs to predict and reconstruct the introduced limb movement.

Conclusions

Neural decoding models that require hundreds of input variables, such as PVA and OLE, not only require considerable computation but also have detrimental effects in the decoding process because of errors in assigning neuronal spikes or non-stationary noise, especially for non-adaptive models. Reducing the neurons that may cause model over-fitting emerges as a significant neural decoding issue. However, with the help of the proposed approached based on SIR, researchers can predict and reconstruct the limb movement of interest with high accuracy using only tens of neurons in a single setting. Furthermore, SIR outperformed other feature selection methods, such as PCA and SFS because of its clustering ability. SIR further achieved more robust performance than NN because there is no random initialization and local optimization problems in SIR. By indexing the contribution of multiple cortical areas with different sizes, it has become feasible to ascertain the importance of selected areas for the motor commands. This will provide valuable insights for follow-up studies in the future.

Ethics Statement

The animal use protocol listed below has been reviewed and approved by the Institutional Animal Care and Use Committee (IACUC) at the National Yang Ming University. Protocol Title: Peripheral Nerve Prostheses: A Paradigm Shift in Restoring Dexterous Limb Function. IACUC Approval No: 1050622. Period of Protocol: Valid From: 08/01/2016. To: 07/31/2020. There is no endangered animal species involved in this study.

Author Contributions

SY, HYL, and YC designed the project, organized the entire research. SL, LL conceived the experiments. CW, PC, TP, HCL, and WH conducted the experiments. HC, HHSL, and YL analyzed the results. SY, HYL, and YC wrote the manuscript. All authors discussed the results and reviewed on the manuscript.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

This research is financially supported by the Ministry of Science and Technology of the Republic of China, Taiwan under Contract numbers of MOST 105-2221-E-010-014-MY2, 103-2320-B-010-014-MY2, 103-2321-B-010-016, 102-2221-E-010-011-MY3, 103-2515-S-035-002 and National Health Research Institutes grant of BN-105-PP-15 and the Zhenjiang University, China under the Fund number of 181110- 193544B01/007.

Supplementary Material

The Supplementary Material for this article can be found online at: http://journal.frontiersin.org/article/10.3389/fnins.2016.00556/full#supplementary-material

References

Aha, D. W., and Bankert, R. L. (1996). “A comparative evaluation of sequential feature selection algorithms,” in Learning from Data, eds D. Fisher and H. J. Lenz (New York, NY: Springer-Verlag), 199–206.

Ashe, J., and Georgopoulos, A. P. (1994). Movement parameters and neural activity in motor cortex and area 5. Cereb. Cortex 4, 590–600. doi: 10.1093/cercor/4.6.590

PubMed Abstract | CrossRef Full Text | Google Scholar

Bishop, W., Chestek, C. C., Gilja, V., Nuyujukian, P., Foster, J. D., Ryu, S. I., et al. (2014). Self-recalibrating classifiers for intracortical brain–computer interfaces. J. Neural Eng. 11:026001. doi: 10.1088/1741-2560/11/2/026001

PubMed Abstract | CrossRef Full Text | Google Scholar

Brown, E. N., Frank, L. M., Tang, D., Quirk, M. C., and Wilson, M. A. (1998). A statistical paradigm for neural spike train decoding applied to position prediction from ensemble firing patterns of rat hippocampal place cells. J. Neurosci. 18, 7411–7425.

PubMed Abstract | Google Scholar

Carmena, J. M., Lebedev, M. A., Crist, R. E., O'Doherty, J. E., Santucci, D. M., Dimitrov, D. F., et al. (2003). Learning to control a brain-machine interface for reaching and grasping by primates. PLoS Biol. 1:e42. doi: 10.1371/journal.pbio.0000042

PubMed Abstract | CrossRef Full Text | Google Scholar

Chapin, J. K., Moxon, K. A., Markowitz, R. S., and Nicolelis, M. A. (1999). Real-time control of a robot arm using simultaneously recorded neurons in the motor cortex. Nat. Neurosci. 2, 664–670. doi: 10.1038/10223

PubMed Abstract | CrossRef Full Text | Google Scholar

Churchland, M. M., Santhanam, G., and Shenoy, K. V. (2006). Preparatory activity in premotor and motor cortex reflects the speed of the upcoming reach. J. Neurophysiol. 96, 3130–3146. doi: 10.1152/jn.00307.2006

PubMed Abstract | CrossRef Full Text | Google Scholar

Collinger, J. L., Wodlinger, B., Downey, J. E., Wang, W., Tyler-Kabara, E. C., Weber, D. J., et al. (2013). High-performance neuroprosthetic control by an individual with tetraplegia. Lancet 381, 557–564. doi: 10.1016/S0140-6736(12)61816-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Cunningham, J. P., Nuyujukian, P., Gilja, V., Chestek, C. A., Ryu, S. I., and Shenoy, K. V. (2011). A closed-loop human simulator for investigating the role of feedback control in brain-machine interfaces. J. Neurophysiol. 105, 1932–1949. doi: 10.1152/jn.00503.2010

PubMed Abstract | CrossRef Full Text | Google Scholar

Dai, J. J., Lieu, L., and Rocke, D. (2006). Dimension reduction for classification with gene expression microarray data. Stat. Appl. Genet. Mol. Biol. 5, 6. doi: 10.2202/1544-6115.1147

PubMed Abstract | CrossRef Full Text | Google Scholar

Dangi, S., Gowda, S., and Carmena, J. M. (2013). Likelihood gradient ascent (LGA): a closed-loop decoder adaptation algorithm for brain-machine interfaces. Conf. Proc. IEEE Eng. Med. Biol. Soc. 2013, 2768–2771. doi: 10.1109/embc.2013.6610114

PubMed Abstract | CrossRef Full Text | Google Scholar

Dethier, J., Nuyujukian, P., Ryu, S. I., Shenoy, K. V., and Boahen, K. (2013). Design and validation of a real-time spiking-neural-network decoder for brain–machine interfaces. J. Neural Eng. 10:036008. doi: 10.1088/1741-2560/10/3/036008

PubMed Abstract | CrossRef Full Text | Google Scholar

Donoghue, J. P. (2002). Connecting cortex to machines: recent advances in brain interfaces. Nat. Neurosci. 5, 1085–1088. doi: 10.1038/nn947

PubMed Abstract | CrossRef Full Text | Google Scholar

Efron, B., and Tibshirani, R. J. (1994). An Introduction to the Bootstrap. New York, NY: CRC press.

Google Scholar

Georgopoulos, A. P., Kettner, R. E., and Schwartz, A. B. (1988). Primate motor cortex and free arm movements to visual targets in three-dimensional space. II. Coding of the direction of movement by a neuronal population. J. Neurosci. 8, 2928–2937.

PubMed Abstract | Google Scholar

Gilja, V., Chestek, C., Diester, I., Henderson, J. M., Deisseroth, K., and Shenoy, K. V. (2011). Challenges and opportunities for next-generation intracortically based neural prostheses. IEEE Trans. Biomed. Eng. 58, 1891–1899. doi: 10.1109/TBME.2011.2107553

PubMed Abstract | CrossRef Full Text | Google Scholar

Gilja, V., Nuyujukian, P., Chestek, C. A., Cunningham, J. P., Byron, M., Yu, B. M., et al. (2012). A high-performance neural prosthesis enabled by control algorithm design. Nat. Neurosci. 15, 1752–1757. doi: 10.1038/nn.3265

PubMed Abstract | CrossRef Full Text | Google Scholar

Gilja, V., Pandarinath, C., Blabe, C. H., Nuyujukian, P., Simeral, J. D., Sarma, A. A., et al. (2015). Clinical translation of a high-performance neural prosthesis. Nat. Med. 21, 1142–1145. doi: 10.1038/nm.3953

PubMed Abstract | CrossRef Full Text | Google Scholar

Hochberg, L. R., Bacher, D., Jarosiewicz, B., Masse, N. Y., Simeral, J. D., Vogel, J., et al. (2012). Reach and grasp by people with tetraplegia using a neurally controlled robotic arm. Nature 485, 372–375. doi: 10.1038/nature11076

PubMed Abstract | CrossRef Full Text | Google Scholar

Hochberg, L. R., Serruya, M. D., Friehs, G. M., Mukand, J. A., Saleh, M., Caplan, A. H., et al. (2006). Neuronal ensemble control of prosthetic devices by a human with tetraplegia. Nature 442, 164–171. doi: 10.1038/nature04970

PubMed Abstract | CrossRef Full Text | Google Scholar

Kennedy, P. R., Bakay, R. A. E., Moore, M. M., Adams, K., and Goldwaithe, J. (2000). Direct control of a computer from the human central nervous system. IEEE Trans. Rehabil. Eng. 8, 198–202. doi: 10.1109/86.847815

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, K. C. (1991). Sliced inverse regression for dimension reduction. J. Am. Stat. Assoc. 86, 316–327. doi: 10.1080/01621459.1991.10475035

CrossRef Full Text | Google Scholar

Lin, H. C., Pan, H. C., Lin, S. H., Lo, Y. C., Shen, E. T.-H., Liao, L.-D., et al. (2016). Central thalamic deep-brain stimulation alters striatal-thalamic connectivity in cognitive neural behavior. Front. Neural Circuits. 9:87. doi: 10.3389/fncir.2015.00087

PubMed Abstract | CrossRef Full Text | Google Scholar

Lu, H. H. S. (2008). “Reconstruction, visualization and analysis of medical images,” in Handbook of Computational Statistics: Vol. 3, Data Visualization, eds C. H. W. Chen, W. Hardle, and A. Unwin (Berlin; Heidelberg: Springer-Verlag), 813–830.

Moran, D. W., and Schwartz, A. B. (1999). Motor cortical representation of speed and direction during reaching. J. Neurophysiol. 82, 2676–2692.

PubMed Abstract | Google Scholar

Naik, P. A., Hagerty, M. R., and Tsai, C. L. (2000). A new dimension reduction approach for data-rich marketing environments: sliced inverse regression. J. Marketing. Res. 37, 88–101. doi: 10.1509/jmkr.37.1.88.18715

CrossRef Full Text | Google Scholar

Nicolelis, M. A. L. (2001). Actions from thoughts. Nature 409, 403–408. doi: 10.1038/35053191

PubMed Abstract | CrossRef Full Text | Google Scholar

Paninski, L., Fellows, M. R., Hatsopoulos, N. G., and Donoghue, J. P. (2004). Spatiotemporal tuning of motor cortical neurons for hand position and velocity. J. Neurophysiol. 91, 515–532. doi: 10.1152/jn.00587.2002

PubMed Abstract | CrossRef Full Text | Google Scholar

Pohlmeyer, E. A., Mahmoudi, B., Geng, S., Prins, N. W., and Sanchez, J. C. (2014). Using reinforcement learning to provide stable brain-machine interface control despite neural input reorganization. PLoS ONE 9:e87253. doi: 10.1371/journal.pone.0087253

PubMed Abstract | CrossRef Full Text | Google Scholar

Polikov, V. S., Tresco, P. A., and Reichert, W. M. (2005). Response of brain tissue to chronically implanted neural electrodes. J. Neurosci. Methods 148, 1–18. doi: 10.1016/j.jneumeth.2005.08.015

PubMed Abstract | CrossRef Full Text | Google Scholar

Reina, G. A., Moran, D. W., and Schwartz, A. B. (2001). On the relationship between joint angular velocity and motor cortical discharge during reaching. J. Neurophysiol. 85, 2576–2589.

PubMed Abstract | Google Scholar

Sadtler, P. T., Quick, K. M., Golub, M. D., Chase, S. M., Ryu, S. I., Tyler-Kabara, E. C., et al. (2014). Neural constraints on learning. Nature 512, 423–426. doi: 10.1038/nature13665

PubMed Abstract | CrossRef Full Text | Google Scholar

Salinas, E., and Abbott, L. F. (1994). Vector reconstruction from firing rates. J. Comput. Neurosci. 1, 89–107. doi: 10.1007/BF00962720

PubMed Abstract | CrossRef Full Text | Google Scholar

Schwartz, A. B. (1993). Motor cortical activity during drawing movements: population representation during sinusoid tracing. J. Neurophysiol. 70, 28–36.

PubMed Abstract | Google Scholar

Schwartz, A. B. (1994). Direct cortical representation of drawing. Science 265, 540–542. doi: 10.1126/science.8036499

PubMed Abstract | CrossRef Full Text | Google Scholar

Schwartz, A. B., Cui, X. T., Weber, D. J., and Moran, D. W. (2006). Brain-controlled interfaces: movement restoration with neural prosthetics. Neuron 52, 205–220. doi: 10.1016/j.neuron.2006.09.019

PubMed Abstract | CrossRef Full Text | Google Scholar

Schwartz, A. B., Taylor, D. M., and Tillery, S. I. (2001). Extraction algorithms for cortical control of arm prosthetics. Curr. Opin. Neurobiol. 11, 701–708. doi: 10.1016/S0959-4388(01)00272-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Serruya, M. D., Hatsopoulos, N. G., Paninski, L., Fellows, M. R., and Donoghue, J. P. (2002). Brain-machine interface: instant neural control of a movement signal. Nature 416, 141–142. doi: 10.1038/416141a

PubMed Abstract | CrossRef Full Text | Google Scholar

Shenoy, K. V., and Carmena, J. M. (2014). Combining decoder design and neural adaptation in brain-machine interfaces. Neuron 84, 665–680. doi: 10.1016/j.neuron.2014.08.038

PubMed Abstract | CrossRef Full Text | Google Scholar

Srinivasan, L., and da Silva, M. (2011). Breaking the fixed-arrival-time restriction in reaching movements of neural prosthetic devices. IEEE Trans. Biomed. Eng. 58, 1555–1564. doi: 10.1109/TBME.2010.2101599

PubMed Abstract | CrossRef Full Text | Google Scholar

Sussillo, D., Nuyujukian, P., Fan, J. M., Kao, J. C., Stavisky, S. D., Ryu, S., et al. (2012). A recurrent neural network for closed-loop intracortical brain–machine interface decoders. J. Neural Eng. 9:026027. doi: 10.1088/1741-2560/9/2/026027

PubMed Abstract | CrossRef Full Text | Google Scholar

Takeda, K., and Funahashi, S. (2004). Population vector analysis of primate prefrontal activity during spatial working memory. Cereb. Cortex 14, 1328–1339. doi: 10.1093/cercor/bhh093

PubMed Abstract | CrossRef Full Text | Google Scholar

Taylor, D. M., Tillery, S. I., and Schwartz, A. B. (2002). Direct cortical control of 3D neuroprosthetic devices. Science 296, 1829–1832. doi: 10.1126/science.1070291

PubMed Abstract | CrossRef Full Text | Google Scholar

Tu, Y., Tan, A., Fu, Z., Hung, Y. S., Hu, L., and Zhang, Z. (2015). Supervised nonlinear dimension reduction of functional magnetic resonance imaging data using Sliced Inverse Regression. Conf. Proc. IEEE Eng. Med. Biol. Soc. 2015, 2641–2644. doi: 10.1109/EMBC.2015.7318934

PubMed Abstract | CrossRef Full Text | Google Scholar

Valle, S., Li, W., and Qin, S. J. (1999). Selection of the number of principal components: the variance of the reconstruction error criterion with a comparison to other methods. Ind. Eng. Chem. Res. 38, 4389–4401. doi: 10.1021/ie990110i

CrossRef Full Text | Google Scholar

Velliste, M., Perel, S., Spalding, M. C., Whitford, A. S., and Schwartz, A. B. (2008). Cortical control of a prosthetic arm for self-feeding. Nature 453, 1098–1101. doi: 10.1038/nature06996

PubMed Abstract | CrossRef Full Text | Google Scholar

Wahnoun, R., He, J., and Helms-Tillery, S. I. (2006). Selection and parameterization of cortical neurons for neuroprosthetic control. J. Neural Eng. 3, 162–171. doi: 10.1088/1741-2560/3/2/010

PubMed Abstract | CrossRef Full Text | Google Scholar

Watanabe, Y., Takeda, K., and Funahashi, S. (2009). Population vector analysis of primate mediodorsal thalamic activity during oculomotor delayed-response performance. Cereb. Cortex 19, 1313–1321. doi: 10.1093/cercor/bhn170

PubMed Abstract | CrossRef Full Text | Google Scholar

Wessberg, J., Stambaugh, C. R., Kralik, J. D., Beck, P. D., Laubach, M., Chapin, J. K., et al. (2000). Real-time prediction of hand trajectory by ensembles of cortical neurons in primates. Nature 408, 361–365. doi: 10.1038/35042582

PubMed Abstract | CrossRef Full Text | Google Scholar

Wodlinger, B., Downey, J., Tyler-Kabara, E., Schwartz, A., Boninger, M., and Collinger, J. (2015). Ten-dimensional anthropomorphic arm control in a human brain−machine interface: difficulties, solutions, and limitations. J. Neural Eng. 12:016011. doi: 10.1088/1741-2560/12/1/016011

PubMed Abstract | CrossRef Full Text | Google Scholar

Wold, S., Esbensen, K., and Geladi, P. (1987). Principal component analysis. Chemometr. Intell. Lab. Syst. 2, 37–52. doi: 10.1016/0169-7439(87)80084-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Wu, H. M. (2008). Kernel sliced inverse regression with applications to classification. J. Comput. Graph. Stat. 17, 590–610. doi: 10.1198/106186008X345161

CrossRef Full Text | Google Scholar

Wu, H. M., and Lu, H. H. S. (2007). Iterative sliced inverse regression for segmentation of ultrasound and MR images. Pattern Recognit. 40, 3492–3502. doi: 10.1016/j.patcog.2007.04.019

CrossRef Full Text | Google Scholar

Wu, W., Black, M. J., Mumford, D., Gao, Y., Bienenstock, E., and Donoghue, J. P. (2004). Modeling and decoding motor cortical activity using a switching Kalman filter. IEEE Trans. Biomed. Eng. 51, 933–942. doi: 10.1109/TBME.2004.826666

PubMed Abstract | CrossRef Full Text | Google Scholar

Yu, B. M., Kemere, C., Santhanam, G., Afshar, A., Ryu, S. I., Meng, T. H., et al. (2007). Mixture of trajectory models for neural decoding of goal-directed movements. J. Neurophysiol. 97, 3763–3780. doi: 10.1152/jn.00482.2006

PubMed Abstract | CrossRef Full Text | Google Scholar