<figure><img src="image_1.png" /></figure>

HAL
open science

# An adaptive closed-loop ECoG decoder for long-term and stable bimanual control of an exoskeleton by a tetraplegic

Alexandre Moly, Thomas Costecalde, Félix Martel, Matthieu Martin, Christelle Larzabal, Serpil Karakas, Alexandre Verney, Guillaume Charvet, Stephan Chabardes, Alim Louis Benabid, et al.

<figure><img src="image_2.png" /></figure>

## ► To cite this version:

Alexandre Moly, Thomas Costecalde, Félix Martel, Matthieu Martin, Christelle Larzabal, et al. An adaptive closed-loop ECoG decoder for long-term and stable bimanual control of an exoskeleton by a tetraplegic. Journal of Neural Engineering, 2022, 19 (2), pp.026021. 10.1088/1741-2552/ac59a0.cea-04132467

HAL Id: cea-04132467

https://cea.hal.science/cea-04132467v1

Submitted on 23 Sep 2024

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_3.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_4.png" /></figure>

CrossMark

OPEN ACCESS

RECEIVED

30 August 2021

REVISED

18 February 2022

ACCEPTED FOR PUBLICATION

1 March 2022

PUBLISHED

30 March 2022

Original content from
this work may be used
under the terms of the
Creative Commons
Attribution 4.0 licence.

Any further distribution
of this work must
maintain attribution to
the author(s) and the title
of the work, journal
citation and DOI.

<figure><img src="image_5.png" /></figure>

PAPER

# An adaptive closed-loop ECoG decoder for long-term and stable bimanual control of an exoskeleton by a tetraplegic

Alexandre Moly¹ ID, Thomas Costecalde¹ ID, Félix Martel¹,*, ID, Matthieu Martin¹ ID, Christelle Larzabal¹ ID, Serpil Karakas¹, Alexandre Verney² ID, Guillaume Charvet¹ ID, Stephan Chabardes³ ID, Alim Louis Benabid¹ ID and Tetiana Aksenova¹,*, ID

¹ University Grenoble Alpes, CEA, LETI, Climatec, Grenoble, F-38000, France

² Université Paris-Saclay, CEA, List, F-91120 Palaiseau, France

³ Centre Hospitalier Universitaire Grenoble Alpes, 38700 La Tronche, France

* Authors to whom any correspondence should be addressed.

E-mail: felix.martel@cea.fr and tetiana.aksenova@cea.fr

Keywords: BCI, brain-computer interface, mixture of experts, adaptive, tetraplegic, exoskeleton control, ECoG

## Abstract

Objective. The article aims at addressing 2 challenges to step motor brain-computer interface (BCI) out of laboratories: asynchronous control of complex bimanual effectors with large numbers of degrees of freedom, using chronic and safe recorders, and the decoding performance stability over time without frequent decoder recalibration. Approach. Closed-loop adaptive/incremental decoder training is one strategy to create a model stable over time. Adaptive decoders update their parameters with new incoming data, optimizing the model parameters in real time. It allows cross-session training with multiple recording conditions during closed loop BCI experiments. In the article, an adaptive tensor-based recursive exponentially weighted Markov-switching multi-linear model (REW-MSLM) decoder is proposed. REW-MSLM uses a mixture of expert (ME) architecture, mixing or switching independent decoders (experts) according to the probability estimated by a 'gating' model. A Hidden Markov model approach is employed as gating model to improve the decoding robustness and to provide strong idle state support. The ME architecture fits the multi-limb paradigm associating an expert to a particular limb or action. Main results. Asynchronous control of an exoskeleton by a tetraplegic patient using a chronically implanted epidural electrocorticography (EpiCoG) recorder is reported. The stable over a period of six months (without decoder recalibration) eight-dimensional alternative bimanual control of the exoskeleton and its virtual avatar is demonstrated. Significance. Based on the long-term (>36 months) chronic bilateral EpiCoG recordings in a tetraplegic (ClinicalTrials.gov, NCT02550522), we addressed the poorly explored field of asynchronous bimanual BCI. The new decoder was designed to meet to several challenges: the high-dimensional control of a complex effector in experiments closer to real-world behavior (point-to-point pursuit versus conventional center-out tasks), with the ability of the BCI system to act as a stand-alone device switching between idle and control states, and a stable performance over a long period of time without decoder recalibration.

## 1. Introduction

Brain-computer interfaces (BCIs) create a new com-
munication pathway between the brain and an
effector without neuromuscular activation. Among
the various potential applications, the functional
compensation/restoration of individuals suffering
from severe motor disabilities has always been a focus

for BCI research. Major milestones have been reached
by the motor-BCI community over the years [1–3].
Nevertheless, many aspects need to be addressed to
translate BCI-driven systems from laboratories into
the patients’ home for daily life applications.

The primary challenge of motor BCIs in clinical
application is the high-dimensional control of effect-
ors using safe, biocompatible and chronic neural
recording systems. Brain signal recordings should remain stable and allow accurate neuronal signal decoding in conditions that are more demanding than in laboratories. The control of many degrees of freedom (DoF), up to ten, based on microelectrode array (MEA) recordings have been reported [1, 3]. However, MEA recording systems are highly invasive, have biocompatibility issues and poor stability. They suffer from a decrease of signal-to-noise ratio over time [4, 5], a high across-day variation in the neural signals [6, 7] and still require wired recording systems despite recent efforts in this domain. Electrocorticography (ECOG) provides a good compromise between invasiveness and signal resolution [8–10]. Numerous preclinical and clinical studies demonstrated the interest in ECOG-based BCIs to control effectors [11–22] and highlighted the good signal-to-noise ratio and the stability of ECOG signals over months and even years [23–26]. Clinical results of high-dimensional (up to 8D) alternative bimanual control of a complex effector by a tetraplegic subject using epidural ECOG (EpiCOG) arrays have been recently reported [19]. This study outperformed the previously reported state-of-the-art ECOG-based BCIs with up to 3D control [11, 12].

Another requirement to bring motor BCI into real life applications is to make the system act as a stand-alone device which can switch between idle/rest state (IS) and multiple active states (AS) [27, 28] (referred to as asynchronous BCI). This is a critical point as the majority of the reported BCIs are synchronous cue-based action-oriented systems providing neuronal control of effector to the user during specific time intervals defined by an operator. For example, continuous BCI performance is often evaluated for single limb center-out experiments which classically reset the cursor position between trials without including the rest period. Besides not being representative of real life applications, it is likely that this may lead to unwanted activations/movements when the user does not intend to control the effector [27].

Moreover, single limb applications are limited compared to daily life actions which commonly require synchronized or alternative bimanual (or generally multi-limb) movements. Despite its clear benefit for patient motor deficit compensation and rehabilitation, multi-limb decoding has been poorly explored in the BCI field. Most breakthroughs involved the control of a single robotic arm or the movement decoding of one hand. Bimanual experiments have only been tested using MEAs with virtual effectors [29] or during ECOG-based movement detection experiments [18] in non-human primates (NHPs). While bimanual and/or asynchronous BCIs are not very common, several decoding strategies have been proposed. In particular, the problem of multi-finger movement trajectory reconstruction from ECOG recordings was studied. In most of the

cases, hybrid models were employed by mixing classifier outputs to detect finger activations and continuous decoders to predict their respective movements [30–32].

A BCI decoder must be sufficiently optimized to enable computation time suitable for real-time application. Despite promising results, translating the off-line trajectory reconstruction algorithm to real-time closed-loop experiments is generally a challenging task. Even if a decoder meets the real-time requirements, drops in the decoding performance have been reported repeatedly when decoders calibrated off-line using open-loop experiments were used for online decoding [4, 33, 34]. Open-and-closed-loop model training lead to distinct decoders [35]. To take the patient feedback into account, BCI studies employ adaptive decoders which integrate the decoding model parameters identification into closed-loop BCI sessions. Adaptive decoders update their parameters in an incremental manner with new incoming data, optimizing the model parameters in real time. While several adaptive linear and nonlinear regression and classification decoders were proposed for MEA [7, 36–40] and electroencephalography (EEG) [41–45] driven BCI, only a few adaptive decoders were developed for ECOG recordings [17]. Most adaptive algorithms are restricted to linear decoders, which may be limiting for complex effector control with a high number of DoF. Closed-loop adaptive model calibration is one strategy to create model stable over time [7, 34] which is a major challenge, considering the non-stationarity and the intra-subject variability of the brain's signals (inattention, habituation). So far, stable long-term BCIs were only achieved using non-adaptive brain switch decoders (1D) for a period of four months with local field potentials [46] and 36 months with ECOG recordings [25].

In the present article, the recursive exponentially weighted Markov-switching multi-linear model (REW-MSLM) is proposed to address the lack of stable asynchronous BCI for bimanual effector control from EpiCOG recordings in tetraplegics [47]. The tensor-based piecewise linear REW-MSLM algorithm is an online adaptive supervised learning algorithm. It updates a decoder in real time in an incremental manner during the calibration period of closed-loop BCI experimental sessions. In this article we report the case study of a closed-loop 8D control performed by a tetraplegic patient in an exoskeleton and with a virtual avatar. These results outperform the 3D control of previous ECOG-based state-of-the-art BCIs [11, 12, 16]. We demonstrated the remarkable stability of the BCI system and the stable performance of the REW-MSLM decoder which was not recalibrated for more than five months when the patient was in the exoskeleton and for more than 6.5 months using the virtual avatar. For both effectors, the patient was able
to switch reliably between discrete states and demonstrated relevant control for continuous movements. The decoding performance outperformed the ECoG-based BCIs state-of-the-art for which such a long-term robustness was never reported before. Compared to the classic center-out experiments the patient was able to perform more complicated tasks such as multiple alternative point-to-point pursuit tasks.

## 2. Methods

### 2.1. REW-MSLM decoder

The REW-MSLM is an online tensor-based fully adaptive mixture of multi-linear experts algorithm (figure 1). The REW-MSLM inherits the MSLM [15] mixture of experts (MEs) structure, generalizing the MSLM model to tensor-input tensor-output variables and introducing the recursive model parameter identification procedure inspired by the N-way partial least squares (REW-NPLS) method [17].

#### 2.1.1. MSLM description

The MSLM [15] is a hybrid discrete/continuous decoder based on a ME model structure. A ME mixes or switches independent decoders, called ‘experts’. Basic assumption of ME is that each expert decodes its own specific region of feature space [48]. Experts are mixed according to the ‘gating’ model which estimates the probability of an expert to be activated or inhibited. This probability is used to compute gating coefficients to weight experts’ outputs. Additionally, MSLM uses dynamic gating assuming a hidden Markov model (HMM) for the state sequence to improve the decoder robustness. Conventional MSLM is a vector-input-vector-output model and employs linear experts. Both experts and gating models are identified offline. Application of MSLM in motor BCI was limited to offline studies: 3D-trajectory decoding of single limb wrist from NHPs ECoG with 2 states separating rest and movement periods, and 1D-trajectory decoding of fingers movement with states associated to individual finger activation and rest periods from ECoG recordings of abled body patients undergoing pre-surgical evaluation [15].

#### 2.1.2. REW-NPLS description

Due to the robustness in the computation of high dimensional data, algorithms of the partial least squares (PLS) family were frequently used in continuous and discrete BCI decoding experiments. Numerous publications which reported offline ECoG-based hand trajectory decoding [15, 16, 18, 20, 49, 50], and EEG-based classification or cursor decoding [51, 52] confirmed the interest for such algorithms. The classical PLS regression algorithm is an offline procedure based on the iterative projection of input $\mathbf{x}_t \in \mathbb{R}^m$ and output $\mathbf{y}_t \in \mathbb{R}^n$ variables into a latent variable space of

dimension $f$ (is referred as the PLS hyperparameter). Projectors are estimated by maximizing the covariance between the input and the output latent variables [53]. The subspace dimension $f$ is typically determined through cross-validation. For online modeling, recursive PLS (RPLS) and REW PLS [54–56] were proposed.

