Online adaptive group-wise sparse Penalized Recursive Exponentially Weighted N-way Partial Least Square for epidural intracranial BCI

Brain–computer interfaces (BCIs) are systems that create a new communication pathway between the brain and an effector without neuromuscular activation. Motor BCIs aim to allow users with severe motor impairments to regain limb mobility by controlling orthoses and prostheses or by recovering motor control over their own limbs, e.g., using electrical stimulation (Mak and Wolpaw, 2009). Most BCI systems include a neural signal acquisition system, a transducer, an effector, and a feedback system (Schwartz et al., 2006). The transducer is typically composed of a neural feature extraction block and a decoder, with optional pre- and post-processing blocks. Decoders are generally user-specific and data-driven. They are created through a supervised tuning of parameters on the training dataset. Once the decoder is established, the transducer can be applied in real time to translate the user's intention into the control command of the effector.

As BCI systems function in real time, computing time and resource management are crucial aspects of such systems (Haufe et al., 2014). High-resolution neuronal activity recording systems are generally required to achieve high-dimensional control of complex effectors. It results in a large volume of data that needs to be processed.

Furthermore, in motor BCI, a high decision rate (8–10 Hz) is necessary to control complex effectors such as robotic arms, an exoskeleton, and so on (Marathe and Taylor, 2015; Shanechi et al., 2017). In addition to high computing power requirements and computing time, the high dimensionality of the feature space presents challenges such as the “curse of dimensionality” during decoder training (Bellman, 1961; Bishop, 2006; Nicolas-Alonso and Gomez-Gil, 2012; Remeseiro and Bolon-Canedo, 2019). The feature space often contains irrelevant and/or redundant features. Moreover, computational load is critical for the potential development of portable BCIs.

Reducing the dimensionality of the neuronal feature space is one approach to addressing these issues. Dimension reduction algorithms have been widely employed in numerous studies on BCI.

Both projection and feature selection methods were applied for dimensional reduction, in both online BCI experiments and offline analysis.

The feature selection family regroups filter-based, wrapper-based, and embedded techniques (Bolón-Canedo et al., 2013; Khaire and Dhanalakshmi, 2019). The filter-based methods rank and select features independently without considering the decoder. Effective in computation time, these methods tend to select highly correlated (redundant) features. The wrapper-based techniques incorporate supervised learning algorithms to evaluate the possible interactions between the features. These techniques add features to the selected subset iteratively and evaluate the subset by combining it with the trained decoder (Lotte et al., 2018). These methods are effective but require a great deal of computing time. On the other hand, embedded techniques integrate the feature selection process directly into the decoding algorithm, thus combining the advantages of both the filter-based and wrapper-based methods (Khaire and Dhanalakshmi, 2019). Embedded feature selection is a promising approach as it is directly performed during the model learning process.

Regularization/penalization is generally performed in a single-wise manner. Features are regularized independently and are not evaluated as belonging to a group of features. However, there are many applications, particularly in BCI, with structurally grouped input features. This approach allows the simultaneous setting of zero of the model coefficients within a group, which can be beneficial in cases where the input features are structurally grouped, such as excluding sensors (Eliseyev et al., 2017). For such applications, a single-wise sparsity-promoting penalty may be suboptimal. Group-wise regularization algorithms perform the feature selection process by grouping and applying penalization to the groups at once (Martínez-Montes et al., 2008; Eliseyev et al., 2012; Giordani and Rocci, 2013; Zhang et al., 2013; Hastie et al., 2015). This grouping can cluster features across various modalities, such as electrodes and frequency bands (van Gerven et al., 2009). Group-wise sparsity is ideal for naturally structured data, allowing for the elimination of variables (such as electrodes or frequency) from the signal processing workflow and reducing computational costs. Additionally, it may simplify the model's interpretation.

