Penalized Multi-Way Partial Least Squares for Smooth Trajectory Decoding from Electrocorticographic (ECoG) Recording

In the current paper the decoding algorithms for motor-related BCI systems for continuous upper limb trajectory prediction are considered. Two methods for the smooth prediction, namely Sobolev and Polynomial Penalized Multi-Way Partial Least Squares (PLS) regressions, are proposed. The methods are compared to the Multi-Way Partial Least Squares and Kalman Filter approaches. The comparison demonstrated that the proposed methods combined the prediction accuracy of the algorithms of the PLS family and trajectory smoothness of the Kalman Filter. In addition, the prediction delay is significantly lower for the proposed algorithms than for the Kalman Filter approach. The proposed methods could be applied in a wide range of applications beyond neuroscience.


Introduction
The Brain-Computer Interface (BCI) is a system for converting the brain's neural activity into commands for external devices [1]. Motor-related BCI aims at improving the quality of life of individuals with severe motor disabilities [2]. Paralyzed or handicapped persons get the opportunity for mental control of wheelchairs, robotic arms, exoskeletons, etc. [3][4][5][6][7][8]. Besides that, motor-related BCI systems can be applied for rehabilitation purposes [9,10].
In this paper, the task of reconstructing continuous 3D trajectories of the subject's hand from ECoG recordings is considered. This problem was studied in both monkeys and humans [22,[29][30][31]. A set of approaches for trajectory decoding from neuronal recordings were proposed [22][23][24]28]. A common approach consists in extraction of informative features from spatial [32], frequency [33], or/and temporal [34] domains of brain signal. Decoder is identified from the trajectory information (response variable) and neural activity features (exploratory variable). Many standard methods for model identification from experimental data are designed for vector exploratory variables which generally represent only one domain (modality) of analysis [26,[35][36][37]. If analysis in one domain does not provide satisfactory results, several modalities could be combined sequentially in this paradigm. A tensor-based (Multi-Way) approach is intended for simultaneous treatment of several domains. Multi-electrodes brain signals recording is mapped to the spatial-temporal-frequency space [38] representing single time epoch by a cube (a third-order tensor). All cubes are stored in a tensor of observations (a fourth-order tensor). Tensor (multi-way array) data representation is convenient, since it save natural structure of the data. Detailed review of the tensor data representation as well as the tensor analysis is given in [39]. Multi-Way Analysis is reported as an efficient tool for the decoder identification [22,40,41].
The peculiarity of the tensor data analysis is their high dimension. The algorithms of the Partial Least Squares (PLS) regression family are known to be particularly appropriate for high dimension tasks [42]. PLS projects both exploratory and response variables to low dimensional spaces of latent variables. Contrary to the widespread Principal Component Analysis (PCA) [43], which projects only explanatory variables, both explanatory and response variables are taken into account by the PLS for projectors identification. Thus, the PLS provides projectors correlated to the task. However, the generic PLS approach is vector oriented and can be applied only for the tensors unfolded to the vectors (Unfolded PLS) [44]. To exploit the advantages of tensor data representation, generalizations of the PLS for the tensor data were invented (N-Way PLS [45], Higher Order PLS [40], etc.) and were used in BCI systems [40,46,47].
The smoothness of the predicted trajectory is an important property of motor-related BCI systems [29,48,49]. The particular importance of smoothness was reported for control of such BCI effectors as wheelchairs [50] and cars [51]. There are two main approaches to manage the problem of a smooth prediction. The first one consists in the smoothing post-processing of the predicted trajectories [29]. However, this post-processing increases the system's latency, depending on the required level of smoothness, the decision rate of the system, etc. Another solution consists in identifying the model, for which the smoothness of the predicted trajectories is an intrinsic property. A Kalman Filter (KF) [52] is used in BCI systems to generate smooth trajectories [53,54]. However, the method is not adjusted for the high dimensional data due to the application of the matrix inversion operation in it. Moreover, the Kalman Filter weights the current prediction with the previous one, which also brings a latency to the system. Minimal system latency is an essential requirement for the comfortable use of real-time BCI systems [55]. Improvement of the prediction smoothness with minimal supplementary delay is the particular goal of the article.
In this paper, a new penalization approach is proposed to improve neural decoding smoothness. By means of penalization, one can impose different penalties related to the particular tasks. The paper novelty resides in two new ways of penalization (new penalization terms), which particularly address solution smoothness. Namely, Sobolev Penalized Multi-Way PLS (SNPLS) and Polynomial Penalized Multi-Way PLS (PNPLS), are proposed. They combine tensor representation of the data with the possibility to control the level of the smoothness of the predicted trajectories. In these methods, the smoothness is not provided by the post-processing or weighting with previous predictions, but is an internal parameter of the identified model. Thus, no supplementary delay is introduced due to post-processing of the response. In the current paper, a set of state of the art BCI methods, namely Kalman Filter, generic PLS, and Multi-Way PLS (NPLS), are considered and compared with the proposed SNPLS and PNPLS on the set of publicly available ECoG recordings. In addition, response latencies are investigated and compared for the proposed approach and state of the art methods.
The proposed methods could be applied in a wide range of applications beyond neuroscience.