A generalization of the conventional PLS algorithm to tensor data: the NPLS algorithm, was proposed by Bro [57, 58]. A tensor is a generalization of a matrix to higher order dimensions, also known as ways or modes. Vectors and matrices are special cases of tensors with one and two modes respectively [59]. Tensor-based algorithms emerged as a promising strategy in the BCI field. They allowed simultaneous processing of high-dimensional data in the temporal, frequency and spatial domains [17, 59]. The NPLS algorithm projects the input and output tensors into low dimensional space of latent variables using a low rank tensor decomposition. It improves the stability and robustness of the model compared to the classic unfold PLS leading to more accurate and interpretable predictions [57, 58] while preserving the structure of the data.

For the online tensor data flow modeling, the recursive N-way PLS (RNPLS) [54] which is a generalization of the RPLS algorithm to tensor variables, and the REW-NPLS [17] inspired from kernel REW PLS [54–56] were proposed. RNPLS still requires fixing the hyperparameter $f$ from offline preliminary study whereas the Recursive-Validation procedure used in REW-NPLS for online optimization of the hyperparameter $f$ enables a fully adaptive algorithm [17]. The decoder is entirely tuned in real time. The Kernel REW-NPLS is also more computationally efficient than the RPLS algorithm.

#### 2.1.3. REW-MSLM description

The REW-MSLM inherits the MSLM MEs structure. A basic assumption of ME is that each expert decode its own specific region of feature space. Given $\mathbf{X}_t \in \mathcal{X} \subset \mathbb{R}^{I_1 \times \dots \times I_m}$ and $\mathbf{Y}_t \in \mathcal{Y} \subset \mathbb{R}^{J_1 \times \dots \times J_n}$ the independent and dependent $m$ and $n$ order tensor variables at time $t$ respectively, the feature space of independent variables is supposed to be partitioned into $K$ possibly overlapping regions $\mathcal{X} = \cup_{k=1}^K \mathcal{X}_k$. It is assumed that the space of input variables is mapped to the space of output variables using $K$ local multilinear functions $\Phi = \{\phi_k : \mathcal{X}_k \to \mathcal{Y}, k = 1, 2, \dots, K\}$. Let $z_t \in [1; K] \subset \mathbb{N}^*$ be a latent state variable which defines the selected local multilinear function at time $t$ (expert): $\mathbf{Y}_t = \phi_{z_t}(\mathbf{X}_t)$. Similar to MSLM, dynamic gating is introduced using a first-order HMM [15]. The latent state variable $z_t$ is assumed to follow the first order Markovian assumption, which states that the dependence of $z_t$ is limited to the past state $z_{t-1}$. $\mathbf{Y}_t$ is estimated as follows:
<figure><img src="image_6.png" /><figcaption>**Figure 1.** Recursive exponentially weighted Markov-switching linear model (REW-MSLM) architecture. The REW-MSLM includes a mixture of experts model, which can be described as the parallel computation of several predictions from different regression models (experts) that are weighted (enhanced or inhibited) according to the input variables using a classifier (gate). We hypothesize that the input feature space X can be divided into several specific local regions $X_k$ and that each sub-space can be fitted using local multilinear functions $\varphi_k$ associated with an expert. Multilinear functions $\varphi_k$ are estimated using k independent REW-NPLS models. The selected expert is determined based on the dynamic gating model. The gating model is a hidden Markov model (HMM) which computes the probability $\gamma_k$ for each expert to be activated. Commands are decoded by the REW-MSLM and sent to the effector to provide visual feedback to the patient.</figcaption></figure>

$$ \hat{\mathbf{Y}}_t = \sum_{k=1}^{K} \gamma_{k,t} (\mathbf{\text{Beta}}_k \mathbf{X}_t + \mathbf{\text{bias}}_k). $$

Here, $\mathbf{\text{Beta}}_k$ and $\mathbf{\text{bias}}_k$ are the $k^{th}$ expert's tensor parameters and its associated bias. $\gamma_{k,t} = p(z_t = k|\mathbf{X}_{1:t})$ is the dynamic gating weight coefficient associated with the $k^{th}$ expert at time t. REW-MSLM models are entirely defined through the experts' parameters $\theta_e = \{\mathbf{\text{Beta}}_k, \mathbf{\text{bias}}_k\}_{k=1}^K$ and HMM parameters $\theta_g = \{\mathbf{A}, \{d_k\}_{k=1}^K, \pi\}$, where $\mathbf{A}$ is the transition matrix, $\mathbf{A} = (a_{ij}) \in \mathbb{R}^{K \times K}$, $a_{ij} = p(z_t = j|z_{t-1} = i)$, $\{d_k\}_{k=1}^K$ is the set of parameters employed to estimate conditional emission probability of the observed variables $p(\mathbf{X}_t|z_t)$, and $\pi \in \mathbb{R}^K$ is the initial state probability vector at $t=0$.

## 2.2. REW-MSLM online/incremental training

The proposed REW-MSLM algorithm recursively estimates $\Theta = \{\theta_g, \theta_e\}$ with a supervised training procedure. At each update $u$, the corresponding block of training dataset $\{\mathbf{X}_u, \mathbf{Y}_u, \mathbf{z}_u\}$ is given with $\mathbf{X}_u \in \mathbb{R}^{\Delta L \times l_1 \times \dots \times l_m}, \mathbf{Y}_u \in \mathbb{R}^{\Delta L \times l_1 \times \dots \times l_n}, \mathbf{z}_u = (z_{t_1}, \dots, z_{t_1+\Delta L})^T \subset \mathbb{N}^{*\Delta L}$ and $\Delta L$ the update block size. The $K$ local multilinear functions $\varphi_k$ are estimated using expert's specific samples. The $k^{th}$ expert's parameter update is performed on the training dataset $\{\mathbf{X}_u^k, \mathbf{Y}_u^k\}$. $\mathbf{X}_u^k$ and $\mathbf{Y}_u^k$ are sub-tensors of $\mathbf{X}_u$ and $\mathbf{Y}_u$ formed by samples labeled as belonging to state $k$. The $k^{th}$ expert's parameters are updated using the REW-NPLS algorithm: REW-NPLS$_e$ = REW-NPLS($\mathbf{X}_u^k, \mathbf{Y}_u^k$) with the forgetting factor $\lambda_k, 0 \le \lambda_k \le 1$.

For online optimization latent variable space dimension (hyperparameter $f$), the REW-NPLS$_e$ algorithm estimates a set of $F$ models for each expert $\{\mathbf{\text{Beta}}_{u,k}^f, \mathbf{\text{bias}}_{u,k}^f\}_{k,f=1}^{K,F}$. $F \in \mathbb{N}^*$ is the fixed highest latent space dimension. The optimal hyper-parameter of the $k^{th}$ expert $f_k^* \le F$ is selected following the Recursive-Validation procedure [17]. For the currently available models, the Recursive-Validation exploits the newly available $\{\mathbf{X}_u^k, \mathbf{Y}_u^k\}$ dataset as testing data to evaluate the best hyperparameters before this dataset is used as a training dataset for the models updating. The best models are chosen independently for each expert: $\{\mathbf{\text{Beta}}_k, \mathbf{\text{bias}}_k\}_{k=1}^K = \{\mathbf{\text{Beta}}_{u,k}^f, \mathbf{\text{bias}}_{u,k}^f\}_{k=1}^K$, and are used for real-time decoding of the neural signals.

Similarly, at each update $u$, the HMM gating parameter are updated using the update block dataset $\{\mathbf{X}_u, \mathbf{z}_u\}$. The HMM transition matrix $\mathbf{A}$ is approximated by counting the successive transition of states in $\mathbf{z}_u$ and is weighted with the transition matrix estimated during the previous updates with the forgetting factor $\lambda_g, 0 \le \lambda_g \le 1$. The HMM conditional emission probabilities $p(\mathbf{X}_t|\mathbf{z}_t)$ are inferred through the combination of $p(z_t|\mathbf{X}_t)$ and class prior $p(z_t)$ using the Bayes' theorem [60]. The REW-NPLS discriminative decoder is embedded into the HMM-based gating process to evaluate $p(z_t|\mathbf{X}_{1:t})$. REW-NPLS was used because of its online adaptive characteristics and its relevance for high dimensional input variable. A discriminative decoder is selected instead of generative ones due to benefits for high dimensional and complex dependencies of features [61, 62]. The
decoder is trained on the observation tensor of input variables $\mathbf{X}_t$ and the latent state dummy variable matrix $\mathbf{Z}_t \in \{0, 1\}^{K \times \Delta L}$ where the column-wise (single) non-zero element depicts the activated state for each sample (one-hot encoding).

The discriminative REW-NPLS decoder computes a set of $F$ multilinear models $\{\mathbf{B}_t^f, \mathbf{b}_t^f\}_{f=1}^F$, where $\mathbf{B}_t^f \in \mathbb{R}^{K \times I_1 \times \dots \times I_m}$ and $\mathbf{b}_t^f \in \mathbb{R}^K$ are the tensor of the gating model parameters and its related bias. The Recursive-Validation procedure selects the best model based on the estimated gating optimal hyperparameter $f_g^* \le F$ and defines the optimal gating model as $\{\mathbf{B}, \mathbf{b}\} = \{\mathbf{B}_t^g, \mathbf{b}_t^g\}$ for the dynamic gating weight $\gamma_{k,t}$ estimation. The output variable $\hat{\mathbf{z}}_t \in \mathbb{R}^K$ determines how likely each hidden state is generated based on $\mathbf{X}_t$. The prediction $\hat{\mathbf{z}}_t$ is computed from the discriminative REW-NPLS decoder. Then, $p(\mathbf{z}_t|\mathbf{X}_t)$ is evaluated with the softmax function [48] to compute $\gamma_{k,t} = p(\mathbf{z}_t|\mathbf{X}_{1:t})$ using HMM forward algorithm.

REW-MSLM uses dynamic HMM gating. The equivalent ME algorithm using static gating (without HMM) is referred as REW-SLM. REW-SLM gating is computed with the REW-NPLS trained on explanatory variables and latent states, using the softmax function but without the HMM forward algorithm.

### 2.3. REW-MSLM application

In real time, each expert $\{\mathbf{Beta}_k, \mathbf{bias}_k\}_{k=1}^K$ output is estimated for each new input data buffer after feature extraction $\mathbf{X}_t$. The dynamic gating coefficients $\gamma_{k,t}$ are estimated using the latent state variable estimator $\hat{\mathbf{z}}_t$ post-processed with a softmax function [48] and the HMM forward algorithm [63]. The forward algorithm evaluates $\gamma_{k,t}$ by considering the past and current observations:

$$\hat{\mathbf{z}}_t = \mathbf{B}\mathbf{X}_t + \mathbf{b},$$

$$p(z_t = k|\mathbf{X}_t) = \frac{\exp(\hat{z}_{k,t})}{\sum_{i=1}^{K} \exp(\hat{z}_{i,t})},$$

$$p(z_t = k, \mathbf{X}_{1:t}) = p(\mathbf{X}_t|z_t = k) \sum_{j=1}^{K} a_{kj} \gamma_{k,t-1},$$

$$\gamma_{k,t} = p(z_t = k|\mathbf{X}_{1:t}) = \frac{p(z_t = k, \mathbf{X}_{1:t})}{\sum_{j=1}^{K} p(z_t = j, \mathbf{X}_{1:t})}.$$

### 2.4. Clinical trial description

The REW-MSLM algorithm was tested and applied as the neural signal decoder during the 'BCI and Tetraplegia' clinical trial (ClinicalTrials.gov, NCT02550522 [64, 65]). The clinical trial was approved by the 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), registration Number 2015-A00650-49, and the ethical Committee for the Protection of Individuals (Comité de Protection des Personnes—CPP), registration number 15-CHUG-19. All research activities were carried out in accordance with the guidelines and regulations of the ANSM and the CPP.

The REW-MSLM was tested on a single patient. The patient signed informed consent prior to surgery. Details of the clinical trial protocol are available in [19]. The subject 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. The patient can perform neck, shoulder and small upper limb movements by contraction of the biceps at the elbow and extensors of the wrists. American Spinal Injury Association Impairment (ASIA) scores the contraction of the biceps at the elbow at four and five for the right and left body side respectively, whereas extensors contractions were scored at 0 and 3 for the right and left wrists, respectively. With the exception of the cited muscles, all others muscles below were scored 0 on the ASIA scale. Moreover, the sensory-motor deficit was complete.