Despite being potentially beneficial for studies on BCI, group-wise penalization has rarely been applied in this field. The natural structure of BCI features space grouping features over modalities, such as electrodes, frequency bands, and time delay, which are heavily exploited (van Gerven et al., 2009; Eliseyev et al., 2012; Motrenko and Strijov, 2018; Wu et al., 2019). For the penalization of such groups, data may be viewed in the form of a tensor (Hastie et al., 2015; Eliseyev and Aksenova, 2016). Tensors, or multiway arrays, are higher-order generalizations of vectors and matrices (for more details, see Ref. Hsu et al., 2016). Tensors have several dimensions, also known as ways of analysis or modes. In BCIs, recorded signals are primarily analyzed in spatial, frequency, and temporal domains. To extract neural features, each epoch of neural signal recording is commonly mapped to temporal-frequency-spatial (or frequency-spatial) space (e.g., Chao et al., 2010; Schaeffer and Aksenova, 2016; Eliseyev et al., 2017; Choi et al., 2018). Tensor-based analyses, which present features in matrix form using tensor unfolding, are reported to be beneficial for BCI decoding as they preserve the natural structure of data (Zhao et al., 2013; Cichocki et al., 2015; Zhang et al., 2015; Eliseyev et al., 2017).

The tensors of data in BCI often benefit from tensor decomposition techniques. Slice-wise data representation can be obtained through sparsity-promoting penalization of tensor decomposition. Regularized PARAFAC and Tucker decomposition are algorithms designed for group-wise tensor penalization. Slice-wise sparsity-promoting penalization is added to the N-way Partial Least Square (NPLS) in van Gerven et al. (2009), Eliseyev et al. (2012), Motrenko and Strijov (2018), and Wu et al. (2019). These techniques have been used in a few offline BCI studies (Martínez-Montes et al., 2008; Eliseyev et al., 2012) and in other fields (Giordani and Rocci, 2013; Kim et al., 2013, 2014; Hervás et al., 2019). In Eliseyev et al. (2012), the L1-Regularized N-PLS algorithm was shown to be superior to its non-penalized version by suppressing noisy/irrelevant electrodes. However, these algorithms were not adapted for online decoder training in closed-loop use.

Most of the presented feature-dimensional reduction algorithms were tested offline. Additionally, feature selection performed in an offline preliminary study (Brunner et al., 2006; Huang et al., 2009; Spüler et al., 2012; Marathe and Taylor, 2013; Bousseta et al., 2018; Cantillo-Negrete et al., 2018; Kim et al., 2018; Nagel and Spüler, 2019) may not be optimal when using a CLDA strategy (Schlögl et al., 2010; Clerc et al., 2016). CLDA involves training the decoder online on data acquired during closed-loop BCI control sessions. Decoders trained in this manner have been reported to outperform decoders trained offline using data from open-loop BCI experiments (Jarosiewicz et al., 2013). Adaptive/incremental learning algorithms are particularly suited for the CLDA strategy. These algorithms continuously update the model, using only the latest data block and relevant statistics on the older signals, while not retaining the whole signals in memory (Schlögl et al., 2010; Brandman et al., 2018; Lotte et al., 2018). Adaptive/incremental learning algorithms are beneficial for the CLDA training strategy as they allow higher decoder update rates and are compatible with long decoder learning periods, which are generally necessary for high-dimensional control. However, only a few adaptive dimensional reduction algorithms were proposed. The adaptive dimensional reduction algorithms applied in the BCI (Zhao et al., 2008; Ang et al., 2011; Song and Yoon, 2015; Woehrle et al., 2015; Hsu et al., 2016; Mobaien and Boostani, 2016; Sannelli et al., 2016; Chen and Fang, 2017; Lotte et al., 2018) and other (Dagher, 2010) fields are primarily based on projection strategies (adaptive CSP, PCA, ICA, and xDAWN algorithms) and were only tested offline. Similarly, a few adaptive feature selection algorithms were proposed. Filter methods were tested on BCI simulations using mutual information (Oliver et al., 2013) or during online BCI experiments based on the Fisher score (Faller et al., 2012). The wapper-based strategy was optimized using parallel computation for the online BCI classifier (Mend and Kullmann, 2012), whereas embedded methods using semi-supervised feature selection (Long et al., 2011) and a weighting feature algorithm (Andreu-Perez et al., 2018) have been designed and used during online BCI applications. An adaptive genetic algorithm was proposed for adaptive channel selection in Moro et al. (2017). However, these methods have been applied to simple online binary classification only.