Prediction performance evaluation
Different criteria for performance evaluation have been applied in BCI studies [26,29,40,56]. In this paper, Pearson correlation (r), Root Mean Squares Error (RMSE), and Mean Absolute Error (MAE) are used to compare the methods. A shortcoming of the above criteria is that they do not reflect the smoothness of the predicted trajectory. They are based on the sum of squares/absolute values of residuals which temporal sequence is not taken into account. Mean Absolute Differential Error (MADE) is a way to assess both prediction accuracy and smoothness [29]. Hence MADE is included to the list of criteria for performance evaluation. The Pearson correlation is a commonly used criterion [22,26,30,46] for the continuedvalue trajectories decoding, Where y = (y(t 1 ),. . .,y(t N )) T = (y 1 ,. . .,y N ) T andŷ ¼ ðŷðt 1 Þ; . . . ;ŷðt N ÞÞ T ¼ ðŷ 1 ; . . . ;ŷ N Þ T are the observed and predicted trajectories. RMSE represents the L 2 -norm distance between the predicted and observed trajectories [22,29,40]: Less sensitive to the outliers, MAE [29,56] is based on the L 1 -norm distance: MADE criterion compares derivatives y' andŷ 0 : MADE ¼ ky 0 Àŷ 0 k 1 =ky 0 À y 0 k 1 ; MADE characterizes the prediction smoothness [29] and, being based on the L 1 -norm, is robust to the outliers.