The patient underwent bilateral implantation of two chronic wireless WIMAGINE implants [19, 47] for EpiCoG signal recording. Two WIMAGINE recording systems were implanted into the skull within a 25 mm radius craniotomy placed in front of the sensory motor cortex (SMC) (figure 2). The electrodes located at the implant lower surface are in contact with the dura mater. Before surgery, the subject's SMC was localized clearly using functional imaging (fMRI and MEG) as the subject imagined virtual movements of his limbs or performed real motor tasks when possible. Details are provided in [19]. WIMAGINE is an active implantable medical device composed of 64 plane platinum iridium 90/10 electrodes with a 2.3 mm diameter and 4–4.5 mm inter-electrodes distance on the lateral and antero-posterior directions [24] dedicated to ECoG neural signal recordings. The WIMAGINE implant was shown to be safe for long-term EpiCoG signal recording [24, 47, 66]. The EpiECoG signals are low and high pass filtered in a bandwidth from 0.5 Hz to 300 Hz using analog bandpass filter and, after digitalization, using a digital low pass FIR filter with a cutoff frequency 292.8 Hz. The filters are directly embedded into the implant [47, 67]. Downsampled to 585.6 Hz, the data are radiotransmitted to a custom designed base station connected to a computer [47]. During the experimental sessions presented in the article, 32 electrodes for each implant were selected in a checkerboard-like pattern (figure 2) because of temporarily limited data rates caused by restricted radio link.

Since the implantation date, the patient was trained using a custom-made BCI platform to control multiple real and virtual effectors [19] (figure 2). The article presents a series of experiments performed
<figure><img src="image_7.png" /><figcaption>**Figure 2.** BCI platform for the ‘BCI and tetraplegia’ clinical research protocol in CLINATEC. Two wireless WIMAGINE implants with a 64-electrode array [47] are used to record EpiCoG signals. Implants were placed into the skull above the motor cortex by a craniotomy. 32 electrodes for each implant were selected in a checkerboard-like pattern for signal recording. BCI decoder translates the neural signals into decisions at 10 Hz frequency to control effectors. In parallel, the EpiCoG data and the movement instructions are used to update the model in real time during the closed-loop BCI calibration sessions. The model is updated at approximately 0.07 Hz (every 150 samples). The exoskeleton is used for training in CLINATEC and the virtual avatar is used for training at home. Both effectors provide a visual feedback to the patient to adapt and respond in a closed-loop fashion to model predictions.</figcaption></figure>

in the laboratory with the enhancing mobility (EMY) [19, 68] exoskeleton and with EMY’s virtual avatar replica for training at home. EMY is a wearable fully motorized four-limb robotic neuroprosthesis (14 joints, 14 actuated DoFs) equipped by a computer station receiving radio-emitted EpiCoG signals. The decoder translates the neuronal signal into the motor commands which activate the limbs and joints to produce movements, mimicking natural limbs movements. During the experiments, effectors were controlled at a 10 Hz frequency rate.

## 2.5. Experimental setup

The patient underwent bilateral implantation of two chronic wireless WIMAGINE implants [19, 47] on 21 June 2017 and since, underwent training for more than 28 months. Experiments were carried out in the laboratory three successive days per month. During these sessions, the patient was strapped into the EMY exoskeleton. For the remaining weeks, experiments were performed in the patient’s home three days a week. The patient was installed in his wheelchair in front of the computer screen (figure 2). During all the experiments, the patient was allowed to move and talk freely during the training and test sessions in order to create models that are robust to artefacts related to muscular activities such as head movements.

All the experiments were online closed-loop BCI experiments. Each experiment was divided into two phases. The training/calibration phase (optional) was designed for online updates of the decoder. Depending on protocol, the decoder was identified from

scratch (initialized by zero), or was updated from the decoder trained during previous session(s). During the test phase, the model was fixed. The test phase was used for the BCI performance evaluation.

A support/assistance system was optionally provided to the patient during the early model calibration phase. The assisted control command $\mathbf{y}_t^{\text{assist}}$, based on the optimal prediction $\mathbf{y}_t$, is $\mathbf{y}_t^{\text{assist}} = \omega_c \hat{\mathbf{y}}_t + \omega_s \mathbf{y}_t$. Here, $\omega_c$ is the weight of support provided by system, and $\omega_s$ is the patient's control weight, $\omega_s = 1 - \omega_c$. Assistance decreased progressively during the calibration phase. A maximum of 30% assistance was provided. The test phase was unassisted (figure 3).

Several series of BCI experiments of the asynchronous alternative controls of the two arms are presented in the article performing 1/ alternative two-handed reaching tasks in 3D with the virtual avatar (6D control), and 2/ alternative two-handed reaching in 3D plus 1D wrist rotation in the exoskeleton or with the virtual avatar (8D control). 6D experiments were performed between March and June 2018. During these experiments, REW NPLS algorithm previously integrated to the experimental chain was employed for real time BCI control. In the current article the 6D datasets was used for offline comparison between the generic REW-NPLS and the newly presented REW-MSLM algorithm. Experimental sessions with the patient performing 8D control were carried out using the proposed REW-MSLM algorithm. The experiments were performed between September, 26th 2018 and May 20th, 2019.
<figure><img src="image_8.png" /><figcaption>**Figure 3.** Examples of 6D alternative multi-limb pursuit tasks. One session is composed of successive tasks. Each active task is composed of several trials in which the 3D cursor must reach the proposed targets. The cursor position is not reset between tasks, during task and during idle state. Assistance optionally provided to the patient during the early model calibration phase gradually decrease during the experiment. The test phase of experiment is unassisted.</figcaption></figure>

Experimental series A

<figure><img src="image_9.png" /></figure>

B

<figure><img src="image_10.png" /></figure>

C

<figure><img src="image_11.png" /></figure>

D1 ( avatar ) / D2 ( exoskeleton )

<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /></figure>

<figure><img src="image_15.png" /><figcaption>**Figure 4.** Experimental series and decoder update strategies. Experimental series (A)–(C) of 6D avatar control: independent trining (A), cross-session training (B), model fixed (C). Experimental series of 8D control of the avatar effector (D1), and the exoskeleton control (D2). Training phases/sessions are shown in blue, testing phase/sessions are depicted in green.</figcaption></figure>

Different decoder update strategies were tested
with 6D control experiments (figure 4). For paradigm
A, the decoder was calibrated from scratch at the
beginning of each session with a small training data-
set. Sessions during the paradigm A (n = 5) were
self-contained experiments. The models were inde-
pendently created (initialized to zero), trained and

tested during the same experiment. Model adapta-
tion with multiple calibration sessions was studied
with the paradigm B. Sessions in series B (n = 4)
were performed to evaluate the importance of cross-
session training. The models were initialized to zero
in the first session. Then, the models created dur-
ing the previous session were used to initialize the
model parameters of the next session. Finally, the last model created during experimental series B were used without adaptation (paradigm C). The C series of experiments (n = 5) were performed to evaluate model robustness across time and were carried out from 9 days to 28 days after model calibration. Series of experiments A–C are closed-loop experiments using virtual avatar as an effector.

Experimental sessions of the 8D control were performed with virtual avatar or exoskeleton (figure 4). A REW-MSLM decoder was recursively trained during 6 closed-loop experiments distributed over two months and was not reupdated since then. The performance of the models were evaluated during 37 avatar experiments distributed over 5–203 days after the last model recalibration session and 15 exoskeleton experiments distributed over 0–167 days after the last model recalibration session. Five exoskeleton experiments conducted between the 62^nd^ and 63^rd^ days were excluded due to patient health issues unrelated to the study. The experimental data bases of 8D control are refered as data sets D1 and D2, and correspond to the virtual avatar and the exosceletone control, respectively.

In all the experimental series, each session was composed by successive series of tasks decided by an experimenter (figure 3). Each task corresponded to one of the available states, the IS or one of the ASs. The AS tasks were translation of the left (AS<sub>LH</sub>) and right (AS<sub>RH</sub>) hand in the 3D space in 6D control experiments. It was completed by the 1D angular rotation of the left (AS<sub>LW</sub>) and right (AS<sub>RW</sub>) wrist for 8D control experiments.

During the IS, no target was presented to the patient and the patient had to remain in a non-active state until the next instruction. Each AS task was a series point-to-point pursuit trials (figure 3). During several successive trials the patient attempted to reach the presented target locations. The target locations were set sequentially with a hand translation or a wrist rotating to a specific angle value. During a session, the hand position was not reset by the system switching between the different tasks and trials. An illustration, a session with the three states IS, AS<sub>LH</sub> and AS<sub>RH</sub> is shown in the figure 3.

The targets and the tasks were decided by the experimenter and was displayed in a random order to avoid habituation/anticipation. The training/test duration and the partition of trials and tasks varied depending on data set. For the data set D1 of 8D virtual avatar control, the total training time of the models was 3 h and 37 min including 189, 194, 181, and 218 trials for the left and right hand translation, and left and right hand rotation, respectively. For the data set D2 of 8D exoskeleton control, the calibration of 3 h and 33 min was performed to train the model including 180, 184, 188, and 226 trials for the left and right hand translation, and left and right hand rotation, respectively.

A total of 22 targets symmetrically distributed in two 3D cubes (11 targets per hand) was proposed to the patient for the left and right hand translation, in both 6D and 8D control experiments. When using exoskeleton as a training platform, the target location for each trial was given with an LED <0.01 m in diameter. For the rotation task, the LEDs were blinking in a rotating pattern to cue the direction of the wrist rotation during experiments with the exoskeleton. When training with the avatar, cubic for the left hand or spherical for the right hand targets of size of 0.1 m were displayed to indicate the destinations of each translation trial. A gauge indicated the current and target wrist positions during experiments with the avatar (figure 2).

The experimenter asked the patient to perform the mental tasks using upper limb motor imagery/ attempted movements. The patient was not instructed to perform specific motor imagery, but he was allowed to define his strategy for each movement. For a given limb, he was limited to motor imagery combining arms, wrists and fingers movements. During idle periods, the patient was asked not performing any specific tasks but take a break in the experiment. The patient was also instructed to avoid real movements using virtual avatar when residual movements were possible for a given MI. He was asked to use the same control strategy for both upper limbs. The patient was urged to maintain it constant for each task through all the experiments.

### 2.5.1. Neuronal feature extraction

During the experimental sessions, at each time step t, 1 s-long epochs of neural signals for all the electrodes $X_t \in \mathbb{R}^{586 \times 64}$ were generated with a 100 ms sliding step [17]. ECoG epochs were mapped to the temporal frequency space using a complex continuous wavelet transform (CCWT) (Morlet) with a frequency range from 10 to 150 Hz (10 Hz step) for all the electrodes. CCWT is a feature extraction strategy that was widely used in the field of BCIs. Its efficiency has previously been demonstrated [14, 15, 17, 18, 20]. The absolute value of CCWT was decimated along the temporal modality to obtain a 10-point description of a 1 s time epoch for each frequency band and for each channel, resulting in the temporal-frequency-spatial neural feature tensor $X_t \in \mathbb{R}^{10 \times 15 \times 64}$.

### 2.5.2. Output features

REW-MSLM is a supervised learning algorithm. Movement (output) features were extracted for model training. At the time t the optimal continuous movement $y_t$, and the discrete state labels $z_t \in [1; K] \subset \mathbb{N}^*$, with K the number of states, were evaluated. The output feature vector was $y_t = ((y_t^{\text{Ltr}})^T, (y_t^{\text{Rtr}})^T)^T$, $y_t \in \mathbb{R}^6$, for alternative two-handed 3D reaching tasks, and $y_t = ((y_t^{\text{Ltr}})^T, y_t^{\text{Lr}}, (y_t^{\text{Rtr}})^T, y_t^{\text{Rr}})^T$, $y_t \in \mathbb{R}^8$, if 1D wrists rotation is additionally considered. Here
<figure><img src="image_16.png" /><figcaption>**Figure 5.** Feature extraction for supervised training. Neural and movement features were extracted during the closed-loop experiments for the adaptive supervised decoder training. The movement features are the discrete state labels $z_t$, and the optimal predicted direction $\mathbf{y}_t$ defined as the 3D cartesian vector between the current position and the target position for the 3D hand translation, and as the 1D angular vector between the current angle and the target angle for 1D wrist rotation. To extract neural features at time $t$, the last second of the neural activity is mapped to the spatial frequency space using a CCWT to create tensor-shaped neural features $\mathbf{X}_t$. Neural and movement features are stored in a buffer until the next update (every 15 s) of the REW-MSLM decoder.</figcaption></figure>