In the study, an adaptive algorithm promoting group-wise model sparsity, Lp-penalized recursive exponentially weighted N-way Partial Least Square (PREW-NPLS), is proposed. Lp, p = 0, 0.5, 1 norm/pseudo-norm penalty is applied to feature groups corresponding to the slices of data tensor related to the mode of analysis (e.g., spatial, frequency, temporal). The algorithm was tested with data recorded during BCI sessions of left/right arm 3D translations of a virtual avatar by a tetraplegic patient during the clinical trial NCT02550522 (ClinicalTrials.gov) conducted at the Grenoble Alpes University Hospital (CHUGA) (Benabid et al., 2019; Moly et al., 2022). The datasets were recorded during online closed-loop experiments using the REW-NPLS decoder that was previously integrated into the BCI system. To reproduce online experimental conditions and to ensure that the proposed algorithm is compatible with real-time CLDA, a pseudo-online simulation was conducted using the same parameters (buffer size, batch training, and so on) and the same model application procedure as it was used in real time. The comparison study was restricted to penalization according to spatial modality, which is the most critical at BCIs due to data recording/transfer limitations. For each type of penalization, a set of models/penalization parameters were evaluated. The PREW-NPLS decoders highlighted equivalent or better decoding performance compared to the generic REW-NPLS algorithm for the majority of the penalization parameters. The sparsest solutions allowed the removal of up to 75% of the electrodes without decreasing performance. L0-PREW-NPLS and L1-PREW-NPLS are sufficiently computationally efficient for online closed-loop decoder adaptation at a high-frequency rate.

Algorithms of the PLS family are widely used in BCI studies due to their stability in the case of high-dimensional data and in the presence of correlated and/or irrelevant variables. In motor BCI, the algorithms of the PLS family were applied in both continuous and discrete BCIs. Offline hand/fingers trajectory decoding (Chao et al., 2010; Chen et al., 2013; Eliseyev and Aksenova, 2014; Bundy et al., 2016; Schaeffer and Aksenova, 2016; Schaeffer, 2017; Choi et al., 2018), real-time hand translation/wrist rotation control (Benabid et al., 2019; Moly et al., 2022), error potential (ERP) detection (Rouanne et al., 2021) from ECoG, and EEG/MEG-based classification (Trejo et al., 2006; Eliseyev et al., 2017; Maleki et al., 2018) using the PLS algorithm were reported in preclinical and clinical studies. In Eliseyev and Aksenova (2014), GAM-PLS (generalized additive model—the partial least square) is reported to outperform the generic PLS in the presence of artifacts for 3D hand trajectory decoding from ECoG data in non-human primates (NHP). In Schaeffer (2017), PLS outperformed principle component regression and demonstrated comparable results with Lasso regression for hand trajectory and 1D finger trajectory decoding from ECoG in preclinical and clinical experiments. In Rouanne et al. (2021), NPLS demonstrated comparable results with Logistic Regression, SVM, MLP, and CNN for ERP detection from ECoG recordings in the sensory-motor cortex of tetraplegics, etc.

The generic PLS is a linear regression algorithm based on the iterative projection of input and output variables into the latent variable spaces of dimension f. The hyperparameter f is generally estimated through cross-validation in the preliminary study. Projectors are set to maximize the covariance between the input and latent output variables. The generic PLS is an offline algorithm. For online data stream modeling, recursive PLS (RPLS) and recursive exponentially weighted PLS (REW PLS) (Helland et al., 1992; Dayal and MacGregor, 1997; Qin, 1998) were developed. All the aforementioned PLS algorithms are vector-input-vector-output algorithms. N-way Partial Least Square (NPLS) is a generalization of the conventional PLS for tensor data (Bro, 1996, 1998). The NPLS algorithm projects the input and output tensors into a low-dimensional space of latent variables using a low-rank tensor decomposition. The recursive N-way PLS (RNPLS) (Helland et al., 1992) and recursive exponentially weighted N-way PLS (REW-NPLS) (Dayal and MacGregor, 1997) are generalizations of the adaptive RPLS and REW PLS algorithms to tensor variables and allow online tensor data stream learning of the regression model. RNPLS still requires fixing the hyperparameter f from the offline preliminary study, whereas REW-NPLS proposes a recursive validation procedure for the online optimization of the hyperparameter, enabling a fully adaptive algorithm (Dayal and MacGregor, 1997). In addition, the REW-NPLS algorithm is more computationally effective than the RPLS algorithm (Dayal and MacGregor, 1997).