Kalman Filter
The KF [52,57] is a vector-oriented method, which recursively operates on streams of input data to produce an estimate of the underlying system state. KF is widely used to decode the neural activity in BCIs [53,[58][59][60][61]. One advantage of KF is the smoothness of the prediction in time [53,59]. Thus, KF was applied for prediction of the hand's trajectory from the brain's neuronal activity, such as ECoG [60] and spikes [53,61]. Brain activity data x k 2 < n , observed at the moment t k , are used to predict the system state y k 2 < m (m = 9 represents 3D coordinates, velocity, and acceleration of the hand at the moment t k ). The generative model of the KF is formulated as a linear dependency between the brain activity and system state: where H k 2 < nÂm is a matrix of linear coefficients, and q k e N ð0; Q k Þ is the zero-mean noise of the observation, Q k 2 < nÂn . The model expects that the system states are linearly propagated in time: y k+1 = A k y k +w k , where A k 2 < mÂm is a matrix of linear coefficients, w k e N ð0; W k Þ, W k 2 < mÂm . If matrices H k = H, A k = A, Q k = Q, and W k = W are constants in time, they can be estimated from the training data X ¼ ðx 1 ; . . . ; x N Þ 2 < nÂN , Y ¼ ðy 1 ; . . . ; y N Þ 2 < mÂN by means of the least squares: The KF algorithm consists of a priori and a posteriori steps. In the first step, the method a priori estimates the state of the system asŷ À k ¼ Aŷ kÀ1 . In the second step, a posteriori estimationŷ k is found as a linear combination of the a priori estimation and the weighted difference between the observed and predicted signals: Here is the a priori estimate error, and P k ¼ ðI À K k HÞP À k is the a posteriori estimate error [52]. More detailed information about the KF can be found in [57].

Generic PLS
A linear regression model between the response variable y 2 < m and the explanatory variable

the Ordinary Least Squares (OLS) gives [63]
: However, the OLS procedure is unstable for the case of the high dimension task [63]. Partial Least Squares regression is a statistical method for linear vector-based data analysis [42]. Due to its efficient dimension reduction technique even in the case of highly correlated variables, PLS is particularly appropriate for the high dimensional data. The PLS algorithm iteratively estimates a linear relation between the matrices of observations of input and output variables X 2 < NÂn and Y 2 < NÂm : Y = XB + D, where B 2 < nÂm is the matrix of linear coefficients, and D 2 < NÂm is the matrix of noise. On each iteration, X and Y are projected in low dimension spaces (latent variables) in such a way as to explain the maximum of their variation simultaneously: . . . ; u F 2 < NÂF are the matrix of the latent variables (scores), P ¼ ½p 1 ; . . . ; p F 2 < nÂF and Q ¼ ½q 1 ; . . . ; q F 2 < mÂF are the matrix of the loading vectors (projectors), E and F are residual matrices, and F is the number of iterations (factors). Estimation of the factors number F can be done through the cross-validation procedure [64]. A detailed description of the PLS approach can be found in [42].
The PLS regression was applied in the BCI systems for hand trajectory reconstruction [22,26,40].

Multi-way PLS
Multi-way (N-way) PLS is a generalization of the PLS method to the tensor input/output data [45]. Tensors (multi-way arrays) allow the representation of data of a high order. Vectors and matrices are tensors of orders one and two, respectively. In this paper, a tensor is denoted as X 2 < I 1 Â...ÂI n , where n is the order of the tensor. More information about the tensors can be found in [39]. NPLS inherits the robustness of PLS and allows simultaneous analysis of data in several domains (e.g. time-frequency-space). Examples of the application of multi-way tensorbased methods in BCI research are given in [40,46,[65][66][67].
Similarly to PLS regression, NPLS builds the linear relation between the dependent and independent tensors of observations X 2 < NÂI 1 Â...ÂI n and Y 2 < NÂJ 1 Â...ÂJ m . The model is constructed iteratively by projection of the tensors to the space of latent variables: where the operator "°" is the outer product [39], NÂf are the matrices of the latent variables after f = 1,. . .,F iterations, w f i 2 < I i and q f j 2 < J j are the projection vectors, B f is the matrix of the linear coefficients, and E and F are the residual tensors.
In contrast to the OLS approach, NPLS, as well as PLS, is a biased estimator. However, it provides the insensitivity of the predictive model to ill-conditioning and noise [68].

Penalized regression
Modification of the cost function (optimization functional) in the model identification procedure could be used to improve the prediction/regression properties. Among other reasons, penalizations are widely applied to achieve desirable characteristics (sparsity, robustness, etc.). Additional terms are added to the optimization task, basically performing the standard L 1norm or L 2 -norm optimization. For linear regression E[y|x] = Bx the least squares estimation is given in (7). To impose the additional restrictions, the regularization (penalization) term is introduced [69][70][71][72][73][74]. (RIDGE) are the most frequently used, λ ! 0 represents the penalization parameter. Whereas the P L 1 penalization encourages sparse solution (9), the P L 2 penalization reduces the variation of the coefficients B, which increases the robustness of the model [71]. In addition, the P L 2 penalization as well as Fused LASSO [75] are often applied to improve the smoothness of the coefficients B [70,73,74]. However, neither P L 1 nor P L 2 penalizations improve the smoothness of prediction of y(t) analyzing the data flow x(t). In this paper, we propose two penalizations to improve the smoothness of the continuous trajectories prediction.
SNPLS. The Sobolev NPLS method is based on the use of the Sobolev norm W s [76] instead of the L 2 -norm in the optimization task: Here ϕ and ϕ (s) are a function and its s-th derivation, respectively. Optimization in Sobolev space corresponds to the optimization problem (9), where the penalization term is added with the parameter λ ! 0 to control the influence of the derivative part: The problem can be represented in the matrix form: where The optimization problem (12) corresponds to the regression estimation task with the matrices of independent and dependent variables e X and e Y, respectively. Moreover, for the tensor case, the same approach can be applied to obtain the tensors e X and e Y. To identify the regression between the tensors e X and e Y, we use the NPLS algorithm, since it is suitable for tensor-based variables and has demonstrated its efficiency for BCI tasks.
PNPLS. Polynomial Penalized NPLS considers the optimization problem: where x i ¼ ðx 1 i ; . . . ; x n i Þ T , and P p;ðx iÀl ;...;x iþl Þ ðÁÞ is a polynomial function of the power p, which provides the minimal least squares approximation of the set of points (x i-l ,. . .,x i+l ), x _ i ¼ P p;ðx iÀl ;...;x iþl Þ ðx i Þ. The procedure of x _ identification can be represented as a result of the power p polynomial filtration in a sliding window of size (2l + 1), where polynomial coefficients are fitted independently for each window.
In the matrix notation, the optimization problem is represented as Penalizing the difference between the results of application of the model B to the initial (X) and "smoothed" (X _ ) explanatory variables, this algorithm is looking for the model which provides similar results being applied to the original and "smoothed" data.
For the new variables e X and e Y: Similar to (12), the task (15) corresponds to the problem of regression estimation with the matrix of the explanatory and response variables e X and e Y. It can also be generalized to the case of the tensors e X and e Y. The NPLS algorithm is used for the following regression estimation.
Both SNPLS and PNPLS algorithms improve solutions without increasing the complexity of the optimization problem (sum-of-squares optimization).

Evaluation and comparison
The BCI problem of the reconstruction of 3D trajectories of a monkey hand from its ECoG brain activity is considered to evaluate the accuracy and smoothness of the proposed approaches. Two proposed methods (SNPLS and PNPLS) are evaluated and compared to three state-of-the-art methods (Kalman Filter, generic PLS, and tensor-based NPLS).
The publicly available dataset recorded and distributed by the Laboratory for Adaptive Intelligence, BSI, RIKEN (http://neurotycho.org) was used to evaluate the methods. During the experiments, the 3D position of the monkey's right hand was recorded simultaneously with the epidural ECoG signal of the brain [26,77]. The dataset used in the current paper consists of 20 recordings, corresponding to two Japanese macaques (10 recordings of each monkey), denoted as 'B' and 'C' in the database [22]. The monkeys were trained to receive pieces of food with their right hands. The pieces were demonstrated by the experimenter at random locations at a distance of about 20 cm from the monkeys. The new trial started when the monkey finished eating the previous piece of food. The ECoG data were acquired with 64 electrodes (Blackrock Microsystems, Salt Lake City, UT, USA) implanted in the epidural space of the left hemisphere (see Fig 1). A sampling rate of 1000 Hz per channel was provided. To record the hand trajectories, an optical motion capture system was used (Vicon Motion Systems, Oxford, UK) with a sampling rate of 120 Hz. The experiment scheme is represented in Fig 2A. Additional description of the experiments and analyzed data is given in [22,77].
The duration of each of 20 experiment sessions was 15 minutes. The first 10 minutes of each recording were used for model identification and the last 5 minutes were employed for testing of the corresponding model.
Feature extraction. The input data tensor X was formed from the ECoG epochs. Each epoch contains 1 second of the signal, mapped to the temporal-spatial-frequency space by The scheme of the experiment. The monkey is following the food with its hand. 3D coordinates of the hand are recorded simultaneously with monkey ECoG brain activity [22]. (B) For each moment t, the response variable y(t) is formed from the 3D coordinates of the monkey's hand at this moment. To form the explanatory variable x(t), the ECoG epoch [t-Δτ, t], Δτ = 1 s is mapped by the Continuous Wavelet Transform (CWT). Then, the logarithm of the absolute values of the wavelet coefficients is taken, the procedure of the artifact filtration is applied, and the tensor is decimated along the temporal modality 100 times (1000 points to 10 points). Then the whole procedure is repeated for the next time moment t i + 1 = t i + Δt, Δτ = 0.1 s. means of the wavelet transform. The complex Morlet wavelet was chosen as a mother wavelet due to its widespread application for BCI data analysis [22,26,30,40,46]. The epochs were taken continuously with time shifts equal to 100 ms. Based on [29], frequencies from 10 to 150 Hz with a step of 10 Hz were chosen. The logarithm of the absolute values of the wavelet coefficients was taken. After formation of the input tensor X, the chewing artifacts [22] were filtered in the way described in [29]. Then, the feature tensor was decimated along the temporal modality 100 times, by averaging the data in 10 sliding windows of 100 ms length. A detailed description of the analyzed data formation procedure is given in [29]. The matrix of the response variables Y was formed from the corresponding 3D coordinates of the monkey's hand. Whereas for PLS, NPLS, SNPLS, PNPLS methods matrix Y represents 3D coordinates, the Kalman Filter uses the extended matrix of the response variables which additionally includes velocity and acceleration Y KF = [Y, Y', Y'']. The explanatory variables tensor X is identical for all the methods. Fig 2B represents the general scheme of the data preparation procedure.
To validate the methods, the explanatory variables tensor X and the matrix of the response variables Y were split into training (70%) and test (30%) subsets for each recording: for test data Parametrization. Three state-of-the-art methods (Kalman Filter, generic PLS, and tensorbased NPLS) and two proposed tensor-based approaches (SNPLS and PNPLS) were compared using a set of optimized parameters. The optimal number of factors F (PLS, NPLS, SNPLS, PNPLS) as well as the smoothing parameter λ (SNPLS, PNPLS) were chosen on the training set by the 10-fold cross-validation procedure [64] for each recording. The preliminary study on one recording gave the polynomic parameters for the PNPLS algorithm: the number of points chosen for the polynomial coefficients estimation was equal to L = 2l + 1 = 9 and the polynomic power was p = 2. Then, these parameters were fixed for all other recordings. Similarly, the optimal derivative order chosen for the SNPLS approach was s = 3.
Results. The models were calibrated separately on each training recording and tested on the corresponding testing data. An example of the application of the compared methods to the same interval of data is illustrated in   Table 1 represents the p-values for all the pairs of the methods. In the present study, the significance level of α = 0.05 is chosen. The algorithms of the PLS family, and in particular the proposed SNPLS and PNPLS, significantly outperform the KF in correlation, RMSE, and MAE. The difference between the generic methods (PLS, NPLS) and the proposed approaches (SNPLS, PNPLS) for these criteria is insignificant (Table 1) and does not exceed 5% (Fig 4). In general, the presented prediction accuracy corresponds to the accuracy reported for this database [22,26,29,46].
KF, SNPLS and PNPLS demonstrated better smoothness than PLS and NPLS, which is reflected by the MADE criterion: the improvement is significant (Table 1) and varies from 30% to 47% (Fig 4). At the same time, although KF showed the best smoothness, the difference between KF and the proposed PNPLS approach is less than 6% and is insignificant. To illustrate PLS-family predictive models, the averaged influence [78] of the frequency, temporal, as well as spatial modalities estimated for the monkey "B" are shown in Fig 5 (the results for the monkey "C" are consistent). For each modality, the influences of its elements are evaluated as the weights of sum of the absolute values of the related predictive model's coefficients.

Prediction delay comparison
System latency is an important characteristic for BCIs [55]. Smoothing generally introduce the delay to the prediction. An additional study was carried out to assess and compare prediction delays of the analyzed methods.
To compare latencies, the predicted trajectories generated by the analyzed approaches were shifted in time to maximize the correlation between predicted and observed trajectories. Prediction delays were estimated as minimal shifts of the predicted signals, which maximizes the correlation between predicted and tested trajectories. The shift interval S = [0, 2] s was studied.

Discussion
Motor-related BCI systems are of great importance since they can help people with severe motor disabilities to recover their functionality. However, decoding the neural activity is a  challenging task [1]. BCI control of external devices like wheelchairs, artificial arms, exoskeletons, etc. in real life is especially difficult, due to safety issues [24]. This imposes additional restrictions on BCI decoders, as well as on the smoothness of the predicted trajectories, in particular. The standard approach for improving smoothness consists of post-filtering of the obtained prediction [55,79]. Such filters (e.g. sliding window averaging) bring a delay to the system response. More complex nonlinear filters, such as KF, are labor-consuming, which decreases the decision rate of real-time BCI systems, especially in the case of high dimensional data, and also introduces a prediction delay. Both a high decision rate and minimal system latency are essential requirements for the comfortable use of real-time BCI systems [55]. The methods proposed in the current paper, namely SNPLS and PNPLS, are based on multi-way data analysis, which is known to be efficient for neural data treatment [40,47,80]. These methods construct the decoders with the possibility of controlling the smoothness of the predicted trajectories. The proposed algorithms try to improve the smoothness of prediction without increasing the system latency. Penalization integrated into the linear decoder identification procedure can be interpreted, as the penalization of informative but noisy variables. The proposed methods improve solutions stability by penalizing unsmooth predictions without increasing the complexity of the optimization problem (sum-of-squares optimization). The proposed methods do not change the predictive model complexity as well. Solutions remains linear, and the resulted model can be efficiently applied for high dimensional data flow decoding in real-time on standard computer (Intel Dual Core, 3.16 GHz; RAM 3.25 Gb) [80]. The proposed methods could be applied in a wide range of applications beyond neuroscience.
To validate the efficiency of SNPLS and PNPLS for BCI tasks, the algorithms were compared with a set of state-of-the-art approaches. The Kalman Filter has been chosen as a method which generates a smooth prediction and is widely applied in BCIs. However, it contains the inversion of the matrix. This limits the application of the KF for high dimensional tasks [53]. Unlike the KF, the algorithms of the Partial Least Squares family are particularly suited for the high dimensional data and are also often used in BCI systems [22,46]. Improvement of the prediction smoothness of these methods, with minimal loss of accuracy as well as supplementary delay, was the particular goal of the article. The prediction accuracy of the compared methods was evaluated by means of Pearson correlation, RMSE, and the more robust to outliers MAE criteria. To evaluate the smoothness, the MADE criterion was chosen. It combines the robustness of the L 1 -norm and allows assessment of the smoothness.
The algorithms of the PLS family significantly outperform the KF in prediction accuracy estimated with Correlation, RMSE, MAE, whereas the difference among the PLS-based approaches is insignificant. Due to matrix inversion, KF could be more sensitive to high dimensional and noisy explanatory variable in comparison to PLS-family algorithms. Additional preand/or post-processing of data or Kalman algorithm's modifications could improve its prediction accuracy for the current task and this is the subject for future researches.
From the smoothness point of view, KF, PNPLS, and SNPLS significantly outperform generic PLS and NPLS. KF demonstrates the best smoothness. However, the results of KF and new methods are comparable. Thus, PNPLS and SNPLS algorithms combined the prediction accuracy of the PLS-family algorithms with the smoothness approaching the KF. In addition, SNPLS and PNPLS provide lower prediction delay in comparison with the KF. The KF approach latency reaches half a second (470±210 ms) that could be noticeable and irritable for subjects. As it was reported in [79] for spikes decoding methods in monkeys experiments, whereas delays up to 200 ms provide minimal influence, every 100 ms of additional delay added about 200 ms to the reaching time. The studies in human also demonstrated that delays about 200 ms slow down reaching tasks and decrease accuracy [81]. For simulated human experiments, increasing of average reaching time by over 5 s due to visual feedback delay of 300 ms is reported in [55]. Contrary to the KF, the PLS, NPLS, SNPLS, and PNPLS methods provide a latency of about 100 ms that is equal to the decision rate of the system. However, additional real-time closed-loop BCI experiments should be carried out in human to clarify the latency's tolerable level.
Generally, criteria values appropriated for the practical applications are derived according to the subject satisfaction in practical experiments. The influence of different parameters such as error magnitude, processing delays, etc. on closed-loop BCI in human was studied, in particular in [55]. The comparison of decoding approaches during the clinical closed-loop BCI trials will be the next step of our research. Fig 5 demonstrates the modality influence analysis results [78] for the PLS-family methods, in frequency, spatial and temporal domains. In the frequency modality, frequencies influence distributions are similar for all the considered approaches. The distributions have two maximums around 20 Hz and 90 Hz. In the temporal modality, the PLS and NPLS models tend to give close weights to all elements of the analyzing interval, with the maximums near 0 moments. Contrarily, coefficient distribution of the SNPLS and PNPLS models have maximums in the middle of the analyzing intervals. It is possible that weight shapes are responsible for the predicted trajectories smoothness. On the other hand, prediction delay evaluation demonstrated that there is no additional latency due to these weights distributions, as one might expect. Additional study is needed to clarify this effect. In the spatial modality, electrode weights are comparable for all the considered approaches. The most informative electrodes (6, 10, and 13) are the same for all the methods. This could be explained by detecting the readiness potentials, i.e. cortical contribution to the pre-motor planning of volitional movement in the supplementary motor area. Additional information about the readiness potentials could be found in [82].
The influences analysis allows estimating the stability of the prediction models by means of assessment of the weights variability in each modality. Standard deviations through all the modalities are σ PLS = 0.005±0.005, σ NPLS = 0.005±0.008, σ SNPLS = 0.005±0.005, σ PNPLS = 0.004±0.005. The differences of the standard deviations are not significant (α = 0.05). Thus, the proposed PNPLS and SNPLS methods keep the robustness of the generic PLS and NPLS approaches.
Comparing the proposed methods, both SNPLS and PNPLS allow improving of the smoothness of the predicted trajectories. The PNPLS prediction smoothness surpasses the SNPLS one, while the prediction latencies for both approaches are comparable. On the other hand, SNPLS provided slightly a better accuracy (Fig 4). The advantage of SNPLS is that it depends on fewer external parameters than PNPLS, thus it is preferable in the case of restricted computation time for model identification.

Limitations of the present study and future work
The main limitation and disadvantage of the proposed approaches is the time-consuming estimation of a set of the external parameters. The smoothness parameter, the polynomial parameters for PNPLS, the derivative order for SNPLS, as well as the number of factors for SNPLS and PNPLS, need to be evaluated. In this study, a time-consuming cross-validation procedure is applied to identify the external parameters of the algorithms. Additional investigation is needed to optimize parameters estimation procedure.
Parametrization and calibration is the offline procedure. Generally, recalibration is not required for ECoG-based BCIs during a long period [80]. Nevertheless, the fast online calibration would be a strong point for the BCI system. Recently, the Recursive NPLS algorithm was proposed [46] for adaptive BCI system calibration. The combination of the SNPLS and PNPLS methods with RNPLS will allow online adjustment of the decoder to the non-stationary brain activity.
As it was reported in [49,83], methods that demonstrate good performance for offline decoding could be less efficient for online applications, while other methods provide better performance in closed-loop conditions. The next step of the study will be closed-loop application and the evaluation of decoding approaches, and in particular the application of the algorithms in clinical BCI to calibrate the BCI decoder. The fully-implantable device WIMAGINE 1 [25] for chronic measurement and wireless transmission of ECoG data, as well as the full body exoskeleton EMY 1 [84], are currently developed within the framework of the CEA-LETI-CLINA-TEC 1 BCI project. Clinical trials of the exoskeleton control by tetraplegic subjects is in preparation.