$\mathbf{y}_t^{\text{Ltr}} \in \mathbb{R}^3$ and $\mathbf{y}_t^{\text{Rtr}} \in \mathbb{R}^3$ are left and right hand translation components of $\mathbf{y}_t$. They are defined as the 3D Cartesian vector between the current hand position at the time moment $t$ and the target position. The left $\mathbf{y}_t^{\text{Ltr}} \in \mathbb{R}^3$ and the right $\mathbf{y}_t^{\text{Rtr}} \in \mathbb{R}^3$ wrist rotation components of output feature vector are defined as a 1D angle between the current angle position and the target angle position (figure 5) [40]. The discrete state $z_t$ labels are determined by the task instruction. $K = 3$ in the 6D control experiments (idle state, left hand translation, and right hand translation states) and $K = 5$ (idle state, left hand translation, right hand translation, left wrist and right wrist rotation states) in the 8D experiments. Output features were extracted during experiments at 10 Hz (figure 5).

The decoder prediction ($\hat{\mathbf{y}}_t \in \mathbb{R}^6$ for 6D experiments, and $\hat{\mathbf{y}}_t \in \mathbb{R}^8$ for 8D experiments) is the cartesian position increments for 3D hands translations and angular increments for 1D wrists rotations.

### 2.5.3. REW-MSLM parameters and structure

REW-MSLM states were associated to particular tasks. In the present study, an ME structure with three states: idle (IS), left (AS<sub>LH</sub>) and right (AS<sub>RH</sub>) hand translation states, was considered in the offline comparison study to decode asynchronous alternative two-hand 3D reaching tasks. An ME structure with 5 states: idle (IS), left (AS<sub>LH</sub>) and right (AS<sub>RH</sub>) hand translation, left (AS<sub>LW</sub>) and right (AS<sub>RW</sub>) wrists rotation states was used during the online closed-loop experimental sessions using the REW-MSLM

algorithm to control the exoskeleton or the virtual avatar.

## 2.6. REW-MSLM decoder integration

The application loop for the online decoding and the adaptation loop for the REW-MSLM decoder update were split and implemented in two independent processes/threads while communicating through shared memory. The application loop received the data from the WIMAGINE implants and decoded the neural signals in order to control the effector. The input and output features were stacked in buffers to be sent to the calibration loop in order to perform the incremental batch update of the gate and expert models of REW-MSLM decoder. The data was not preselected, all the data recorded was used in the online calibration/test of the model. As series of experiment reported in the article were performed once the patient had chosen a strategy, a forgetting factor was set to keep all the data (no forgetting). In this study, neural signals were decoded at 10 Hz (a prediction was generated every 0.1 s) while the model was updated at approximately a 0.07 Hz update rate (every 15 s). Each incremental update was based on $\Delta L = 150$ samples.

The decoder prediction, the Cartesian position increments for 3D hand translation and the angular increments for 1D wrist rotation, was sent to the exoskeleton/virtual avatar after post-processing. In both cases the decoder output was saturated to a maximum of $0.1 \text{ m.s}^{-1}$ for hand translation and $1 \text{ rad.s}^{-1}$ for wrist rotation, and post-processed by the exoskeleton
control system using inverse kinematics to transform the cartesian predictions into joint movements.

Every analysis and online experiment, including training and decoding, were performed with Matlab2017b using an Intel Xeon E5-2620v3 computer with 64 GB RAM.

## 2.7. REW-MSLM performance evaluation

Experimental sessions of alternative 3D two-hand reaching tasks (6D) of the virtual avatar effector (Series A–C) were used for the offline/pseudo-online models comparison. The datasets were recorded during online closed-loop experiments using REW-NPLS decoder previously integrated to the BCI system. They were re-computed with the different algorithms in a pseudo online manner using sample-wise indicators for the performance evaluation. Pseudo-online simulation was conducted using the same training data, the same parameters (buffer size, batch training etc) and the same model application procedure as the one used for online real-time experiments to reproduce the online experiment conditions. Pseudo-online comparison is not fully generalizable for the online case due to the lack of patient's feedback. Nevertheless, it allows characterizing to some extent the algorithms before an integration into the clinical BCI decoding platform. Both discrete and continuous decoders of REW-MSLM were evaluated offline using the series of experiments A–C.

### 2.7.1. Comparison study

#### 2.7.1.1. State discrimination

We first highlighted the discrete multi-state decoding performance evaluating the accuracy of switching between all AS and, especially, the robustness of IS support. The REW-MSLM discrete decoding was compared to the REW-NPLS algorithm thresholded in post-processing (referred to as REW-NPLST) to label the continuous decoding results as discrete IS and AS states. Such a comparison underlined the benefits of computing an additional discrete decoder to inhibit the experts continuous outputs. The threshold in REW-NPLST was applied to the norm of groups of DoFs to transform the continuous output into a discrete one. The DoFs were regrouped to correspond to the states decoded by the state decoder in the REW-MSLM. The threshold value was chosen to maximize the REW-NPLST classification accuracy. Next, the REW-MSLM was compared to its own variant without HMM (referred to as REW-MSLM) to determine the benefits of dynamic HMM gating.

#### 2.7.1.2. Continuous decoding

The REW-MSLM algorithm benefits from the ME structure which splits the neural space into state-related subsets associated to independent expert decoders. The training data are divided into subsets associated with particular experts, allowing

independent expert learning. However, continuous decoder-experts are trained on a smaller specific subset of the training dataset. This may affect regression performance. The expert-specific subset training strategy was evaluated by comparing the continuous decoding performance of piece-wise linear REW-MSLM to the state-of-the-art adaptive linear regression which was trained on the entire dataset. REW-MSLM experts trained on specific subsets of the training dataset were compared to the REW-NPLS model trained on the entire dataset. Performances were compared to the REW-NPLS algorithm because it is a state-of-the-art online adaptive tensor-input tensor-output algorithm which has been previously used for closed-loop ECoG-based BCI [17, 19]. A single multilinear decoder was identified by REW-NPLS. An ME structure with three states: idle (IS), left (ASLH) and right (ASRH) hand translation states, was considered using REW-MSLM.

Finally, the REW-MSLM algorithm was integrated into the BCI platform to carry out 8D control of the virtual avatar and the exoceletone (Series D1, D2). Real-time closed-loop experiments with both effectors were carried out to evaluate the long term REW-MSLM decoder performance.

All the performance evaluations studies were performed using out of update and unassisted test experiments/periods.

### 2.7.2. Performance criteria

#### 2.7.2.1. Discrete decoding performance

Discrete decoding performance was evaluated using the accuracy (acc) and the F-score (fsc) for the multi-class case (IS versus ASLH versus ASRH versus ASLW versus ASRW if presented) for all series of experiments. Additionally, the two binary cases: IS versus ASLH and ASRH combined (named AS) and the classification between active states (ASLH versus ASRH) were considered for offline Series A–C studies. Accuracy is generally used in BCI for binary and multi-state classification [15, 16, 30, 69, 70] and is useful for performance comparison due to its ease of computation and interpretation. Nevertheless, as accuracy presents weaknesses in the case of highly imbalanced datasets, F-score is also computed. The indicators are computed from the confusion matrix, which summarizes the number of correctly classified samples from one state (true positives, tp), incorrectly labeled samples in one state (false negatives, fn), correctly classified samples not belonging to the state (true negatives, tn) and incorrectly labeled samples not belonging to the state (false positives, fp):

$$ \text{Accuracy} = \frac{1}{K} \sum_{k=1}^{K} \frac{\text{tp}_k + \text{tn}_k}{\text{tp}_k + \text{tn}_k + \text{fp}_k + \text{fn}_k}, $$

$$ \text{Fscore} = \frac{1}{K} \sum_{k=1}^{K} \frac{(\beta^2 + 1) \text{ Precision}_k \text{ Recall}_k}{\beta^2 \text{ Precision}_k + \text{Recall}_k}, $$
$$Precision_k = \frac{tp_k}{tp_k + fp_k}, \quad Recall_k = \frac{tp_k}{tp_k + fn_k}.$$

The weighting coefficient $\beta$ was set to one, the true positives $tp_k$ are considered for samples labeled as belonging to state $k$, and the true negatives $fn_k$ include those from all the other states (one versus all analysis).

Accuracy and F-score indicators are sample-based performance estimators and do not reflect the dynamic behavior of misclassified samples. Therefore, supplementary indicators were introduced. Consecutive misclassified samples were counted to evaluate the error block rate (ErrB<sub>rate</sub>) and the error block durations (ErrB<sub>time</sub>). A sequence of misclassified samples is referred as an error block. Error block durations present mean duration of error blocks in seconds. The error block rate represents the occurrence of blocks of wrong detections per minute. Additionally, the latency (lat) between the task instruction initiated by the experimenter and initiation of the movement by the patient was computed to evaluate the response time variation introduced by the HMM. The computed latency includes the patient's reaction time and the decoder latencies.

### 2.7.2.2. Continuous decoding performance

The sample-based cosine similarity (CosSim) was applied to compare the predictions of several algorithms for hand translation. At each time moment CosSim<sub>t</sub> compares the predicted $\hat{y}_t$ and the optimal prediction $y_t$. The optimal prediction is a 3D Cartesian vector between the current position and the target position for 3D translation. CosSim<sub>t</sub> is defined as the normalized scalar product of $y_t$ and $\hat{y}_t$.

$$CosSim_t = \frac{y_t \cdot \hat{y}_t}{\|y_t\| \|\hat{y}_t\|}.$$

CosSim is evaluated for a specific task as the expectation of CosSim<sub>t</sub> through related trials. CosSim varies from -1 to 1 evaluating the algorithm's global static prediction performance. CosSim was evaluated for the left hand (CosSim<sub>LH</sub>) and the right hand (CosSim<sub>RH</sub>) translation tasks. Wrist rotation performance was evaluated using Pearson correlation coefficient $R_{LW}$ and $R_{RW}$ computed for the left and the right wrist rotation tasks.

### 2.7.2.3. Online closed-loop performance evaluation

In addition to indicators previously described, the closed-loop experiment performance was evaluated using the success rate (SR) [1, 3, 19] and the R-ratio [19]. The SR is defined as the percentage of targets hit, while the R-ratio is defined as the ratio between the distance realized by the effector to reach a target and the distance from the initial position of the effector to the target location. Computed post-experiment, targets are considered to have a diameter of 10 cm with the exoskeleton as with the avatar. R-ratio [19] is also named as the distance ratio [12] and is equivalent to

the inverse of the individual path efficiency [2, 3] of each task. The SR and R-ratio performance indicators are defined in the same way for the evaluation of wrist rotation performance.

Finally, we evaluated the evolution of the performance indicators across experiments. The linear fit with a 95% confidence interval was computed for each indicator to test the zero slope hypothesis and evaluate the performance stability across time. Supplementary videos (SV1, SV2 and SV3 available online at stacks.iop.org/JNE/19/026021/mmedia) present examples of sessions 36, 106, 167 days after the last model calibration using the exoskeleton.

### 2.7.3. Chance level study

To control for potential experimental biases, the chance level of the SV performance indicators was computed and the quality of the neuronal signal recorded during the experimental sessions was evaluated.

Discrete states are not uniformly distributed, with a higher prior probability for idle and hand movements than wrist rotations (for exoskeleton-based experiments: idle, left and right hand, left and right wrist states represented 26%, 36%, 27%, 6%, 5% of the discrete state distribution, respectively). For the SR and R-ratio, $n = 100$ random hit experiments were repeated. Random movement reaching tasks were performed with the same target locations as those used during the exoskeleton-based experiments. A 3D randomly moving cursor must reach a randomly selected target within a fixed duration (defined as the 99% of the cumulative distribution of the experimental time used in the exoskeleton-based experiments). At each time step, the cursor moved in a 3D random direction with a speed fixed to the maximal speed of the exoskeleton. A target was considered a hit when the distance between the cursor and the target was less than 5 cm. These random sessions resulted in an averaged SR of 7.1 ± 5.5% (R-ratio:24 ± 14) for the left hand translation, 9.5 ± 6.6% (R-ratio: 33 ± 19) for the right hand translation, 40 ± 7.1% (R-ratio: 15 ± 4.6) for the left hand rotation and 33 ± 4.9% (R-ratio: 12 ± 2.7) for the right hand rotation tasks.