The adaptive REW-NPLS algorithm has been tested offline in BCI studies for trajectory decoding from ECoG signals and for classification from MEG data, demonstrating similar or better results compared to other algorithms from the PLS family designed for offline use (Eliseyev et al., 2017). It was applied in real time for closed-loop adaptation of 3D hand translation/wrist rotation decoders in tetraplegics (Benabid et al., 2019; Moly et al., 2022). Finally, REW-NPLS was tested in the simulation of auto-adaptive continuous (bi-directional cursor control) and discrete multiclass motor imagery (MI) BCI in a tetraplegic patient (Rouanne et al., 2022). In the (Sliwowski et al., 2022) offline study, ANN algorithms were reported to outperform the REW-NPLS decoder. However, these algorithms cannot be applied in real time under CLDA.

Fully adaptive REW-NPLS is compatible with CLDA. However, it may be further improved by integrating real-time adaptive dimension reduction and promoting group-wise decoder sparsity using regularization. Group-wise sparsity, e.g., in the spatial dimension, may allow the elimination of irrelevant or highly correlated electrodes, decreasing computational time and memory consumption at the BCI use stage. This may be critical for portable BCI systems. A sparse solution may be advantageous for small training data sets, preventing overtraining. As BCI decoders require regular updates due to neuronal signal non-stationarity, reducing the decoder training time is desirable for real-life scenarios.

NPLS (Bro, 1996, 1998) estimates a linear relationship between a tensor of independent (input) and a tensor of dependent (output) variables. Given Xt∈ ℝI1×…×Im and Yt∈ℝJ1×…×Jn the m and n order tensors of the input and output variables at time t, Yt = Betat+bias+Dt, where Beta and bias are the tensors of parameters and their associated bias, D_t∈ℝJ1×…×Jn is the tensor of noise. The parameters are estimated from the training dataset {X_,Y_}, X_∈ℝL×I1×…×Im, Y_∈ℝL×J1×…×Jn, L is the training dataset size. NPLS constructs the linear regression iteratively by projecting tensors of observation X_ and Y_ to the space of latent variables using tensor decomposition: X_=∑fi=1frfi∘ wfi1 ∘…∘ wfim+E_x, Y_=∑fi=1fufi∘ qfi1∘…∘ qfin+E_y. Here, “∘” is the outer product operator (Eliseyev et al., 2017), Rf=[r1,…,rf]∈ℝL×f and Uf=[u1,…,uf]∈ℝL×f are matrices of the latent variables, Wfi=i=1m, wfii∈ℝIi and Qfi=i=1n, qfii∈ℝJi, fi = 1, …, f are the projection vectors constructed at iteration fi, and Ex, Ey are the residual tensors. The set of projectors is constructed iteratively, increasing latent variable space dimensions.

REW-NPLS (Eliseyev et al., 2017) update a set of F (F ∈ ℕ* is the fixed upper bound latent space dimension) models f=1F using current blocks of tensors of observation and previously computed models. Here, UI ∈ ℕ is the updated iteration number,  ​​ ​Beta_UIf∈ℝ(I1×…×IM)×(J1×…×JN),  ​​ ​​ bias_UIf∈ℝJ1×…×JN are current update of model coefficients, and X_UI∈ℝΔL×I1×…×IM, Y_UI∈ℝΔL×J1×…×JN are the current input and output tensors of observations. ΔL ∈ * is the number of samples recorded between the update blocks UI − 1 and UI.

The REW-NPLS algorithm is a generalization of the REW-PLS algorithm to tensor variables and belongs to the family of kernel PLS algorithms (Eliseyev et al., 2017). It evaluates a set of projectors and model coefficients from covariance tensors XY_=X_ ×1Y_, XY_∈ℝ(I1×…×Im)×(J1×…×Jn), and XX_=X_×1X_, XX_∈ℝ(I1×…×Im)×(I1×…×Im). Here, “ × k” is the k-mode tensor product (Eliseyev et al., 2017). First, a set of input variable projectors W are evaluated from the covariance tensor XY. The projectors are estimated using a rank one decomposition of the tensor V_= XY_∈ℝI1×…×Im in the case of single output. For higher dimensions, the eigenvector with the largest eigenvalue is computed from the covariance tensor XY to decrease the dimension of the tensor to decompose: e = eig(XYTXY), V_=reshape (XY·e),V_∈I1×…×Im. Here, XY∈ℝ(I1·…·Im)×(J1·…·Jn) is the unfolded tensor XY, and e = eig() is an eigenvector with the largest eigenvalue. Output projectors and the model parameters Beta and bias are computed using W and covariance tensors XX and XY. The projectors sets and model parameters are evaluated sequentially, increasing the latent variables' space dimension in the internal REW-NPLS iterations f.

