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.

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.

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.

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.

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.

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.