### 2.7.4. Neuronal signal recording quality evaluation

The ECoG recorded at rest prior to each experiment was analyzed to assess the signal quality over the sessions performed with an avatar or in the exoskeleton. Because of recording issues, the rest sessions recorded on day = 168 and day = 167 after the last model calibration were removed for the virtual and exoskeleton sessions respectively. A 90 s time window (from + 20 s to + 110 s post-recording onset) was used to calculate the power spectral density on the demeaned 64 electrodes using a 4th order Butterworth, IIR filter. Band-power values (dB) were computed for the whole frequency range used in the study (10–150 Hz) and for
the two frequency bands which are generally used in ECoG-driven BCI studies: 20–40 Hz and 60–110 Hz [11, 16, 71]. For each frequency band, the bandpower values were fitted to a linear regression to estimate the corresponding slope and its error-estimate with a 95% confidence interval.

### 2.7.5. Model convergence evaluation

The convergence of the models created during the online closed-loop asynchronous alternative 8D experiments using the avatar and the exoskeleton were studied. The Frobenius distance was evaluated between consecutive update of the models during the calibration period for the expert models dedicated to the 3D left and right hand translation decoding and 1D left and right wrist rotation decoding. The Frobenius distance is the generalization of the Euc- lidian distance applied to tensors.

### 2.7.6. Activation spectrograms

Motor activation spectrograms are computed from time frequency feature maps in the [-4, 4] s time interval related to the each motor state onset (t = 0) by averaging trough corresponding trials. The aver- aged (along time direction) idle state time frequency feature map is used as a baseline. The first and the last 2 s intervals of idle state epochs, potentially related to state modification, are excluded from the base line evaluation. The intervals with false activation, or erroneous state activation were removed from the analysis. The activation spectrograms of each motor state onset are computed on the sessions performed between 26 September 2018 and 20 May 2019 with the exoskeleton.

## 3. Results

### 3.1. Pseudo-online REW-MSLM 6D decoding performance evaluation

#### 3.1.1. Discrete decoding performance

The REW-MSLM demonstrated strong discrim- inative abilities (figure 6(a)) between all states (acc = 93 ± 1.8%, fsc = 86 ± 3%), between IS and AS (acc = 91 ± 3%, fsc = 84 ± 5%) and between AS<sub>LH</sub> and AS<sub>RH</sub> (acc = 99 ± 0.8%, fsc = 99 ± 0.8%) regardless of the experimental paradigm. The same performance indicators lead to acc = 87 ± 2%, fsc = 76 ± 3% between all states, acc = 86 ± 2%, fsc = 75 ± 3% between IS and AS and acc = 93 ± 0.3%, fsc = 93 ± 0.2% between AS<sub>LH</sub> and AS<sub>RH</sub> for REW-SLM algorithm whereas REW-NPLS performs acc = 62 ± 2%, fsc = 36 ± 5% between all states, acc = 70 ± 7%, fsc = 49 ± 0.6% between IS and AS and acc = 59 ± 8%, fsc = 57 ± 9% between AS<sub>LH</sub> and AS<sub>RH</sub>. The REW-MSLM strongly discrim- inated each state with a particularly robust distinction between the left and right hand. Significant improve- ments compared to REW-NPLS<sub>T</sub> and REW-SLM were evident in the majority of the performance

indicators (figure 6(a)). No significant differences between the performance in the experimental ses- sions B and C were found (p > 0.1), indicating model stability in session C, even though the model was not recalibrated in these experiments.

The latency of the switching state averaged over the three experimental paradigms (A, B and C) was higher for the REW-MSLM than for the REW-SLM: lat = 2.05 ± 0.059 s versus lat = 1.46 ± 0.31 (figure 6(b)). Similarly, the error block duration increased with the REW-MSLM decoders. The HMM state decoder error lasted ErrB<sub>time</sub> = 4.31 ± 0.88 s, whereas the discrete static decoder error duration of the REW-SLM was ErrB<sub>time</sub> = 0.49 ± 0.024 s. How- ever, the error block rate decreased considerably with the REW-MSLM decoders: the error block rate for the REW-SLM was high (ErrB<sub>rate</sub> = 20.7 ± 1.95 error blocks per minute), whereas that of the REW-MSLM was reduced to ErrB<sub>rate</sub> = 1.6 ± 0.26 blocks per minute.

#### 3.1.2. Continuous decoding performance

To evaluate expert-specific subset training strategy, piece-wise linear continuous REW-MSLM predic- tions were compared to those of the REW-NPLS decoder trained on the entire dataset. No statist- ical differences in the decoding performance were highlighted between REW-MSLM and REW-NPLS. For the paradigm A, CosSim<sub>LH</sub> = 0.095 ± 0.05, and CosSim<sub>RH</sub> = -0.03 ± 0.16 in aver- age for the REW-MSLM decoder com- pared to CosSim<sub>LH</sub> = -0.03 ± 0.14, and CosSim<sub>RH</sub> = -0.04 ± 0.1 for the REW-NPLS model (figure 7). Left hand decoding of experimental sessions B and C demonstrated equivalent average decoding performance: CosSim<sub>LH</sub> = 0.21 ± 0.06 and CosSim<sub>LH</sub> = 0.23 ± 0.13 for exper- imental sessions B and C for REW-MSLM decoder and CosSim<sub>LH</sub> = 0.18 ± 0.05 and CosSim<sub>LH</sub> = 0.18 ± 0.11 for experimental sessions B and C for the REW-NPLS model. Right hand decoding average performance of REW-MSLM (B: CosSim<sub>RH</sub> = 0.15 ± 0.07 and C: CosSim<sub>RH</sub> = 0.2 ± 0.03) is sim- ilar to the decoding performance of REW- NPLS (B: CosSim<sub>RH</sub> = 0.14 ± 0.09 and C: CosSim<sub>RH</sub> = 0.19 ± 0.03) (figure 7). Signific- ant improvements in performance between dataset A and datasets B and C highlighted the benefits of cross-session training for increasing both the train- ing data length and robustness to signal variability. No statistically significant performance differences were observed between datasets B and C, stressing the model robustness.

### 3.2. Online closed-loop REW-MSLM 8D control performance evaluation

When considering the whole frequency range used in this study (10–150 Hz), the ECoG analysis performed
<figure><img src="image_17.png" /><figcaption>**Figure 6.** State decoding results obtained during pseudo online experiments. (a) Average accuracy and F-score over datasets A, B and C for 3 different analyses: all states (idle state IS, left hand translation active state AS<sub>LH</sub> and right hand translation active state AS<sub>RH</sub>) considered independently, IS versus AS (both hand translation states merged) performance and AS<sub>LH</sub> versus AS<sub>RH</sub>. (b) Time dynamic performance indicators. Latency duration is evaluated as the time required to reach the desired state. Error block durations show the average time of the consecutive misclassified samples. The error block rate represents the occurrence of blocks of wrong detections per minute. Standard deviation is represented for each algorithm and each dataset using a vertical bar. Significance of the differences between the three decoders are computed for datasets A and C (B is excluded because of the sample size) using the Mann-Whitney U test with Bonferroni corrections ($\alpha_{multi-class} = 0.0167$) in the multi-class comparisons. Otherwise, $\alpha = 0.05$. Significant values are indicated by an asterisk.</figcaption></figure>

at rest showed a stable bandpower with a slope of -0.84% (CI = ±0.61%) and -0.99% (CI = ±0.84%) for the avatar and the exoskeleton experiments respectively. A similar trend was observed for the two frequency ranges which were mostly used by the decoder: the 20-40 Hz band with the respective slopes of -0.97% (CI = ±0.58%) and -0.75% (CI = ±0.59%) for the avatar and the exoskeleton experiments, and the 60-110 Hz band with the respective slopes of -0.23% (CI = ±0.26%) and -0.13% (CI = ±0.65%) for the avatar and the exoskeleton (figure 8).

### 3.2.1. Virtual avatar control performance

High classification decoding performance discriminating five states (idle, left and right hands translation

and left and right wrist rotation) was demonstrated with the REW-MSLM algorithm across all the experiments (figure 9(a)) with an average (across states and experiments) F-score of fsc = 76 ± 9% and accuracy of acc = 93 ± 3% (figure 10).

Mean continuous performances are given in table 1. It can be observed that the results are very close for both hands in the case of hand translation. In the case of wrist rotation SR is very close to 100% for both hands, R-ratio was higher in the case of left wrist (4.5 ± 2.2 vs 2.8 ± 1.4) and, coherently, correlation was lower for left wrist (0.391 ± 0.135 vs 0.450 ± 0.0966).

The zero slope hypothesis was not rejected for 20 of the 22 indicators (table TS1 in supplementary materials). It was rejected for the left wrist rotation
<figure><img src="image_18.png" /><figcaption>**Figure 7.** Continuous decoding performance for each hand for datasets (A)–(C). Statistics of the scalar product between the predicted hand directions and the optimal prediction (defined as the target-cursor oriented distance) averaged over time and the experiments of each dataset. The scalar products are represented for the left and right hand performance. The performance indicators are shown in blue for the state-of-the-art REW-NPLS model and in yellow for the new REW-MSLM. On each box, the central mark indicates the median, and the bottom and top edges of the box indicate the 25th and 75th percentiles, respectively. The whiskers extend to the extreme data.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>**Figure 8.** Evolution of the mean bandpower values (dB ± sd) over days for the avatar (a) and the exoskeleton (b) experiments, with day = 0 being the last day the model calibration was updated. The bandpower values were computed for three frequency ranges of interest: 10–150 Hz, 20–40 Hz and 60–110 Hz.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>**Figure 9.** Confusion matrices corresponding to the evaluation of the avatar (a) and exoskeleton (b) gate models on their respective test set.</figcaption></figure>
<figure><img src="image_21.png" /><figcaption>**Figure 10.** Online 8 DoF experiment performance (for each state: idle, left and right hand translation and rotation) using virtual avatar effector across 203 days after last model calibration (a) or using exoskeleton effector across 167 days after last model calibration (b). Discrete performances are represented for each state with F-score and accuracy. Continuous performances are given by the success rate (SR) (percentage of targets hit), the R-ratio (ratio between the distance travelled by the effector to reach a target and the distance from the initial position of the effector to target location), Pearson correlation and cosine similarity. Standard deviation is represented for each algorithm and each dataset using a vertical bar. In the case of cosine similarity, for clarity purposes, the vertical bars correspond to standard deviation divided by 2.</figcaption></figure>
<table><caption>**Table 1.** Mean continuous performance indicator in the case of avatar.</caption><thead><tr><th></th><th>SR (%)</th><th>R-ratio</th><th>Cosine Similarity</th><th>Pearson Correlation</th></tr></thead><tbody><tr><td>AS<sub>LH</sub></td><td>55 ± 18</td><td>5.2 ± 3.1</td><td>0.146 ± 0.108</td><td></td></tr><tr><td>AS<sub>RH</sub></td><td>53 ± 15</td><td>5.4 ± 3.5</td><td>0.137 ± 0.0599</td><td></td></tr><tr><td>AS<sub>LW</sub></td><td>95 ± 9.5</td><td>4.5 ± 2.2</td><td></td><td>0.391 ± 0.135</td></tr><tr><td>AS<sub>RW</sub></td><td>96 ± 6.8</td><td>2,8 ± 1,4</td><td></td><td>0.450 ± 0.0966</td></tr></tbody></table>

<table><caption>**Table 2.** Mean continuous performance indicator in the case of exoskeleton.</caption><thead><tr><th></th><th>SR (%)</th><th>R-ratio</th><th>Cosine Similarity</th><th>Pearson Correlation</th></tr></thead><tbody><tr><td>AS<sub>LH</sub></td><td>69 ± 13</td><td>6.7 ± 5.4</td><td>0.129 ± 0.0865</td><td></td></tr><tr><td>AS<sub>RH</sub></td><td>65 ± 29</td><td>13 ± 4.5</td><td>0.105 ± 0.140</td><td></td></tr><tr><td>AS<sub>LW</sub></td><td>90 ± 14</td><td>3.4 ± 1.2</td><td></td><td>0.556 ± 0.110</td></tr><tr><td>AS<sub>RW</sub></td><td>97 ± 6.2</td><td>2.4 ± 0.73</td><td></td><td>0.562 ± 0.109</td></tr></tbody></table>