Finally, at the UIth update, the covariance tensors XXUI and XYUI are computed from the previous XXUI−1 and XYUI−1 tensors, and the current block of observations :

μ1 is a forgetting factor. A set of projectors is evaluated using a rank-one tensor decomposition (Eliseyev et al., 2017) of the current tensor VUI.

In the REW-NPLS algorithm, only the covariance tensors, the normalization coefficients, and the current model are stored together with the current block of observations collected since the previous update.

Several tensor decomposition strategies were designed: the Parallel factor analysis (PARAFAC), Tucker, multilinear SVD decomposition, and so forth (Cichocki et al., 2015). Similar to generic NPLS, the REW-NPLS algorithm employs PARAFAC tensor decomposition (Eliseyev et al., 2017). PARAFAC or CANDECOMP/PARAFAC (CP), also known as polyadic decomposition (PD), can be considered the generalization of principal component analysis (PCA) and singular value decomposition (SVD) to the tensor case (Sheikhattar et al., 2015; Sharghian et al., 2019). This method represents a M-order tensor V_∈ℝI1×…×Im as the linear combination of vectors' outer products (rank-one tensors) such as follows:

Here, 1 ≤ i ≤ m corresponds to the ith mode/dimension of the tensor variable, “∘” is the (vector) outer product of the decomposition factors (projectors) wri∈ℝIi, R ∈ ℕ is the number of rank-one tensors used for decomposition, ρr is the weight associated with each rank-one tensor of the decomposition and E_∈ℝI1×…×Im is the residual tensor (Kolda and Bader, 2009). PARAFAC evaluates the projectors, minimizing the residuals.

Similar to generic NPLS, only one step of PARAFAC (R = 1) is applied to the current tensor decomposition at each internal iteration fi = 1, …, f of REW-NPLS:

To solve the optimization problem, the alternating least squares (ALS) algorithm is employed in REW-NPLS. ALS optimizes one projector iteratively at a time and fixes others, reducing, at each iteration of ALS, the optimization problem to a least-squares linear regression (Kolda and Bader, 2009; Cichocki et al., 2015; Pereira Da Silva et al., 2015).

REW-NPLS includes several iterations inside other iterative procedures: ALS iterations for PARAFAC (R = 1) decomposition, internal REW-NPLS iterations increasing latent variables space dimension, and, finally, the update iterations UI.

Sparse input variables projectors may result in a sparse model. Variable excluded from all the projectors Wfi, fi = 1, …, f is excluded from consideration. As for the tensor data, the projection is made according to the mode of analysis (e.g., special, frequency, or temporal). The sparsity of the projectors may allow excluding the slices of data from the model (e.g., exclude non-informative or redundant electrodes). To achieve sparse projectors, the proposed PREW-NPLS algorithm employed a penalized version of the PARAFAC (R = 1), introducing sparsity-promoting penalization. Lp, p = 0, 0.5, 1 penalization, being the classic lasso regularization (L1) or less conventional L0 and L0.5 penalization is studied. This section describes Lp-Penalization in PARAFAC (R = 1) and its integration into the REW-NPLS algorithm to build a new sparsity-promoting adaptive PREW-NPLS algorithm.

To simplify the notations without losing generality, a case of a three-order tensor is considered in the section. Given a three-order tensor, one step of PARAFAC (R = 1) applied at each iteration of REW-NPLS solves the optimization problem:

V_∈ℝI1×I2×I3 and wi∈ℝ*Ii, i = 1, 2, 3. As the norms of the projectors are arbitrary values in (1), the decomposition vectors are evaluated, minimizing V_−w1 ∘ w2 ∘ w32 before being normalized. The optimization problem (1) is solved using the ALS procedure. At each step, ALS fixes two of the three vectors w1, w2, w3 reducing the problem to a linear least-squares optimization, e.g. w2, w3 are fixes to approximate w1 and then w2 is evaluated fixing w1, w3 etc. until convergence (Uschmajew, 2015):

Here, V_(i)=(v11…v1I1)∈RI1×I2I3 is the tensor V_ unfolded according to i-the direction and ⊗ is the Kronecker product. Taking into account that (w2⊗w1)T∈RI1I2, (w3⊗w1)T∈RI1I3 and (w3⊗w2)T∈RI2I3 are vectors, the optimization tasks are separated into element-wise optimizations:

Here, w1=(w11,…,  wI11)T∈ℝ*I1, w2=(w12,…,  wI22)T∈ℝ*I2, and w3=(w13,…,  wI33)T∈ℝ*I3. The least square (LS) solutions of (2)–(4) are:

In this study, sparse promoting penalization using Lp, p = 0, 0.5, 1 norm/pseudo norm is proposed to be integrated into the cost function of PARAFAC (R = 1) to provide a slice-wise sparsity to the solution. The optimization task (1) is replaced by:

0 < λi ≤ 1 are regularization coefficients, the Kronecker delta δ0,wki=1 if wki=0, and δ0,wki=0 otherwise. Notably, only a part of the indices defined by a set ℒi ⊂, i = 1, 2, 3 are penalized, resulting in selective penalization (Lutay and Khusainov, 2021). A set ℒi  may vary depending on the REW-NPLS iteration. Selective penalization is introduced for the integration of Lp. -Penalized PARAFAC (R = 1) into the REW-NPLS algorithm.

Similar sparse promoting penalization in PARAFAC (R = 1) was considered in Eliseyev et al. (2012). However, this study was limited to L1-norm penalization. Integrated with the conventional non-adaptive NPLS algorithm, the optimization problem was solved using time consuming numerical optimization, which was then applied offline. In the current manuscript, a more general case of Lp, p = 0, 0.5, 1 norm/pseudo norm penalization is considered. In addition, selective penalization is applied for efficient integration of penalized PARAFAC into the REW-NPLS algorithm. Finally, an efficient optimization procedure compatible with real-time online applications is proposed.

The same ALS strategy is applied to complete the optimization task (8). ALS fixed all projectors at each step except one, leading to three successive optimization tasks. The least square (LS) sotion wLSi of the non-regularized problem is used as an initial approximation. Notably, unlike in non-regularized optimization, due to the penalization terms, the norms of the projectors are not arbitrary values anymore. Therefore, the normalization of a current estimate is added to ALS optimization iterations:

It should be noted that all considered regularization functions are decomposed as a sum of element-wise functions. Consequently, similarly to (3)–(5) optimization tasks, (9)–(11) are split into element-wise optimizations:

gp(wji)=f=1F. For f = 1, all projector elements can be potentially penalized: ℒj,1=, j = 1, …m. After the first set of projectors' extraction, non-zero elements of the projectors correspond to tensor slices already included in the decoding model with non-zero coefficients. For the next iterations, corresponding indexes were removed from a set of indexes to be penalized, resulting in a sequence ℒj,2⊂ℒj,1. A scheme representing the PREW-NPLS algorithm in the case of penalizing one data tensor direction is shown in Figure 1. www.frontiersin.org

Figure 1. Penalized REW-NPLS (PREW-NPLS) algorithm. For the model update iteration UI, an updated training data set is collected. The update loop includes the recursive validation procedure for the evaluation of optimal latent variable space dimension f (number of factors) and penalized NPLS loop, which includes the sequential tensors XY and XX decomposition. Each step of decomposition uses Rank 1 PARAFAC, which employs penalized alternating least square procedure. Selective sparsity-promoting penalization with penalization parameter λ is applied, e.g., to the spatial dimension. Part of the weights is set to zero due to penalization (shown in red). This results in zero slices in the decoding model (shown in red). A set of updated models for all parameters f are stored up to the maximum f value, 100 in the current setting.

2.5. Experiments