R-ratio, which increased by 0.014 daily, and the right hand translation SR, which reduced daily by 0.07%. These results highlight the stability of the REW-MSLM over 6.5 months using a virtual avatar effector during 8D experiments.

### 3.2.2. Exoskeleton control performance

The discrete decoding performances of 8D exoskeleton control experiments yielded relevant and stable results across the 167 days. The REW-MSLM's gating yielded an average F-score of 75 ± 12% and accuracy of 92 ± 4% with high distinctiveness between the classification of the left and right sides of the body (less than 1% of misclassified samples) and strong idle state decoding with an average of 85% accurately classified idle state samples. The confusion matrix is depicted at figure 9(b).

Mean continuous performances are given in table 2. In comparison to left hand translation, right hand translation had a similar SR but a higher R-ratio (13 ± 4.5 vs 6.7 ± 5.4) and a lower cosine similarity (0.129 ± 0.0865 vs 0.105 ± 0.140), which means control of the right arm was less stable. Performances for both wrist indicates a very high SR and fast target reaching.

It is worth to note, that for the period 0–37 days after the last decoder calibration session, the online sessions using the exoskeleton yielded a decoding accuracy of 94% averaged across the five classes. Additionally, an average SR for both hands of 83% and 97% with an average R-ratio of 6.4 ± 2.3 and 3.3 ± 1.7 for the 3D hand translation and 1D wrist rotation were reported for 8D control on the same period. This period corresponds or over-passes the time interval reported generally in ECoG-based BCI studies. Commonly, ECoG based clinical trials last from several days to 1 or two weeks of research with an implantation from 3 to 35 days [8, 10–13, 16, 21, 22, 72–74].

The decoding stability was evaluated with zero slope hypothesis, which was not rejected for 17 of the 22 indicators (table TS2 in supplementary materials). The right side of the body seemed to have a

slow performance decrease across experiments, gathering 5 of the 6 diminishing indicators. The linear fits demonstrated significant decreases in right limb performance for the discrete right wrist rotation indicators (−0.25% F-score and −0.04% accuracy per day) and for the right hand translation F-score (−0.17%) and SR (−0.48%). Significant decreases were found in the left hand SR (−0.18% per day). The left hand SR seemed to decay in the first experiments before stabilizing.

All the 22 performance indicators had higher values than those obtained by chance level studies for all the experiments. Chance level studies highlighted an averaged SR of 7.1 ± 5.5% (R-ratio: 24 ± 14) and 9.5 ± 6.6% (R-ratio: 33 ± 19) for the left and right hand translation respectively whereas left and right hand rotation tasks chance level was evaluated at, SR = 40 ± 7.1% (R-ratio: 15 ± 4.6) and SR = 33 ± 4.9% (R-ratio: 12 ± 2.7).

Figure 11 illustrates the convergence through the model update iterations during 6 calibration session of coefficients of the expert models for the left and right hand translation and left and right wrist rotation decoding.

Examples of hand trajectories performed on the session 106 days after the model calibration are presented in figures 12(a) and (b) for the left and right hand translation, respectively. Additional trajectories are proposed in the supplementary materials (figure S1) and supplementary videos (SV1, SV2 and SV3) present examples of sessions 36, 106, 167 days after the last model calibration. The entire session of the 106th days is represented in figure 12(c). This session is composed of successive tasks with a total of two right hand translation tasks and three idle, left hand translation, left and right hand rotation tasks. Each tasks is composed of several trials. Trajectories from figures 12(a) and (b) are trials from the first left hand and second right hand translation tasks.

Gating models used for exoskeleton control is represented on the spatial, frequency and temporal modality in figure 13. For comparison purposes, the avatar model can be found in the supplementary materials,
<figure><img src="image_22.png" /><figcaption>**Figure 11.** Convergence through the model update iterations of coefficients of the expert decoding models of left hand translation (a), right hand translation (b), left wrist rotation (c), right wrist rotation (d). The Frobenius distance between consecutive coefficients update is depicted by coloured lines for a set of updated models including *f* factors, *f* = 1, ..., 100. Models selected by online validation procedure are depicted by the bold purple line. Calibration sessions are separated by the black vertical dotted lines.</figcaption></figure>

in figure S2. Variability of the exoskeleton model coefficients according to the different modalities are presented in the supplementary materials in figure S3.

The presentation of decoding model, using the sum of the absolute values of coefficients according to tensor dimensions, allows characterizing at some extend the decoding model. In the same time, they are not fully representative. As a complement, the motor activation spectrograms were computed for each motor state (figure 14).

## 4. Discussion

Based on the first successful long-term (more than 36 months) chronic exploitation of bilateral epidermal ECoG recordings in a tetraplegic individual [19], we developed the REW-MSLM decoder to address the poorly explored field of asynchronous multi-limb effector control. This decoder was designed to overcome the major issues related to the translation of BCIs from the laboratory to real-life applications, such as the high-dimensional control of effectors based on chronic neural recordings, experiments closer to real life behavior, and the ability of the BCI system to act as a stand-alone device switching between IS and AS phases. ME architecture was employed to handle numerous dimensions and to decode the robust idle state. We developed an adaptive/incremental learning algorithm which allows cross-session training of the decoder with multiple recording conditions during closed loop BCI experiments taking into account neuronal feedback. The incremental learning algorithm is able to handle in real-time high-dimensional tensor data flow. Tensor-based algorithms emerged as a promising strategy for brain signal processing allowing

simultaneous treatments of high-dimensional data in the temporal, frequency and spatial domains [17, 59]. Dynamic gating procedure using a HMM was added to ME decoder to ensure the robustness of states. The proposed algorithm is fully adaptive. It learns the model in an incremental manner, including the hyperparameters.

The REW-MSLM was integrated into a custom-made BCI adaptive brain signal decoder (ABSD)[19] software platform to provide a tetraplegic patient with the control of a virtual avatar and an exoskeleton. Volitional alternating rotation and 3D translation movements of both hands could be executed. This performance was achieved using EpiCoG recordings which are less invasive than the subdural ECoG recordings reported in most of the BCI studies. High dimensional control of complex effectors was stable for several months. These results were obtained in the period from 14 to 23 months after the recording device implantation, demonstrating the longevity of the implant and the longevity of high quality neural data recording [66].

### 4.1. ME structure benefits for multi-limb effectors control

The REW-MSLM architecture fits the multi-limb paradigm. Each expert can be associated to a particular limb or action while the HMM gating model aims to establish the state selection and handle robust idle state detection for complex asynchronous state decoding [15]. Moreover, we hypothesized, based on well-established neurophysiological knowledge [18, 29], that neural data associated with each limb or specific action can be partitioned into different neural regions/patterns. Consequently, each expert was only trained on a small subset of the entire
<figure><img src="image_23.png" /></figure>

training dataset. Such training allowed individual
update of the experts, and incrementally appending
new experts to control new dimensions without full
re-training of other experts. To demonstrate the rel-
evance of the ME model structure and the importance
of dynamic vs. static gating, we compared the REW-
MSLM to the state-of-the-art adaptive algorithms.

The comparison of several algorithms is a con-
ventional tool to conclude on the improvements
obtained with the proposed algorithms. However,
comparing several online algorithms during closed-
loop experiments is a complicated task as, during

such experiments, the predicted trajectories are
related to the current decoding model and patient's
feedback. Consequently, it is not possible to com-
pare in online closed-loop experiments several
algorithms that produce different predictions and
feedbacks. Several series of online closed loop ses-
sions are particularly time consuming. In the current
study, offline comparison study was undertaken in
pseudo-online manner with 3 databases. The data-
sets were recorded using 3 different experimental
paradigms (single session decoder training, cross ses-
sion decoder training, fixed decoder) during online
<figure><img src="image_24.png" /><figcaption>**Figure 13.** Example of a gating model for exoskeleton control. Gating parameter weights (discrete decoding) of the REW-MSLM according to the spatial (A), frequency (B) or temporal (C) modalities for each state: rest state (IS), left hand 3D translation and rotation state (AS<sub>LH</sub> and AS<sub>LW</sub>) or right hand 3D translation and rotation state (AS<sub>RH</sub> and AS<sub>RW</sub>). The sensory sulcus (SS) and motor sulcus (MS) are represented in the spatial domain in yellow and red curves respectively. The spatial modality shows, as expected, strong parameter weights on the contralateral electrode array for left and right hand (translation and rotation) states; in addition, translation and rotation from the same hand seem to activate nearby but distinct electrodes. The frequency modality highlights the beta and high gamma frequency bands as relevant frequencies for state decoding, whereas the temporal modality emphasis parameters from 0.5 s to 0.1 s. The variability of the state decoding model coefficients according to the time and frequency modalities is represented in the supplementary materials (figure S3).</figcaption></figure>

closed-loop experiments using conventional REW-NPLS algorithm as the decoder. While offline studies gave an initial overview of the potential REW-MSLM decoding performance and benefits, they were not generalizable due to lack of the appropriate user feedback. Nevertheless, it allowed characterizing the studied algorithms before an integration into the clinical BCI decoding platform.

For discrete decoding, the REW-MSLM outperformed alternative approaches in discrete classification regardless of the dataset and paradigm with an averaged F-score improvement across all paradigms of 39 ± 4% and 8.3 ± 2% compared

to the REW-NPLS and the REW-SLM respectively. These results confirmed the benefits to train a specific model dedicated to state classification and the improvements related to dynamic classification. The switching state latency study related to the state transition delay between the instruction and the discrete decoding response demonstrated an average increase in duration by 0.45 s, 0.87 s and 0.38 s (over 3 datasets) between the discrete decoder with and without dynamic HMM processing. However, the REW-MSLM results show a drastic 92% decrease in the error block rate between the discrete decoder with and without dynamic HMM processing, overcoming
<figure><img src="image_25.png" /><figcaption>**Figure 14.** Activation spectrograms. Modulation of the neural signal activity between idle state and each active state for the electrodes selected with a checkerboard pattern. The sensory sulcus (SS) and motor sulcus (MS) are shown with yellow and red curves respectively. The dotted squares highlight the electrodes enlarged on the bottom right.</figcaption></figure>

the high frequency misclassified sample issue of static classifier. For physical effectors, such as an exoskeleton, which are in direct contact with the patient and has a latency of mechanical activation/deactivation

of up to a few seconds, false activation should remain exceptionally rare events.

For continuous control, REW-MSLM experts highlighted slight improvement or similar
performance compared to REW-NPLS, whereas the training datasets were different. REW-MSLM allows experts training using independent data sets. This may be highly profitable for progressive BCI decoder training increasing the tasks complexity.

In addition, the developed REW-MSLM and the state of the art REW-NPLS algorithms demonstrated similar decoding performance. However, numerous non-desired movements of the other limbs are observed using REW-NPLS. Unintended movements of a limb that the patient does not want to move impede the control of complex effectors such as an exoskeleton, and especially the asynchronous control with an idle state to decode. In contrast, REW-MSLM demonstrated similar decoding performance for the limb to be activated without unintentional movements from the other limbs thanks to accurate state classification. The suppression of the unintended movements leads to better visual feedbacks and concentration of the patient which may induce better model calibration. This is illustrated in the supplementary video SV4.

## 4.2. EpiCoG based neural decoder for complex tasks completion