This study relies on the neural signal dataset recorded during the online closed-loop BCI clinical experiments. The ≪ BCI and Tetraplegia ≫ clinical trial (NCT02550522, ClinicalTrials.gov) (University Hospital, Grenoble, 2015) was approved by French authorities: National Agency for the Safety of Medicines and Health Products (Agence nationale de sécurité du médicament et des produits de santé, ANSM) with the registration number 2015-A00650-49 and the Ethic Committee for the Protection of Individuals (Comité de Protection des Personnes, CPP) with the registration number 15-CHUG-19. The research activities were carried out in accordance with the guidelines and regulations of the ANSM and the CPP. The patient provided informed consent for the clinical trial and publication as well as to publish the information/image(s) in an online open-access publication. Details of the clinical trial protocol are available in Benabid et al. (2019).

The participant was a 29-year-old right-handed male with traumatic sensorimotor tetraplegia caused by a complete C4–C5 spinal cord injury two years prior to the study. He underwent bilateral implantation of two chronic wireless WIMAGINE implants (Benabid et al., 2019) for ECoG signal recording on 21 June 2017. Two WIMAGINE recording systems were surgically implanted into the skull near the sensory-motor cortex (SMC) through a 25-mm radial craniotomy. Before the surgery, the patients' SMC was clearly localized using functional imaging. Details are provided in Benabid et al. (2019). The WIMAGINE device is made up of an active implantable medical device composed of 64-plane platinum-iridium 90/10 electrodes with a 2.3 mm diameter and a 4–4.5 mm inter-electrode distance (Sauter-Starace et al., 2019). The recorded signals were low- and high-pass filtered, with a bandwidth range of 0.5Hz to 300Hz, using analog low-pass filters as well as a digital low-pass FIR filter embedded into the implant hardware. The digitized ECoG data from 32 electrodes from each implant (Figure 2) were radio transmitted to a custom-designed base station at a 586 Hz sampling rate.

www.frontiersin.org

Figure 2. Experimental design. (A) The patient underwent bilateral implantation of two chronic wireless WIMAGINE implants. The recording systems were epidurally implanted within a 25 mm radial craniotomy in front of the sensory-motor cortex. The WIMAGINE implant is composed of 64 plane electrodes with a 2.3 mm diameter and a 4–4.5 mm inter-electrode distance. The digitized ECoG data from 32 electrodes from each implant (shown in red) are radio transmitted to a base station at the 586 Hz sampling rate. Reference electrodes are shown in green. (B) The database was recorded during the online closed-loop BCI experiments of alternative upper limbs BCI control. The virtual avatar was used as visual feedback. (C) The data sets were recorded in the period of 468 to 666 days after implantation and includes 43 sessions in general. The first six experimental sessions were used for incremental real-time decoder training. Separate decoders were identified to control different tasks: alternative 3D reaching tasks, and wrists rotation. They were mixed together using the Markov Mixture of expert approach (Moly et al., 2022). The decoders were tested without re-calibration with 37 experimental sessions distributed over 5–203 days after the last model update. Only 3D reaching tasks of both hands are included in the current study. The 3D reaching tasks are composed of a series of trials in which the sequentially proposed targets must be reached (pursuit tasks). The cursor position is not reset during the task, between tasks, and during the idle state.

Since the implantation, the patient had been trained to control multiple real and virtual effectors, such as games that were specially created for BCI training, a wheelchair, an exoskeleton, an avatar, and so on, using a custom-made BCI platform (Benabid et al., 2019). The database used in the study was recorded during the online closed-loop BCI experiments of upper limbs BCI control. The BCI sessions included alternative active states (AS) of the 3D reaching task for each hand, 1D wrist rotation of each hand, and the idle state (IS) (Benabid et al., 2019). The patient aimed to reach the proposed targets or rotate the wrist to specific angles following pursuit tasks. A pursuit task session was composed of successive tasks, e.g., a left-hand 3D reaching task, a right-hand 3D reaching task, IS, and so on. Each task is composed of several trials in which the cursor must reach the proposed targets. The cursor position is not reset between tasks, during the task, or during the idle state. Twenty-two targets were symmetrically distributed in two cubes in front of the patient. The virtual avatar was used as visual feedback (Benabid et al., 2019). In this study, only the left- and right-hand 3D translation trials were used for algorithm evaluation. The experimental paradigm is illustrated in Figure 2.

Data were recorded in the period of 468–666 days after implantation and included 43 sessions. The average session duration was 29 ± 8 min. During the experiments, the first six experimental sessions were used for incremental real-time decoder updates. REW NPLS decoders for each active state task were integrated into the BCI clinical trial platform. The decoders were initialized with a zero. The total training time of the decoding models was 3 h and 37 min, including a total of 189 and 194 trials of the left- and right-hand translation, respectively. The decoders were calibrated during experimental sessions in late September 2018 and were tested without re-calibration in experiments from early October 2018 to mid-March 2019, with 37 experimental sessions distributed over 5–203 days after the last model update session (468–666 days after implantation; Figure 2).

2.6. Application, feature extraction

During the experimental sessions, neural signal epochs of 1s for 64 recording channels Xt∈ℝ586x64 with sliding steps of 100 ms used for feature extraction. The epochs were mapped into the feature space using a complex continuous wavelet transform (CCWT) (Morlet) with a frequency range of 10–150 Hz and a 10-Hz step. CCWT is a common feature extraction strategy that is widely used in the field of BCIs (Chao et al., 2010; Shimoda et al., 2012; Schaeffer and Aksenova, 2016; Eliseyev et al., 2017; Choi et al., 2018). The absolute value of CCWT was decimated by averaging along the temporal modality to obtain a 10-point description of 1s time epoch for each frequency band and for each channel, resulting in the temporal-frequency-spatial neural feature tensor X_t∈ℝ10x15x64. The observations of neural features Xt and the movement features yt recorded during the decoder update/calibration sessions were used for the decoding model identification. The movement features (optimal prediction) at time t is defined as the 3D Cartesian vector between the current effector position and the target position (Benabid et al., 2019). The control command ŷt for 3D hand translation is defined as the cartesian increments of the current position and is sent to the effector at 10 Hz.

2.7. Comparison study

To compare the proposed PREW-NPLS and the generic, non-penalized REW-NPLS algorithms, the decoding of the 3D reaching task for each hand is considered. The comparison study was restricted to penalization according to spatial modality. Three types of penalization Lp, p = 0, 0.5, 1 were tested. For Lp, p = 0.5, 1 penalties, a set of models were evaluated with the penalization parameter λ going from 0 to 0.5 with 0.02 steps. In the case of the L0-PREW-NPLS, preliminary results highlighted that the studied λ range was not relevant. Therefore, the models with the penalization parameter λ going from 0 to 0.05 with 0.002 steps were estimated. The cases with λ= 0 correspond to the generic, non-penalized REW-NPLS algorithm.

The performance of algorithms was evaluated offline in a pseudo-online manner. Pseudo-online simulation uses the procedure and parameters of online stream data processing. The dataset was recorded during the online closed-loop BCI experiments using the REW-NPLS decoder previously integrated into the BCI system. As the recorded data integrate the neuronal feedback into real-time decoding, the presented offline comparison study is not fully generalizable to the online case. Nevertheless, it allows the characterization of, to some extent, the studied algorithms.

To be as close as possible to the settings of the online experiments, in the comparison study, the penalized models were calibrated on the same experiments that were used for decoder training during the online closed-loop experiments.

The predicted trajectories performed during the online closed-loop experiments are related to the decoding model used during the experiments. Consequently, sample-based indicators were used to compare the predictions of the tested algorithms in the offline study. The cosine similarity indicator was based on the comparison between the predicted directions ŷt and the optimal direction yt was employed:

CosSimt=yt·y^t‖yt‖‖y^t‖.

Here, “·” defined the dot product yt and ŷt, and they are the optimal and predicted output. The optimal 0 output is defined as the 3D Cartesian vector between the current position and the target. Notably, CosSimt[−1, 1] represent how direct the movement is to the final destination. A mean cosine similarity of 1 corresponds to a direct and short trajectory. CosSim performance criterion is evaluated as the median, 25th (Q1) and 75th (Q3) percentiles of the CosSimt over samples and is notated as CosSim = median (Q1−Q3).

The PREW-NPLS algorithms converge to sparse solutions, fixing non-relevant model coefficients to exactly zero. The decoding performance, therefore, is not the only relevant indicator. Considering a penalized model with the penalization restricted to the ith dimension, the model Sp

留言 (0)

沒有登入
gif