The control tasks, proposed to the patient during the experiments are more challenging than the usual state-of-the-art control tasks. Center-out tasks require to go from the center of a workspace to one of the targets localized at equal distances. Moreover, after each trial (succeeded our failed), the position is reset to the initial position after few seconds of rest. In the point-to-point pursuit task experimental paradigm reported in this article, the patient controls the effector all along the session and without resetting the hand position. This control task is more complex because the initial position of the hand changes constantly, and decoding mistake/drifting of the hand from one trial affect the following trial. A point-to-point pursuit task is more complex compared to conventional center-out tasks in terms of explored space due to multiple (arbitrary under the constraints of control region) possible starting points and numerous targets. In the current study, 22 target positions, 11 for each hand, are proposed to the user, while a majority of center-out experiments consider 8 equally distanced targets [11, 12, 16, 75]. Point-to-point pursuit tasks are more representative of daily life applications, and have less restricted experimental conditions. In addition, asynchronous and alternative bimanual point-to-point pursuit experiments support rest period as well as asynchronous switch between active control tasks without external intervention. All dimensions of control (eight in general) are available to the user at any moment. While not all DoF may be active simultaneously, any point in the control region (8D) may be achieved by user at his own intention.

In the beginning of the experiments (out of experimental sessions presented in the manuscript), the patient optimized the motor imagery strategy to allow controlling the effectors. He reported that he was able after several months of training to control the effectors unconsciously, without focusing on motor imagery.

## 4.3. Closed-loop decoder stability using EpiCoG recordings

Generally ECoG-based BCI studies are performed using temporary ECoG subdural grid implantation from 3 to 35 days post-surgery [8, 10–13, 16, 21, 22, 72–74]. In our experiments, the online closed-loop SR for both effectors realized from 0 to 37 days after the last model calibration (figure 10) are similar or above current ECoG-based state of the art performance for 3D decoding. Importantly, compared to these subdural ECoG studies we did not perform any model recalibration during this time period even though we used a system which is less invasive [12].

The online closed-loop results presented a high stability level and were far above the realized chance level study across all experiments for both effectors. For the exoskeleton experiments, the left hand translation SR seemed to decay between the 37th and the 104th day and stabilize until the end, whereas the right hand translation SR showed higher variability in the performance (between 17% and 100%). For discrete decoding, switching from left arm control to right arm control (and vice versa) produced less than 1% of the errors. Most of the decoding misclassifications were related to two issues. First, the majority of the mistakes were related to false positive idle state activation. Second, the decoders struggled to differentiate between rotation and translation from the same limb. These difficulties may be related to the similarity of both tasks and may consequently lead to brain neural signal pattern activations within close areas.

Our results seem to demonstrate higher average performance to control the exoskeleton than the virtual avatar. This could be explained by the fact that the exoskeleton provides a more realistic feedback to the patient than the virtual avatar. However, it is difficult to make any conclusion due to the small number of experiments considered in this study.

The online control of both effectors was maintained, without recalibration, over six months of clinical experiments (for 167 days and during 203 days for the exoskeleton and virtual avatar effectors, respectively), indicating the stability of both the REW-MSLM decoder and the neural activity recording method with the two WIMAGINE EpiCoG recording implants [47]. These results show that this system outperforms the state-of-the-art ECoG-based BCIs, and outperforms both the state-of-the-art ECoG and MEAs-based BCIs in terms of decoder stability.

The pseudo online study induces the benefits of cross-session training for obtaining a better decoder,
more robust to brain and experimental condition variability. Indeed, continuous performance was low for dataset A (single session decoder training). Results from dataset C (fixed decoder) showed stable performance whereas the model was trained on the basis of cross-session calibration procedure from dataset B recorded 9–28 days before. In the online study, the REW-MSLM was trained for each effector based on cross-session calibration procedure for 6 experiments over six days, distributed over two months, for approximately 3.5 h (with in averaged 195 trials). The duration of the model training periods seems moderate, considering the high number of dimensions to control and performance obtained compared to those in similar studies [11, 12]. More training data may lead to a more generalized model and thus, better results.

Figure 13 illustrates gating model weights in the frequency, temporal and spatial modalities for exoskeleton control. Spatial modality presents heavy parameter weights on the contralateral electrode array for left and right hand (translation and rotation) states. Both translation tasks present similar model with dominant frequency band between 20–30 Hz (β-band) and 80 Hz–120 Hz (γ-band). The same frequency band are relevant for rotation and idle state model, nevertheless, lower frequency band (<20 Hz) significantly contribute to the decoding, especially for idle state decoding. Parameter weights in the temporal modalities are similar for all states, emphasizing parameters between 0.5 s and 0.1 s. In the frequency modality the model coefficients are consistent with the previous studies which highlighted the significance of β and high γ-band to decode movements from direct neural signals [16, 71, 73]. As expected, spatial weights were higher in the contralateral electrodes of the realized movement for both left and right hand translation and rotation which is corroborated by previous studies [71, 76, 77]. Compared to the exoskeleton gating model, the avatar gating model focuses on the β-band and presents a less clear contralateral importance of the electrodes.

## 4.4. Limitations and perspectives

The current paper reports the long-term stability of high dimensional (8D) control of bimanual exoskeleton and its avatar. While the study demonstrates promising results, they were demonstrated for a single patient. The implantation of 3 more patients is planned in the ‘BCI and tetraplegia’ clinical trial protocol and would provide more data to support the conclusions of this article.

In the paper, the alternative bimanual point-to-point pursuit experiments supporting asynchronous switch between control tasks without external intervention is demonstrated. The study reports an experimental paradigm less restrictive in term of experimental conditions, compared to traditional

center-out tasks, with a wider exploration of the control space. In the same time, experiments closer to domestic, urban, and professional environments are still needed to move the technology from clinical trials to daily life applications.

During experiments, the patient was able to switch reliably between idle/rest state and multiple active states. During idle periods, the patient was not asked to perform specific tasks but rather take a break in the experiment, interacting with people around him for instance, hereby approaching to real-life situation. In the same time, due to restriction in experimental session duration, the rest periods were relatively short. Further studies with longer idle state periods are still needed to demonstrate the performance required for the real-life asynchronous BCI applications.

Only alternative bimanual control was performed due to experimental paradigm. However, simultaneous bimanual control is theoretically possible thanks to REW-MSLM soft gating strategy: the gating is not a selection of one limb among the others but the mixing of all of them depending on the probability of limb activation computed by the HMM gating. Simultaneous bimanual effector control is a nearest perspective of the study.

The REW-MSLM is fully adaptive/incremental algorithm including online optimization of its hyper-parameters such as latent space dimension. In the same time, feature extraction step is fixed from previous research and may be further optimized. In particular, uniform frequency resolution used in the paper may be suboptimal. An evaluation of the potential benefits of using a higher frequency resolution is proposed in the supplementary materials (figure S4). We investigated offline the possible benefits of using a doubled resolution in frequencies <50 Hz where low resolution may mostly affect the results. Multilinear models were optimized with features computed with wavelets whose central frequency varied of 5 Hz from 5 to 50 Hz and of 10 Hz from 50 to 150 Hz, and compared models computed using features extracted in frequency bands from 10 to 150 Hz with a step of 10 Hz as it was used in real time online experiments. Results showed that doubled resolution in frequency band <50 Hz might slightly improve the cosine similarity for hand translations. However, no improvement in correlation was observed for wrists rotation tasks. Then, a higher frequency resolution might be beneficial for complex tasks such as 3D hand translations. For more simple actions models could struggle to find relevant features using lower frequency resolution. Optimization of feature extraction for each task will be further explored.

Faster and straighter reaching trajectories are likely to be particularly profitable for patients. Various post-processing strategies will be investigated in future studies to provide better control and feedback to the patient. As the drop of decoding performance
in the target neighborhood is regularly observed in BCI studies, alternative ME architecture with states associated to movement phases will be explored [61].

A restricted dataset was used for decoders training. Decoding models were fixed without determining an optimal training time. More training data may lead to a more generalized model and better results. The optimization of training paradigm, evaluation of the impact of a larger dataset on decoding performance, is one of the perspectives of the presented research. Model interpretation and convergence will be further investigated. In addition, patients' adaptations and improvements will be analyzed to evaluate the impact of experiment frequency on performance stability.

The REW-MSLM benefits from an MEs architecture, which splits the dataset to train particular experts. Continuous decoders are responsible for a single or group of dimensions. This structure allows us mixing experts from different training sets and different models or adding new dimensions without retraining all the experts. Increasing control complexity by adding dimensions sequentially is highly profitable for patient training. The mixture of experts decoder architecture favors further increase of control dimensions. Doubling the resolution of the recording system is expected in near future and may allow an increase in number of degrees of freedom, which is crucial for BCI medical use.

## Data availability statement

The data generated and/or analysed during the current study are not publicly available for legal/ethical reasons but are available from the corresponding author on reasonable request.

## Acknowledgments

CLINATEC is a Laboratory of CEA-Grenoble and has statutory links with the University Hospital of Grenoble (CHUGA) and with University Grenoble Alpes (UGA). This study was funded by CEA (recurrent funding) and the French Ministry of Health (Grant PHRC-15-15-0124), Institut Carnot, Fonds de Dotation CLINATEC. Fondation Philanthropique Edmond J Safra is a major founding institution of the CLINATEC Edmond J Safra Biomedical Research Center.

We thank Marion Mainsant, Antoine Lassauce, Gaël Reganha, Vincent Rouanne, Benoit Milville, Jean-Claude Royer, Stéphane Pezzani and the CHUGA members for the help in the experimental setup and clinical trial experiments.

## Author contributions

A M and T A developed and implemented the algorithm and software and wrote the manuscript.

T C did the experiments. A M, F M, M M, and C L analyzed the data. S K, A M and F M implemented the algorithm to the exoskeleton platform, G C supervised the design and utilization of the recording platform, A V designed the exoskeleton, S C and A L B did the surgery, A L B supervised the Brain Computer Interface CLINATEC program, all authors reviewed the manuscript.

## Declaration of interests

We declare no competing interests.

## ORCID iDs

Alexandre Moly https://orcid.org/0000-0003-1781-3537

Thomas Costecalde https://orcid.org/0000-0003-2216-4447

Félix Martel https://orcid.org/0000-0003-4877-1771

Matthieu Martin https://orcid.org/0000-0001-5954-8087

Christelle Larzabal https://orcid.org/0000-0003-4651-9265

Alexandre Verney https://orcid.org/0000-0001-5093-4918

Guillaume Charvet https://orcid.org/0000-0003-4938-9419

Stephan Chabardes https://orcid.org/0000-0002-7930-1476

Alim Louis Benabid https://orcid.org/0000-0002-4479-1807

Tetiana Aksenova https://orcid.org/0000-0003-4007-2343

## References

[1] Hochberg L R et al 2012 Reach and grasp by people with tetraplegia using a neurally controlled robotic arm *Nature* **485** 372–5

[2] Collinger J L et al 2013 High-performance neuroprosthetic control by an individual with tetraplegia *Lancet* **381** 557–64

[3] Wodlinger B et al 2015 Ten-dimensional anthropomorphic arm control in a human brain–machine interface: difficulties, solutions, and limitations *J. Neural Eng.* **12** 016011

[4] Murphy M D, Guggenmos D J, Bundy D T and Nudo R J 2016 Current challenges facing the translation of brain computer interfaces from preclinical trials to use in human patients *Front. Cell Neurosci.* **9** 497

[5] Ward M P, Rajdev P, Ellison C and Irazoqui P P 2009 Toward a comparison of microelectrodes for acute and chronic recordings *Brain Res.* **1282** 183–200

[6] Perge J A et al 2013 Intra-day signal instabilities affect decoding performance in an intracortical neural interface system *J. Neural Eng.* **10** 036004

[7] Sussillo D, Stavisky S D, Kao J C, Ryu S I and Shenoy K V 2016 Making brain–machine interfaces robust to future neural variability *Nat. Commun.* **7** 13749

[8] Leuthardt E C, Schalk G, Wolpaw J R, Ojemann J G and Moran D W 2004 A brain–computer interface using electrocorticographic signals in humans *J. Neural Eng.* **1** 63

[9] Leuthardt E C, Miller K J, Schalk G, Rao R P N and Ojemann J G 2006 Electrocorticography-based brain
computer Interface-the seattle experience *IEEE Trans. Neural Syst.* Rehabil. Eng. **14** 194–8

[10] Schalk G and Leuthardt E C 2011 Brain-computer interfaces using electrocorticographic signals *IEEE Rev. Biomed. Eng.* **4** 140–54

[11] Wang W *et al* 2013 An electrocorticographic brain interface in an individual with tetraplegia *PLoS One* **8** e55344

[12] Degenhart A D *et al* 2018 Remapping cortical modulation for electrocorticographic brain–computer interfaces: a somatotopy-based approach in individuals with upper-limb paralysis *J. Neural Eng.* **15** 026021

[13] Schalk G *et al* 2008 Two-dimensional movement control using electrocorticographic signals in humans *J. Neural Eng.* **5** 75

[14] Shimoda K, Nagasaka Y, Chao Z C and Fujii N 2012 Decoding continuous three-dimensional hand trajectories from epidural electrocorticographic signals in Japanese macaques *J. Neural Eng.* **9** 036015

[15] Schaeffer M-C and Aksenova T 2016 Switching Markov decoders for asynchronous trajectory reconstruction from ECoG signals in monkeys for BCI applications *J. Physiol.* **110** 348–60

[16] Bundy D T, Pahwa M, Szrama N and Leuthardt E C 2016 Decoding three-dimensional reaching movements using electrocorticographic signals in humans *J. Neural Eng.* **13** 026021

[17] Eliseyev A *et al* 2017 Recursive exponentially weighted N-way partial least squares regression with recursive-validation of hyper-parameters in brain-computer interface applications *Sci. Rep.* **7** 16281

[18] Choi H *et al* 2018 Improved prediction of bimanual movements by a two-staged (effector-then-trajectory) decoder with epidural ECoG in nonhuman primates *J. Neural Eng.* **15** 016011

[19] Benabid A L *et al* 2019 An exoskeleton controlled by an epidural wireless brain–machine interface in a tetraplegic patient: a proof-of-concept demonstration *Lancet Neurol.* **18** 1112

[20] Chao Z C, Nagasaka Y and Fujii N 2010 Long-term asynchronous decoding of arm motion using electrocorticographic signals in monkey *Front. Neuroeng.* **3** 3

[21] Yanagisawa T *et al* 2012 Electrocorticographic control of a prosthetic arm in paralyzed patients *Ann. Neurol.* **71** 353–61

[22] Nakanishi Y *et al* 2013 Prediction of three-dimensional arm trajectories based on ECoG signals recorded from human sensorimotor cortex *PLoS One* **8** e72085

[23] Nurse E S *et al* 2018 Consistency of long-term subdural electrocorticography in humans *IEEE Trans. Biomed. Eng.* **65** 344–52

[24] Sauter-Starace F *et al* 2019 Long-term sheep implantation of WIMAGINE, a wireless 64-channel electrocorticogram recorder *Front. Neurosci.* **13** 847

[25] Pels E G M *et al* 2019 Stability of a chronic implanted brain-computer interface in late-stage amyotrophic lateral sclerosis *Clin. Neurophysiol.* **130** 1798–803

[26] Vansteensel M J *et al* 2016 Fully implanted brain–computer interface in a locked-in patient with ALS *New Engl. J. Med.* **375** 2060–6

[27] Williams J J, Rouse A G, Thongpang S, Williams J C and Moran D W 2013 Differentiating closed-loop cortical intention from rest: building an asynchronous electrocorticographic BCI *J. Neural Eng.* **10** 046001

[28] Müller-Putz G R, Scherer R, Pfurtscheller G and Rupp R 2006 Brain-computer interfaces for control of neuroprostheses: from synchronous to asynchronous mode of operation/brain-computer interfaces zur steuerung von neuroprothesen: von der synchronen zur asynchronen Funktionsweise *Biomed. Tech.* **51** 57–63

[29] Iftt P J, Shokur S, Li Z, Lebedev M A and Nicolelis M A L 2013 A brain-machine interface enables bimanual arm movements in monkeys *Sci. Trans. Med.* **5** 210ra154

[30] Hotson G *et al* 2016 Individual finger control of a modular prosthetic limb using high-density electrocorticography in a human subject *J. Neural Eng.* **13** 026017

[31] Elgharabawy A and Wahed M A 2016 Decoding of finger movement using kinematic model classification and regression model switching 2016 8th Cairo Int. Biomedical Engineering Conf. (CIBEC) pp 84–89

[32] Flamary R and Rakotomamonjy A 2012 Decoding finger movements from ecog signals using switching linear models *Front. Neurosci.* **6** 29

[33] Lebedev M A and Nicolelis M A L 2006 Brain–machine interfaces: past, present and future *Trends Neurosci.* **29** 536–46

[34] Orsborn A L *et al* 2014 Closed-loop decoder adaptation shapes neural plasticity for skillful neuroprosthetic control *Neuron* **82** 1380–93

[35] Jarosiewicz B *et al* 2013 Advantages of closed-loop calibration in intracortical brain–computer interfaces for people with tetraplegia *J. Neural Eng.* **10** 046012

[36] Dangi S *et al* 2014 Continuous closed-loop decoder adaptation with a recursive maximum likelihood algorithm allows for rapid performance acquisition in brain-machine interfaces *Neural Comput.* **26** 1811–39

[37] Brandman D M *et al* 2018 Rapid calibration of an intracortical brain–computer interface for people with tetraplegia *J. Neural Eng.* **15** 026007

[38] Shanechi M M *et al* 2017 Rapid control and feedback rates enhance neuroprosthetic control *Nat. Commun.* **8** 13825

[39] Li Z, O'Doherty J E, Lebedev M A and Nicolelis M A L 2011 Adaptive decoding for brain-machine interfaces through bayesian parameter updates *Neural Comput.* **23** 3162–204

[40] Gilja V *et al* 2012 A high-performance neural prosthesis enabled by control algorithm design *Nat. Neurosci.* **15** 1752

[41] Vidaurre C, Kawanabe M, von Bünau P, Blankertz B and Müller K R 2011 Toward unsupervised adaptation of LDA for brain–computer interfaces *IEEE Trans. Biomed. Eng.* **58** 587–97

[42] Nicolas-Alonso L F, Corralejo R, Gomez-Pilar J, Álvarez D and Hornero R 2015 Adaptive semi-supervised classification to reduce intersession non-stationarity in multiclass motor imagery-based brain–computer interfaces *Neurocomputing* **159** 186–96

[43] Lotte F *et al* 2018 A review of classification algorithms for EEG-based brain–computer interfaces: a 10 year update *J. Neural Eng.* **15** 031005

[44] Hazrati M K and Erfanian A 2010 An online EEG-based brain–computer interface for controlling hand grasp using an adaptive probabilistic neural network *Med. Eng. Phys.* **32** 730–9

[45] Rong H, Li C, Bao R and Chen B 2018 Incremental adaptive EEG classification of motor imagery-based BCI *2018 Int. Joint Conf. on Neural Networks (IJCNN)* pp 1–7

[46] Milekovic T *et al* 2018 Stable long-term BCI-enabled communication in ALS and locked-in syndrome using LFP signals *J. Neurophysiol.* **120** 343–60

[47] Mestais C S *et al* 2015 WIMAGINE: wireless 64-Channel ECoG recording implant for long term clinical applications *IEEE Trans. Neural Syst.* Rehabil. Eng. **23** 10–21

[48] Yuksel S E, Wilson J N and Gader P D 2012 Twenty years of mixture of experts *IEEE Trans. Neural Networks Learn. Syst.* **23** 1177–93

[49] Chen C *et al* 2013 Prediction of hand trajectory from electrocorticography signals in primary motor cortex *PLoS One* **8** e83534

[50] Eliseyev A and Aksenova T 2014 Stable and artifact-resistant decoding of 3D hand trajectories from ECoG signals using the generalized additive model *J. Neural Eng.* **11** 066005

[51] Maleki M, Manshouri N and Kayikcioğlu T 2018 Fast and accurate classifier-based brain-computer interface system using single channel EEG data 2018 26th Signal Processing and Communications Applications Conf. (SIU) pp 1–4

[52] Trejo L J, Rosipal R and Matthews B 2006 Brain-computer interfaces for 1-D and 2D cursor control: designs using
volitional control of the EEG spectrum or steady-state visual evoked potentials *IEEE Trans. Neural Syst. Rehabil. Eng.* **14** 225–9

[53] Wold S, Ruhe A, Wold H and Dunn W J III. 1984 The collinearity problem in linear regression. The partial least squares (PLS) approach to generalized inverses *SIAM J. Sci. Stat. Comput.* **5** 735–43

[54] Eliseyev A and Aksenova T 2013 Recursive N-way partial least squares for brain-computer interface *PLoS One* **8** e69962

[55] Dayal B S and MacGregor J F 1997 Recursive exponentially weighted PLS and its applications to adaptive control and prediction *J. Process Control* **7** 169–79

[56] Dayal B S and MacGregor J F 1997 Improved PLS algorithms *J. Chemom.* **11** 73–85

[57] Bro R 1998 Multi-way analysis in the food industry: models, algorithms, and applications. *Academish proefschrift*. Dinamarca

[58] Bro R 1996 Multiway calibration. *Multilinear PLS J. Chemom.* **10** 47–61

[59] Cichocki A et al 2015 Tensor decompositions for signal processing applications: from two-way to multiway component analysis *IEEE Signal Process. Mag.* **32** 145–63

[60] Bishop C M 2006 Introduction *Pattern Recognition and Machine Learning* (New York: Springer)

[61] Schaeffer M-C 2017 ECoG signal processing for brain computer interface with multiple degrees of freedom for clinical application Université Grenoble Alpes

[62] Sutton C and McCallum A 2012 An introduction to conditional random fields *Found. Trends Mach. Learn.* **4** 267–373

[63] Rabiner L R 1989 A tutorial on hidden Markov models and selected applications in speech recognition *Proc. IEEE* **77** 257–86

[64] Brain computer interface: neuroprosthetic control of a motorized exoskeleton clinicaltrials.gov (available at: https://clinicaltrials.gov/ct2/show/NCT02550522)

[65] ICTRP clinical trial NCT02550522 (available at: https://apps.who.int/trialsearch/Trial2.aspx?TrialID=NCT02550522)

[66] Larzabal C, Bonnet S, Costecalde T, Auboiroux V, Charvet G, Chabardes S, Aksenova T and Sauter-Starace F 2021

Long-term stability of the chronic epidural wireless recorder WIMAGINE in tetraplegic patients *J. Neural Eng.* **18** 056026

[67] Robinet S, Audebert P, Regis G, Zongo B, Beche J F, Condemine C and Charvet G 2011 A low-power 0.7 muV 32-channel mixed-signal circuit for ECoG recording *IEEE J. Emerging Sel. Top. Circuits Syst.* **1** 451–60

[68] Morinière B, Verney A, Abroug N, Garrec P and Perrot Y 2015 EMY: a dual arm exoskeleton dedicated to the evaluation of brain machine interface in clinical trials 2015 *IEEE/RSJ Int. Conf. on Intelligent Robots and Systems (IROS)* pp 5333–8

[69] Nguyen C H, Karavas G K and Artemiadis P 2019 Adaptive multi-degree of freedom brain computer interface using online feedback: towards novel methods and metrics of mutual adaptation between humans and machines for BCI *PLoS One* **14** e0212620

[70] Vidaurre C, Schlogl A, Cabeza R, Scherer R and Pfurtscheller G 2006 A fully on-line adaptive BCI *IEEE Trans. Biomed. Eng.* **53** 1214–9

[71] Waldert S et al 2009 A review on directional information in neural signals for brain-machine interfaces *J. Physiol.* **103** 244–54

[72] Schalk G et al 2007 Decoding two-dimensional movement trajectories using electrocorticographic signals in humans *J. Neural Eng.* **4** 264–75

[73] Volkova K, Lebedev M A, Kaplan A and Ossadtchi A 2019 Decoding movement from electrocorticographic activity: a review *Front. Neuroinf.* **13** 74

[74] Nakanishi Y et al 2017 Mapping ECoG channel contributions to trajectory and muscle activity prediction in human sensorimotor cortex *Sci. Rep.* **7** 1–13

[75] Young D et al 2019 Closed-loop cortical control of virtual reach and posture using Cartesian and joint velocity commands *J. Neural Eng.* **16** 026011

[76] Fukuma R et al 2015 Closed-loop control of a neuroprosthetic hand by magnetoencephalographic signals *PLoS One* **10** e0131547

[77] Jerbi K et al 2011 Inferring hand movement kinematics from MEG, EEG and intracranial EEG: from brain-machine interfaces to motor rehabilitation *IRBM* **32** 8–18