Repositorio Institucional de Documentos

2014 23 Iñaki Asier Iturrate Gil Robot Learning and Control Using ErrorRelated Cognitive Brain Signals Departamento Director/es Informática e Ingeniería de Sistemas Mínguez Zafra, Javier Montesano del Campo, Luis Director/es Tesis Doctoral Autor Repositorio de la Universidad de Zaragoza – Zaguan http://zaguan.unizar.es UNIVERSIDAD DE ZARAGOZA Departamento Director/es Iñaki Asier Iturrate Gil ROBOT LEARNING AND CONTROL USING ERROR-RELATED COGNITIVE BRAIN SIGNALS Director/es Informática e Ingeniería de Sistemas Mínguez Zafra, Javier Montesano del Campo, Luis Tesis Doctoral Autor 2014 Repositorio de la Universidad de Zaragoza – Zaguan http://zaguan.unizar.es UNIVERSIDAD DE ZARAGOZA Departamento Director/es Director/es Tesis Doctoral Autor Repositorio de la Universidad de Zaragoza – Zaguan http://zaguan.unizar.es UNIVERSIDAD DE ZARAGOZA Robot learning and control using error-related cognitive brain signals Iñaki Asier Iturrate Gil Thesis submitted in fulllment of the requirements for the degree of Doctor of Philosophy in Computer Science and Systems Engineering Supervised by Javier Minguez and Luis Montesano Instituto de Investigación en Ingeniería de Aragón (I3A) Departamento de Informática e Ingeniería de Sistemas (DIIS) Escuela de Ingeniería y Arquitectura (EINA) Universidad de Zaragoza November, 2013 ii Acknowledgements First of all, I would like to thank my advisors Javier Minguez and Luis Montesano for their patience, assistance and guidance during the good and the bad moments. If I had to choose only one of the previous three terms, I would surely choose patience. I strongly hope the professional and personal links forged keep going on in the future. Carlos, Edu, Mauricio, my present for you is that you begin this paragraph. This work would not have been possible without my thesis brothers (the best anyone could have) Carlos Escolano, Edu López and Mauricio Antelis. They have been my friends long before I started it, they have always been there during it, and I am pretty sure they will be there afterwards. Jason Omedes, our last vcizar mate, I do not forget you either. Thanks for re-shaping the lab and being such a good guy. I am very happy to have met you. And of course, I wanted to thank all the workers from Bitbrain, especially Dani. María López, thank you for always trusting in me and giving me the opportunity to enter in this amazing world. I am strongly indebted to my advisors in Lausanne José del R. Millán and Ricardo Chavarriaga. Thanks to them, I was able to grow and see things from a dierent point of view. I could not forget thanking at all the teammates from the CNBI, and specially Andrea Biasiucci and Michele Tavella for taking care of me. Fabiola, I could not forget you either. I am also in debt with Manuel Lopes and Jonathan Grizou, which helped me (and are helping right now) in the completion of my thesis. My family has always been there when I needed it for giving me support and their shoulder. Dad, you have always been my example to follow. Mom, you are simply the best. My incredible brothers, Iker (and Maria), Joseba (and Miriam), thank you for all your counsel and caring. To Ramón and Elena, the best parents-in-law one could ever wish. My friends, nothing could be possible without you. TPM (and ladies), this goes for you: Buka and Nahikari (ponme otro que lo he tirado), Gastón (vamos a ver a los ciclos), Ichi and Irene (colectivo viernes), Latorre and Yoli (pirotecnia?), Motilva and Isa (a ver cuándo nos pegamos) and Oscar (y no pocou). Motilva, Ichi, Gastón, por los momentos duros. Popy, Tamara, Miguel, Estefania, César, Ana, you became closer than I would have never imagined. Matthias osaesmamado, mein lieber Deutsch, thank you n . Next time I will throw your shoe further. Romeo, R, Jarmok, decrack, Morfeo, siempre nos quedará Juslibol. Sinu, Pablo, Raul, you had to be here. To the array of pointers to Informáticos. To the Redima F7 team. Thank you all. This thesis is for you, María. You are my light in the dark, I love you. iii iv Resumen Durante los últimos años, el campo de los interfaces cerebro-máquina (BMIs en inglés) ha demostrado cómo humanos y animales son capaces de controlar dispositivos neuroprotésicos directamente de la modulación voluntaria de sus señales cerebrales, tanto en aproximaciones invasivas como no invasivas. Todos estos BMIs comparten un paradigma común, donde el usuario trasmite información relacionada con el control de la neuroprótesis. Esta información se recoge de la actividad cerebral del usuario, para luego ser traducida en comandos de control para el dispositivo. Cuando el dispositivo recibe y ejecuta la orden, el usuario recibe una retroalimentación del rendimiento del sistema, cerrando de esta manera el bucle entre usuario y dispositivo. La mayoría de los BMIs decodican parámetros de control de áreas corticales para generar la secuencia de movimientos para la neuroprótesis. Esta aproximación simula al control motor típico, dado que enlaza la actividad neural con el comportamiento o la ejecución motora. La ejecución motora, sin embargo, es el resultado de la actividad combinada del córtex cerebral, áreas subcorticales y la médula espinal. De hecho, numerosos movimientos complejos, desde la manipulación a andar, se tratan principalmente al nivel de la médula espinal, mientras que las áreas corticales simplemente proveen el punto del espacio a alcanzar y el momento de inicio del movimiento. Esta tesis propone un paradigma BMI alternativo que trata de emular el rol de los niveles subcorticales durante el control motor. El paradigma se basa en señales cerebrales que transportan información cognitiva asociada con procesos de toma de decisiones en movimientos orientados a un objetivo, y cuya implementación de bajo nivel se maneja en niveles subcorticales. A lo largo de la tesis, se presenta el primer paso hacia el desarrollo de este paradigma centrándose en una señal cognitiva especíca relacionada con el procesamiento de errores humano: los potenciales de error (ErrPs) medibles mediante electroencefalograma (EEG). En esta propuesta de paradigma, la neuroprótesis ejecuta activamente una tarea de alcance mientras el usuario simplemente monitoriza el rendimiento del dispositivo mediante la evaluación de la calidad de las acciones ejecutadas por el dispositivo. Estas evaluaciones v vi se traducen (gracias a los ErrPs) en retroalimentación para el dispositivo, el cual las usa en un contexto de aprendizaje por refuerzo para mejorar su comportamiento. Esta tesis demuestra por primera vez este paradigma BMI de enseñanza con doce sujetos en tres experimentos en bucle cerrado concluyendo con la operación de un manipulador robótico real. Como la mayoría de BMIs, el paradigma propuesto requiere una etapa de calibración especíca para cada sujeto y tarea. Esta fase, un proceso que requiere mucho tiempo y extenuante para el usuario, diculta la distribución de los BMIs a aplicaciones fuera del laboratorio. En el caso particular del paradigma propuesto, una fase de calibración para cada tarea es altamente impráctico ya que el tiempo necesario para esta fase se suma al tiempo de aprendizaje de la tarea, retrasando sustancialmente el control nal del dispositivo. Así, sería conveniente poder entrenar clasicadores capaces de funcionar independientemente de la tarea de aprendizaje que se esté ejecutando. Esta tesis analiza desde un punto de vista electrosiológico cómo los potenciales se ven afectados por diferentes tareas ejecutadas por el dispositivo, mostrando cambios principalmente en la latencia la señal; y estudia cómo transferir el clasicador entre tareas de dos maneras: primero, aplicando clasicadores adaptativos del estado del arte, y segundo corrigiendo la latencia entre las señales de dos tareas para poder generalizar entre ambas. Otro reto importante bajo este paradigma viene del tiempo necesario para aprender la tarea. Debido al bajo ratio de información transferida por minuto del BMI, el sistema tiene una pobre escalabilidad: el tiempo de aprendizaje crece exponencialmente con el tamaño del espacio de aprendizaje, y por tanto resulta impráctico obtener el comportamiento motor óptimo mediante aprendizaje por refuerzo. Sin embargo, este problema puede resolverse explotando la estructura de la tarea de aprendizaje. Por ejemplo, si el número de posiciones a alcanzar es discreto se puede pre-calcular la política óptima para cada posible posición. En esta tesis, se muestra cómo se puede usar la estructura de la tarea dentro del paradigma propuesto para reducir enormemente el tiempo de aprendizaje de la tarea (de diez minutos a apenas medio minuto), mejorando enormemente así la escalabilidad del sistema. Finalmente, esta tesis muestra cómo, gracias a las lecciones aprendidas en los descubrimientos anteriores, es posible eliminar completamente la etapa de calibración del paradigma propuesto mediante el aprendizaje no supervisado del clasicador al mismo tiempo que se está ejecutando la tarea. La idea fundamental es calcular un conjunto de clasicadores que sigan las restricciones de la tarea anteriormente usadas, para a continuación seleccionar el mejor clasicador del conjunto. De esta manera, esta tesis presenta un BMI plugand-play que sigue el paradigma propuesto, aprende la tarea y el clasicador y nalmente alcanza la posición del espacio deseada por el usuario. Abbreviations Ne Error-related negativity Pe Error-related positivity ACC Anterior cingulate cortex BCI Brain-computer interface BMI Brain-machine interface CAR Common-average reference EEG Electroencephalography ERN Error-related negativity ERP Event-related potentials ErrP Error-related potentials fMRI Functional magnetic resonance imaging FRN Feedback-related negativity GA Grand average IRL Inverse reinforcement learning LP Low pass RL Reinforcement learning SCI Spinal cord injury SMA Supplementary motor area SNR Signal-to-noise ratio STF Spatio-temporal lter xiii xiv Contents 1|Introduction Brain-machine interfacing (BMI) is an emergent technology developed to provide a communication channel between a human and a device using only brain activity [1]. These systems have been successfully used in four dierent applications: communication and control, motor restoration, entertainment and games, and motor rehabilitation [2]. The most signicant feature that denes a BMI is the method used for the electrophysiological recordings [3]. On one hand, invasive BMIs, such as single-neuron or multiple-neuron recordings, or electrocorticography (ECoG), require surgical procedures as they rely on intracranial electrodes for the signal acquisition, but provide a high quality signal allowing a complex control of dierent devices. On the other hand, non-invasive BMIs, primarily the electroencephalography (EEG), do not require any surgical intervention but usually suer from low signalto-noise ratios. In applications involving motor control, these systems have mainly focused on the operation of neuroprosthetic devices such as mobile robots, robotic wheelchairs or robotic arms [2]. Over the last years, the eld of BMI has witnessed impressive demonstrations of how humans and animals can control the operation of the aforementioned neuroprosthetic devices directly from a voluntary modulation of their brain signals, using both invasive and non-invasive approaches [414]. All these BMIs share a common approach, namely the control paradigm (see Figure 1.1). In this paradigm, the user conveys a variety of information related to the operation of the neuroprosthesis, ranging from continuous velocity/position [4,8,12,1419] and muscular activity [11] to rather discrete states such as directions [5,10,13], targets [6,9], hand opening/closing [4,7], and self-initiation of movements [2022]. This information is recorded from the user's brain activity, and translated into control commands to operate a neuroprosthetic device. Whenever the device executes a command, the user receives a feedback of the system performance, closing this way the loop between user and device [23]. Most of these BMIs decode control parameters from cortical areas in order to generate the sequence of movements for the neuroprosthesis. This 1 2 Chapter 1. Introduction DECODER LOW-LEVEL CONTROLLER Control signals Control commands Brain commands Feedback Figure 1.1: Usual `control BMI' approach: the user delivers mental commands that specify the next state of the neuroprosthesis. approach closely resembles normal motor control in that it links neural activity to motor behavior or to motor execution [24]. Motor control, however, is the result of the combined activity of three dierent levels: the cerebral cortex, the brainstem, and the spinal cord [24]. The lowest level corresponds to the spinal cord, which handles many elements of basic or stereotypical movements, mainly via rhythmic and oscillatory outputs commonly known as central pattern generators (CPGs) [25,26]. Following the spinal cord, the brainstem is mainly responsible for the selection, enhancement or variation of the dierent patterns present on the lower levels [27]; nally, the cerebral cortex represents the highest level, playing a wide range of roles such as sensory-motor transformations, action understanding, and planning and executing of goal-directed and skilled motor tasks [24,28,29]. Therefore, the involvement of and communication between these three levels will vary together with the motor task being executed. In fact, many elements of skilled movements, from manipulation to walking, are mainly handled at the spinal cord level with cortical areas providing just goals and movement onset [24,26]. But, could a BMI also mimic this strategy? As mentioned before, some studies have shown the feasibility to decode such a kind of cognitive information associated to voluntary, goal-directed movements [6,2022,30]. Nevertheless, this approach requires that the intelligent controller, emulating the roles of the lower levels of motor control, knows how to (i) learn optimal motor executions (i.e. optimal trajectories), and (ii) store and execute them as needed. 3 −200 0 200 400 600 800 1000 −4 −3 −2 −1 0 1 2 3 Time (ms) Amplitude (mV) Error Correct Difference Correct Error T T Figure 1.2: (Left) The device objective is to reach the target position (marked with T), and executes either correct or wrong actions while the user assesses them. (Right) Averaged error-related potential generated from the user's assessments, together with their topographic scalp interpolation (red is positive, blue is negative). The dierence average (error minus correct averages) are used to study these potentials to remove those eects common to both potentials. This dierence average is composed of three main components: a negativity at around 250 ms, and two broader positive and negative peaks at 350 ms and 500 ms respectively. The main objective of this thesis is to propose and develop an alternative BMI paradigm to neuroprosthetics that, instead of decoding and executing low level control commands, makes use of higher-level cognitive information related to the task being performed. This way, the intelligent controller of the device could use this data as a way of dealing with the lack of information from subcortical levels. In other words, the paradigm tries to emulate the role of subcortical levels during motor control. However, entirely emulating this role would require extracting a large amount of cognitive signals, each one dealing with a specic piece of information about the task. In this thesis, we present the rst step towards developing this paradigm by focusing on a specic cognitive signal related to the user's error processing. Speci- cally, this thesis studies human brain signals that carry cognitive information associated to decision-making processes that arise during goal-directed move- 4 Chapter 1. Introduction ments, and whose low-level implementation is handled at the subcortical and spinal cord level. In this approach, the neuroprosthesis is actively performing a task, while the user simply monitors the performance of the neuroprosthesis by assessing the quality of its actions. These assessments are extracted from the user's brain, and translated into feedback signals that help the intelligent controller of the neuroprosthesis how to solve the task. Thus, the development of this new paradigm was guided by its two main modules: the extraction and decoding of error signals from the user's brain, and the use of this information as input for an intelligent controller associated to the device. Regarding the cognitive signals used as feedback for the paradigm, this thesis focuses on EEG signals related to the error processing system of the brain. In the eld of cognitive neuroscience, these event-related potentials (ERP) are usually dened as a negative voltage deection elicited in the user's brain around 200 ms after noticing his/her expected outcome diered from the actual outcome. Several works have associated them to the dopaminergic neural system [3135], which is thought to be richly projected into the anterior cingulate cortex (ACC) region of the brain [36]. This cortex region would be in charge of conveying reinforcement learning signals (i.e. rewards) about the ongoing events [31]. In fact, many works have linked the dopamine, reinforcement learning and error processing as a single neural framework [31,32,35,3741]. A large number of works have demonstrated the existence of these signals in very dierent situations, always by averaging across hundreds of trials to increase the EEG signal-to-noise ratio: when a subject performs a choice reaction task under time pressure and realizes that he/she has committed an error [42] (the so-called error-related negativity ERN or Ne , usually followed by a sharp positivity, Pe ); after a user is given feedback about a performed task [43] (the feedback-related negativity, FRN); and when the subject perceives an error committed by another person (observation ERN) [44]. A few years ago, the eld of BMI became interested in these signals [45] nding that a variant of them, the error-related potentials (ErrP), were also present when the user delivered an order to the machine and the machine executed another one [46] (interaction ErrP); or when a user observed a device committing an error [47]. In spite of resembling a similar pattern to the ERN and FRN potentials, the error-related potentials are always followed by a sharp positivity (probably the Pe [42]) and by a second negativity, which could be related to a visual semantic mismatch [48] (see Figure 1.2). Despite the existence on average of these signals in controlled situations, detecting them on single trial poses a harder problem, due to the non-stationary behavior of the EEG and its poor signal-to-noise ratio. Furthermore, the event-related potentials have additional sources of variability 5 that can aect the amplitude or the latency of their components, such as the arousal [49], the probability of occurrence of the expected stimulus [50] or the stimulus evaluation time (i.e., the amount of time required to perceive and categorize a stimulus) [50,51]. To solve this problem, the BCIs need a calibration phase to train a classier prior to the device control. Recently, several works have successfully decoded these signals online, and used them to adapt the BMI classier [5254] or to prevent the execution of missclassi- ed commands [5560]. In this thesis, once these signals are decoded thanks to the calibration phase, they can be used as feedback for an intelligent controller. Given that these signals are allegedly related to a human reinforcement learning model, the eld of reinforcement learning (RL) [61] provides a natural decision for the development of this module. In reinforcement learning, the goal of the device is to learn a policy mapping from situations to actions that maximize the expected reward [62]. Since its beginnings over twenty years ago, many approaches have been proposed and demonstrated over a wide range of applications, either in discrete spaces [61], continuous spaces under virtual scenarios [6366], and for learning complex tasks in robotics [67,68]. Throughout this thesis, the reinforcement learning module played an important role for the development of this paradigm, closing the BMI loop during the device control. However, for this rst attempt of the proposed BMI paradigm, the closed loop worked in a dierent way: rather than sending commands and receving feedback from the device as in the control paradigm, the proposed system sends feedback in both directions (user to device, and device to user), allowing for a constant co-adaptation between the two. Once the paradigm is dened in a general context, it presents the rst challenge to resolve: it is necessary to prove the advantages of the paradigm in realistic scenarios, such as the real-time use of a neuroprosthetic device. In this thesis, we show how can be achieve this objective. A second challenge comes from the fact that, under this paradigm, the device needs a certain amount of time to learn the task being executed. Despite it would be desirable to have a system that is able to learn from the beginning of the experiment, in practice it is also needed for each task a calibration phase to train and detect online the EEG signals used, a timeconsuming operation that adds to the task learning time and hinders the deployment of BMIs out of the lab. In the specic case of our paradigm, it would be thus convenient to train classiers able to work irrespectively of the learning task being performed. Many works have tried to reduce this calibration time in several ways: Adapting the classier with supervised techniques incorporating labeled examples of subsequent sessions [69,70] or with unsupervised techniques [71,72]; initializing the classier model with 6 Chapter 1. Introduction data from a pool of subjects [73, 74]; or nding time-invariant features of the EEG to control the BCI [75, 76]. However, it is still unclear whether (and how) the error-related potentials vary across dierent tasks, and how to create classiers able to generalize among dierent tasks. In this thesis, we study transfer learning across dierent tasks as a way of reducing the calibration time. The third important challenge of the proposed paradigm is related to scalability. Firstly, it is impractical to learn the behavior each time a task has to be performed. Therefore, once optimal motor behaviors are learned, it is necessary to store and execute them on demand similarly to human motor control. In this way, learned behaviors are reused and can be part of more complex behaviors. Indeed, learning more complex behaviors also pose several problems for the proposed paradigm since BMIs have in general low information transfer rates. As the learning time grows exponentially as the task space increases, learning the optimal motor behavior via RL using brain signals may require too much time and eort from the user. Instead, it would be useful to pre-compute optimal motor behaviors rather than learning them from the user's feedback. This pre-computation can be done by exploiting the learning task structure (e.g. using external sensors to analyze the environment and localize tasks/objects of interest). The policies can be computed from scratch or adapted from previously learned policies taking into account the current context. In this thesis, this concept is illustrated using a reaching task with a discrete number of possible nal positions. The corresponding optimal policies are computed oine, assuming that an optimal policy is the one that reaches the target following the shortest path. Similar strategies have been exploited by the so-called shared-control systems, where the device does not only execute the decoded commands, but is also involved in performing the task [2] (e.g. by taking into account the environment while reaching a target or avoiding an obstacle [9,10]). The thesis shows how it is possible to store a collection of motor behaviors, and choose among them via the error potentials the one that most closely resembles the user expectations. In other words, the user can exploit previously learned policies to select which one has to be executed. An alternative interpretation is that by reducing the task space to a set of predened policies, the time needed to learn and execute the task is drastically reduced (from tens of minutes to even half a minute) and, more importantly, it improves the scalability of the system to larger and more complex tasks. Finally, given a suciently informative task structure, it could be possible not only to learn the task, but also to learn online and unsupervisedly the classier while the task is being executed [74,77]. This way we will show how, using all the lessons learned throughout the thesis, we can have a plug- 7 and-play BMI that follows the described paradigm, learns the task and the classier and nally reaches the desired position by the user. 8 Chapter 1. Introduction 1.1 Structure and publications The contents of the thesis are organized as follows. Chapter 2 focuses on the rst objective of the thesis: the development of an alternative paradigm for BMIs based on reinforcement learning and brain-decoded reward signals. Specically, this chapter demonstrates how it is possible to learn optimal motor executions using users' brain feedback. This chapter was presented in [7881]. Chapter 3 performs a deep analysis of the brain signals used for the paradigm: the error-related potentials. In this chapter, we demonstrate how the error-related potentials vary from one experimental protocol to another, and how this impedes the classier generalization among dierent tasks. This work has been presented in [82,83]. Chapter 4 shows how to reduce the calibration time of a given ERP experimental protocol by compensating the ERP latency variations across experiments so as to re-use data from previous experiments. This analysis was performed in two dierent ERPs: the error-related potentials and the P300 potentials [84]; and presented in [85,86]. Chapter 5 extends and improves the paradigm presented in chapter 2 exploiting the task structure by choosing among a set of optimal motor behaviors given the user feedback. These ndings were presented in [8789]. Chapter 6 presents a novel way of removing the calibration phase for the proposed paradigm by estimating at the same time the optimal motor behaviors and the EEG meanings, also exploiting the task structure. This work is under review [90]. Finally, chapter 7 presents the conclusions of this thesis and summarizes some points of future work. 2.2. Methods 15 was applied. Then, eight fronto-central channels (Fz, FCz, Cz, CPz, FC1, FC2, C1, and C2) within a time window of [200,800] ms were downsampled to 64 Hz and concatenated to form a vector x of 312 features per trial. These vectors were normalized and decorrelated using principal component analysis (PCA). A feature selection process based on the r2 score then retained the f -most discriminant features, using a ve-ten-fold cross validation. The f features of all labeled trials were used to train a linear classier (linear discriminant analysis, LDA), which mapped the input feature vector xf of a given trial into a binary output y∈ {-1,+1}. On average, 36 features were selected across all subjects ( ± 13). We report the online single-trial accuracies during the reaching task. To assess the statistical signicance of the obtained accuracies, we compute the chance levels according to the available number of trials using the binomial cumulative distribution [30]. The estimated chance levels at α= 0.05 were 56% for Experiment 1 and 54% for Experiments 2 and 3. 2.2.5 Reinforcement Learning with ErrPs The RL strategy [61] was modeled by a Markov decision process, denoted by the tuple {S, A, r, γ} with S the state space (the possible positions of the device), and A the action space (the possible actions of the device). The reward function r:S×A→ R represented the goodness of the executed action at a given state. The goal of RL was to obtain a policy π:S→A mapping the state space into the action space (i.e., which action had to be performed at each state) so as to maximize the expected return R= P∞ k=0 γkrk+1 at time k . The RL implementation was the Q-learning iterative algorithm [61]: Qk+1(sk, ak) = Qk(sk, ak)+ αhrk+1(sk, ak) + γmax a0∈AQk(sk+1, a0)−Qk(sk, ak)i where k is the current step, γ is a discount factor, and α is the learning rate (for the designed experiments, γ and α were set empirically to 0.4 and 0.1, respectively). During the iterative process at time k , the device executed an action ak from state sk to state sk+1 , receiving a reward rk+1(sk, ak) . This reward was used to update the reinforcement learning policy after each action. All the Q-values were set to zero at the beginning of each run ( k= 0 ). At the end of the run, the nal policy π was computed as the policy that, at each state s , always followed the action a0 with the maximum Q-value, π= arg maxa0∈AQπ(s, a0) . An ε -greedy policy was used to select the next action ak to be executed at each step k of the iterative process. This policy selected the action with 16 Chapter 2. Reinforcement learning using brain signals highest Q-value (best action) for ( 100 −ε) % of the times in the state sk , while a random action was selected the remaining times. The experiments started with a completely exploratory behavior ( ε= 100%), and every time an exploratory action was chosen ε was decreased by a small factor (5%) until reaching a minimum value (20%) to always maintain a small percentage of exploration. The output of the ErrP classier at each step k was used as the reward function rk+1 : after an action was executed, the reward was -1 or +1 depending on whether or not an ErrP was detected. 2.3 Results Figure 2.2A summarizes the three experiments executed by 12 subjects. Subjects' task was to monitor the performance of a neuroprosthetic device that had to reach a target. This target was known by the subject, but not by the neuroprosthesis. Online operation of the closed-loop system was tested on several runs where the goal is to repeatedly reach a xed target location from dierent starting points. At the beginning of each run the device controller was initialized to a random behavior (i.e., equiprobable actions for all states). After each action the controller was updated by means of RL based on the online decoding of the ErrPs. Each of the runs had a dierent target location which remained constant during its entire length (100 device actions). Whenever the device reached the target, it was randomly reset to a new location. On average (for all runs and targets), subjects reached the target 12, 13 and 13 times for experiments 1, 2 and 3, respectively. 2.3.1 Analysis of error-related potentials We analyzed variations in the grand average EEG potentials for both conditions (see Fig. 2B and Fig. 2.3). To this end we performed a statistical analysis on the dierence ERP (error minus correct condition) for all electrodes in the time window [-200, 1000] ms, t= 0 ms being the instant when the device starts to move. Only signals from the training runs having a constant error-rate (20%) are used in this analysis. These runs yielded between 200 and 400 trials for each subject. A 3 (brain area: frontal, central or centro-parietal electrode locations) x 3 (left, midline or right locations) x 3 (experiments) within-subjects ANOVA was performed on the peak amplitudes and latencies of the dierence average [50]. Each group is summarized in Table 2.1. When needed, the GeisserGreenhouse correction was applied to assure sphericity. Pairwise post-hoc 2.3. Results 17 EXPERIMENT 1 293 ms 484 ms B A 352 ms 551 ms EXPERIMENT 2 380 ms 598 ms EXPERIMENT 3 −3 0 0 5 Figure 2.3: Topographic interpolation of the most prominent (A) positive peak and (B) negative peak of the dierence average for each experiment, together with the time point of each peak (in milliseconds). 18 Chapter 2. Reinforcement learning using brain signals tests with the Bonferroni correction were computed to determine the dierences between pairs of experiments. Table 2.1: Electrode locations used as factors for the ERP statistical analysis. Each row corresponds to the brain areas, while columns correspond to the laterality. Left Midline Right Frontal FC3, FC1 Fz, FCz FC2, FC4 Central C3, C1 Cz C2, C4 Centro-parietal CP3, CP1 CPz CP2, CP4 We mainly found signicant eects on the latency but not the amplitude of the dierence potential. The type of experiment signicantly aected the latencies of both the positive ( F2,22 = 41.594, p = 3 ×10−8 ) and negative peaks of the dierence ERP ( F2,22 = 7.522, p = 0.003 ). The brain area also aected the latencies of the positive peak ( F1.32,14.55 = 14.175, p = 0.001 ) but not the negative one F2,22 = 0.911, p = 0.417 ). Similarly, the hemisphere aected the latency of the positive peak ( F2,22 = 5.279, p = 0.013 ), but not the latency of the negative one ( F2,22 = 1.711, p = 0.204 ). No signicant interactions were found. For the positive peak latencies, post-hoc pairwise tests revealed significant dierences between experiments 1 and 2 ( p= 0.0001 ), and between experiments 1 and 3 ( p= 0.0001 ), but not between experiments 2 and 3 ( p= 0.068 ). For the negative peak latencies, there were signicant dierences between experiments 1 and 3 ( p= 0.009 ), but not between experiments 1 and 2 ( p= 0.357 ) nor experiments 2 and 3 ( p= 0.122 ). In contrast, the amplitude of the positive and negative peaks were not signicantly aected by the experiment ( F2,22 = 0.124, p = 0.884 and F2,22 = 2.304, p = 0.123 , respectively) nor the brain area ( F1.19,13.08 = 1.227, p = 0.737 and F1.32,14.47 = 0.071, p = 0.857 ). The laterality signicantly aected the positive peak amplitude ( F2,22 = 4.556, p = 0.022 ), but not the negative peak ( F2,22 = 3.425, p = 0.051 ). As in the case of the peak latencies, no signicant interactions were found. 2.3.2 Analysis of ocular artifacts We assessed the possibility of EEG signal contamination by movement-related saccades. In this study, we compute the grand average ERPs (correct and 2.3. Results 19 EXPERIMENT 1EXPERIMENT 2EXPERIMENT 3 CORRECT ASSESSMENT FC3 FC4Fz −200 0 200 400 600 800 1000 −5 0 5 Time (ms) Amplitude (µV) Left Right −200 0 200 400 600 800 1000 −5 0 5 Time (ms) Amplitude (µV) Left Right −200 0 200 400 600 800 1000 −5 0 5 Time (ms) Amplitude (µV) Left Right −200 0 200 400 600 800 1000 −5 0 5 Time (ms) Amplitude (µV) Left Right Up Down −200 0 200 400 600 800 1000 −5 0 5 Time (ms) Amplitude (µV) Left Right Up Down −200 0 200 400 600 800 1000 −5 0 5 Time (ms) Amplitude (µV) Left Right Up Down −200 0 200 400 600 800 1000 −5 0 5 Time (ms) Amplitude (µV) Left Right Up Down −200 0 200 400 600 800 1000 −5 0 5 Time (ms) Amplitude (µV) Left Right Up Down −200 0 200 400 600 800 1000 −5 0 5 Time (ms) Amplitude (µV) Left Right Up Down FC3 Fz FC4 ERROR ASSESSMENT FC3 FC4 −200 0 200 400 600 800 1000 −5 0 5 Time (ms) Amplitude (µV) Left Right −200 0 200 400 600 800 1000 −5 0 5 Time (ms) Amplitude (µV) Left Right −200 0 200 400 600 800 1000 −5 0 5 Time (ms) Amplitude (µV) Left Right Fz −200 0 200 400 600 800 1000 −5 0 5 Time (ms) Amplitude (µV) Left Right Up Down −200 0 200 400 600 800 1000 −5 0 5 Time (ms) Amplitude (µV) Left Right Up Down −200 0 200 400 600 800 1000 −5 0 5 Time (ms) Amplitude (µV) Left Right Up Down −200 0 200 400 600 800 1000 −5 0 5 Time (ms) Amplitude (µV) Left Right Up Down −200 0 200 400 600 800 1000 −5 0 5 Time (ms) Amplitude (µV) Left Right Up Down −200 0 200 400 600 800 1000 −5 0 5 Time (ms) Amplitude (µV) Left Right Up Down FC3 Fz FC4 EXPERIMENT 1EXPERIMENT 2EXPERIMENT 3 Figure 2.4: Grand averages in channels FC3, Fz, and FC4 for correct and error assessments (left and right respectively) of each movement direction (left, right, up or down), averaged for all subjects ( N= 12 ). Rows correspond to each experiment. The results show that the averaged signals of the same assessment (either correct or error) are very similar across all directions, reducing the possibility of having a systematic inuence on the ErrP classication process. 20 Chapter 2. Reinforcement learning using brain signals error) of all channels separately for each dierent action (moving left, right, up, or down). No substantial dierences were found among these ERPs, suggesting little inuence of eye movements. This is illustrated in Figure 2.4 that show the averages of the three most frontal electrodes: FC3, Fz and FC4. As can be seen, the dierences among assessments (correct or error) were larger than the dierences among directions (left, right, up or down). This is consistent with previous studies that found no inuence of this type for experiment 1 [46,47]. To evaluate the existence of statistical dierences due to both assessments and movement directions, we performed 2 (factor assessment: error or correct) x 4 (factor movement direction: left, right, up or down) within-subjects ANOVAs on the values of the most prominent positive and negative peak amplitudes of the grand averages (note that for experiment 1 the ANOVA was 2 x 2 since there were only two possible movement directions). When needed, the Geisser-Greenhouse correction was applied to assure sphericity. The assessment and direction main eects and the assessment x direction interaction were studied. Regarding the main eects, statistical dierences were found for the assessment for all the experiments, for the positive ( F1,11 = 17.277, p = 0.002 , F1,11 = 15.567, p = 0.002 , and F1,11 = 14.202, p = 0.003 for experiments 1 to 3) and negative ( F1,11 = 10.087, p = 0.009 , F1,11 = 14.658, p = 0.003 , and F1,11 = 11.581, p = 0.006 ) peaks. On the contrary, no signicant differences were found for the direction main eect ( p > 0.1 ). Regarding the assessment x direction interaction, signicant dierences were found during experiment 2 ( F3,33 = 3.721, p = 0.02 and F3,33 = 3.903, p = 0.02 for the positive and negative peak); and during experiment 3 for the negative peak ( F3,33 = 3.461, p = 0.03 ) but not for the rest of the cases ( p > 0.35 ). These results indicated that the largest dierences on the potentials were due to the dierent assessments (error / correct), whereas the movement directions of the device aected less the potentials. To further discard the inuence of artifacts on the assessments, data used to train the classier included all possible movements for each class, thus reducing the possibility that classication was biased by their directions. For instance, during experiment 1, both targets are used for the classier training, thus the error and correct assessments are not likely to be correlated with left or right eye movements. Moreover, results obtained in the generalization test for experiments 2 and 3 further support the fact that classication depends on the movement evaluation and not on its direction. Indeed, the training set contained samples where the target locations were Up and Down, while the BMI was tested on targets Left and Right. Finally, to test whether the trained classiers discriminated dierent directions rather than assessments, 2.3. Results 21 we computed for each subject the accuracy of decoding the dierent pairs of movement directions (e.g. left versus right, up versus left, ...) from a xed assessment (either correct or erroneous) with the same features and classier used during the experiments. The mean accuracies obtained were of 52.16±5.22 , 50.07±5.07 , and 49.48±6.41 for experiments 1 to 3, and thus did not reach the chance levels of 56%, 54% and 54% (see `ErrP classier' in Materials and Methods), indicating that the classier was not trained to distinguish saccadic artifacts, but user's assessments. 2.3.3 Event-related potentials and single-trial classication Figure 2.2B shows the grand average potential on electrode FCz elicited by the user's assessment of both types of movements (correct and erroneous), as well as the dierence ERP (error minus correct). As reported by previous studies [47,78], the dierence potential in the three experiments is mainly composed of two prominent positive and negative peaks at around 300 and 500 ms, respectively. A within-subject statistical analysis shows signicant eects of the experiments on the latency but not the amplitude of the dierence ERPs (see Supplementary Information). Overall, this suggests that these signals reect a common phenomenon (i.e., error-processing) across the three experiments, where the complexity of the experimental protocol mainly aects the temporal characteristics of the brain response. Online single-trial decoding of the ErrPs (error vs. correct) is shown in Fig. 2.2C. Classier accuracy was comparable for all experiments. Performance was above chance level (except for one subject in experiment 2) a necessary condition when using stochastic rewards for reinforcement learning [61]. This demonstrates the feasibility of extracting useful reward-related signals from human brain activity. 2.3.4 Reaching task: Learning to reach the training targets Figure 2.5A summarizes the reaching performance when the closed-loop system is tested on the same target locations used for training the ErrP classier. It shows the optimal action per state (arrow) and the number of subjects that successfully learned those actions (dark red corresponds to all subjects). In experiment 1, the optimal actions were learned in 94% of the cases indicating that the users were able to teach the device a quasi-optimal control strategy. 22 Chapter 2. Reinforcement learning using brain signals EXPERIMENT 1 EXPERIMENT 3 EXPERIMENT 2 291,48 mm 0 50 100 0 20 40 60 80 100 Number of actions executed Actions correctly learned (%) 0 50 100 0 20 40 60 80 100 Number of actions executed Actions correctly learned (%) 0 50 100 0 20 40 60 80 100 Number of actions executed Actions correctly learned (%) TT 0 50 100 0 20 40 60 80 100 Number of actions executed Actions correctly learned (%) T T 0 50 100 0 20 40 60 80 100 Number of actions executed Actions correctly learned (%) TT 0 50 100 0 20 40 60 80 100 Number of actions executed Actions correctly learned (%) 0 50 100 0 20 40 60 80 100 Number of actions executed Actions correctly learned (%) T T 0 50 100 0 20 40 60 80 100 Number of actions executed Actions correctly learned (%) 0 50 100 0 20 40 60 80 100 Number of actions executed Actions correctly learned (%) T T 0 50 100 0 20 40 60 80 100 Number of actions executed Actions correctly learned (%) B 12 11 10 9 8 7 6 V 12 11 10 9 8 7 6 V 12 11 10 9 8 7 6 V A Figure 2.5: Online reinforcement learning. [A] Reaching the same targets used to train the ErrP classier (Left and Right locations for experiment 1, and Up and Down for experiments 2 and 3). [B] Reaching new targets (Left and Right for experiments 2 and 3). For each experiment the top rows in [A] and [B] show, for each state, the optimal action (arrows), and the number of subjects that correctly learned that action (color coded, dark red corresponds to all subjects). Target locations are denoted by green squares. The bottom row shows, in red traces, the percentage of correctly performed actions (mean ± SEM, straight ± dashed lines, N= 12 subjects) as a function of the number of performed actions within a run (100 actions, starting from a random device controller). Blue traces correspond to chance performance level calculated using random rewards (100 repetitions). [C] Percentage of correctly performed actions (y-axis) versus the online ErrP decoding accuracy (x-axis) at the end of each run, for the training (red circles) and new targets (green squares). The straight lines show the best-tting regression line for training (red) and new (green) targets, and each point in the plot corresponds to the mean accuracy of one subject. 2.3. Results 23 In experiment 2, where a more complex state-action space has to be learned, the performance decreased. Nonetheless, 85% of the trajectory motion was learned. Interestingly, no major dierences were found between the performance in experiments 2 and 3. The correct actions were learned for 87% of the states in experiment 3. For all experiments and runs we observe a steady acquisition of a suitable control policy. As shown in Fig. 2.5A(bottom), the device quickly learns the correct action after visiting a state (see Methods). For experiment 1, the number of actions correctly learned was always signicantly above chance level (one-tailed unpaired t-tests, p < 0.05 ). The number of actions correctly learned consistently increased as more actions were performed (correlation between time and actions learned of r= 0.70, p < 1×10−8 ). For experiment 2 and 3, the convergence was slower due to the higher number of states and actions. However, after 9 and 6 actions (for experiment 2 and 3, respectively) the results were signicantly above the chance level (one-tailed unpaired ttests, p < 0.05 ). In both experiments performance does not appear to have reached asymptotic values after 100 actions. This suggests that given more time, the device would have been able to learn the entire policy. There was also a high correlation between time and learned actions ( r= 0.79, p < 1×10−8 for both experiments). The number of actions correctly learned was signicantly correlated with the mean performance of the ErrP classier (Fig. 2.5C). In particular, for experiments 2 and 3 ( r= 0.75, p = 2×10−5 and r= 0.76, p = 1×10−5 , respectively). A lower correlation was found in experiment 1 ( r= 0.13, p = 0.54 ), indicating that learning the 1D policy is less sensitive to misclassication of the error-related signal. 2.3.5 Reaching task: Learning to reach new targets The system was also able to rapidly learn suitable policies to reach target locations not used during the training phase (Fig. 2.5B). For both experiments 2 and 3, no signicant dierences in the number of correctly performed actions were observed with respect to the training targets (paired t-test, p= 0.26 and p= 0.65 , respectively). In experiments 2 and 3, 81% and 85% of the actions were learned, respectively. This shows that the ErrP does not depend on targets, as the ErrP classier does not need to be retrained for unseen targets. As before, the number of actions correctly learned across time for the new targets (Fig 2.5B) was above chance level after the rst movements (onetailed unpaired t-tests, p < 0.05 ), and increased along the entire run. Furthermore, the correlation between time and the number of correctly learned 24 Chapter 2. Reinforcement learning using brain signals actions was also very high ( r= 0.70, p < 1×10−8 for both experiments). When comparing the acquisition of the control policy across time between training and new targets, the learning rate was signicantly equal for more than 19 and 8 actions for experiments 2 and 3, respectively (two-tailed unpaired t-tests, p > 0.05 ). As before, the number of actions correctly learned across time for the new targets (Fig 2.5B) was above chance level after the rst movements (onetailed unpaired t-tests, p < 0.05 ), and increased along the entire run. Furthermore, the correlation between time and the number of correctly learned actions was also very high ( r= 0.74, p = 3 ×10−5 and r= 0.78, p = 8 ×10−6 for experiments 2 and 3, respectively). When comparing the acquisition of the control policy across time between training and new targets, the learning rate was signicantly equal for more than 19 and 8 actions for experiments 2 and 3, respectively (two-tailed unpaired t-tests, p > 0.05 ). There was also a high correlation between the mean classier accuracy and the number of actions correctly learned (Fig. 2.5C) ( r= 0.74, p < 1×10−5 and r= 0.78, p < 1×10−5 for experiments 2 and 3, respectively), with no signicant dierences between the accuracies of training and new targets (paired t-test, p= 0.90 and p= 0.85 for experiments 2 and 3, respectively). 2.4 Discussion This paper describes an alternative and complementary BMI paradigm to neuroprosthetics that decodes cognitive brain signals associated to decisionmaking processes relevant for achieving goals. This new approach is demonstrated with a BMI where the subject monitors the robot actions, and it decodes the brain correlate of the subject's assessment of those decisions (erroneous or correct). This cognitive signal is exploited by an RL algorithm to learn which is the sequence of actions that solves a task. In the experiments reported here, it is assumed that the neuroprosthesis owner wishes to initiate a voluntary, goal-directed movement whose lowlevel execution is delegated to subcortical, spinal cord and musculoskeletal structures. In our case, this lower level of motor control is emulated by an intelligent controller. But, where do the movement goal and onset come from? This kind of information can be decoded from the subject's cortical activity [6,2022]. As a result, the combination of all these sorts of cognitive brain signals would be sucient for the operation of any neuroprosthesis, no matter its complexity and number of degrees of freedom. Indeed, as the experimental results reported in this paper illustrate, the burden of learning the (quasi) optimal trajectories is on the controller side rather than on the 3.2. Methods 31 3.2.3 Analysis of error-related potentials and their taskdependent variations The denition of the observation error-related potentials encompasses the appearance of three main and distinct components on the dierence (error minus correct) average: an N2, a P3, and an N4 component [46]. Regarding the N2 and P3 components, several studies suggest that they may actually be the error-related negativity (ERN or N e ) and the following positivity (P e ) [42], but there is still an open discussion about it [111]. Regarding the N4 component, previous studies have suggested that its generation could be due to a visual semantic mismatch [48]. Nonetheless, for the observation error-related potentials, the three components are originated in the anterior cingulate cortex (ACC, Brodmann Areas 24 and 32) and the pre-supplementary motor area (pre-SMA, Brodmann area 6) [46], suggesting an activation of an errorprocessing system on the brain [32]. The electrophysiology of the error-related potentials and their associated variations were studied through the analysis of the raw EEG, and with a ltered EEG where those components not originated in the brain sources involved in the error-related potentials were removed. The lter eliminated several types of artefacts, including electromyographic activity (such as that provoked by scalp and neck muscles), ocular activity (such as eye movements and blinks) and brain activity not originated within the error processing brain sources [112] (such as spatial attention components [113]). The lter is constituted by two main steps: (i) application of independent component analysis (ICA) [114]; and (ii) isolation of the independent brain components related to error processing, with a posterior re-projection of this information to the sensor space. Note that while ICA techniques have been widely used for the characterization of brain sources [115117] and the removal of artefacts [118], the lter proposed herein focused on the isolation of the brain process of interest (see [119] for a similar approach for P300 classication). The ICA spatial lter is a statistical model dened as x=As , where x are the input data, and A and s are the mixing matrix and independent components estimated by maximizing the temporal independence among the components. Each column vector ai of A is the spatial pattern associated with the component si . While there are many ways to compute the ICA model [114], its computation has two diculties in the EEG context: the number of independent components to estimate [117] and the non-reliable nature of the estimation process [120]. The number of components d was estimated using m -fold cross-validation principal component analysis (PCA) in the sensor space [121] (the number of folds was xed to m= 5 in the experiments, leading to d≃15 dimensions retained from the original 32 32 Chapter 3. ErrP task-dependent variations channels). The ICASSO technique was used to address the non-reliability of ICA [120]. ICASSO estimates N ICA lters using the FastICA algorithm [114] under changes of the initial conditions, and then performs clustering on the obtained estimations ( N was xed empirically to 100 ). Once the ICA model was computed, the DipFit source localization [117] was used to estimate the neural source of each component. Those independent components whose brain source was in the ACC or the pre-SMA were selected, as these areas are believed to be the main generators of error-processing brain activity [32,42,116]. For each subject and task/subtask, the number of components selected was between one and four, which were re-projected back to the sensor space to obtain the ltered data. The analysis of the shape and timing of the potentials (with and without the lter) for each task and subtask was carried out through the computation of time-locked averaged potentials for the error and non-error potentials in channel FCz, through the dierence average (error minus non-error averages) [46,79] and by an r2 discriminability test [1]. A topographic interpolation of the potentials was obtained at the time of the main peaks of the dierence average. A source localization analysis was also performed with sLoreta [122] at the N4 component of the error averages of each operational task [46]. Additionally, the peak amplitudes and latencies of the most prominent negativity were extracted from single-trial signals as the minimum value within the time window [320,600] ms in channel FCz. The latency-sorted single-trial potentials were plotted as a colour-encoded image with a smoothing window of 50 trials. Finally, to assess the statistical dierences among tasks/subtasks of the error-related potentials, one-way within-subjects ANOVAs (factor: tasks or subtasks) were conducted over the latencies and amplitudes of each component (N2, P3 and N4) of the dierence average, averaged from channels Fz, FCz and Cz [50,123]. When needed, the Geisser-Greenhouse correction was applied to the data to assure sphericity. 3.2.4 Analysis of the impact of task-dependent signal variations Once the existence of signal variations was studied, its impact on BCI performance was analyzed at two levels: (i) changes in the features distributions used for error detection and (ii) the corresponding classication accuracy. This analysis was performed using the ltered EEG data to avoid the in- uence of activity not originated in the error-processing brain areas (i.e., artefacts). 3.2. Methods 33 Feature extraction Previous studies have demonstrated that amplitude values of the error-related potentials from several fronto-central channels are suitable features for their discrimination (error vs non-error) [46, 47, 79]. In this study, features are constructed as linear combinations of amplitudes of channels and time points that best separate these two classes [80]. Given a set of n labelled trials of the two classes, for each trial, eight fronto-central channels (Fz, FC1, FCz, FC2, C1, Cz, C2, and CPz) within a time window of [200,800] ms were subsampled to 64 Hz and concatenated as a vector of 312 features. The feature vectors of all trials were normalized, and then decorrelated using PCA, retaining 95% of the explained variance. The k -most discriminant features were selected based on a robust variant of the Fisher score [72]: FS(fi) = |med(fi 1)−med(fi 2)| medad(fi 1) + medad(fi 2), (3.1) with med(fi j) and medad(fi j) being the median and the median absolute deviation of feature fi for class j∈ {1,2} . The number of features to retain k was determined by a ten-fold cross validation. The eect of signal variations was measured by the statistical signicance between the features distributions for each class, between operational tasks (inter-task) and between subtasks (intra-task). One-way within-subjects ANOVAs (factor: tasks or subtasks) were conducted on the single-trial features of each class to assess the statistical signicance. In addition, the inter/intra-task similarity of the features' distributions was quantied by the Kullback-Leibler (KL) divergence. The KL divergence from P∼ N(µP,ΣP) to Q∼ N(µQ,ΣQ) is: DKL(P||Q) = 1 2 tr(Σ−1 QΣP) + vTΣ−1 Qv−ln |ΣP| |ΣQ|−k! (3.2) with v= (µQ−µP) . High values of DKL(P||Q) entail large dierences between distributions. The KL divergences were computed using the k= 10 most discriminant features (according to equation 3.1) to compare the intra/inter features distributions. Single-trial classication The classier used in the analysis was a regularized version of the linear discriminant analysis (LDA) [124]. The LDA discriminant function D(f) is the hyperplane that maximally separates the feature distributions corresponding 34 Chapter 3. ErrP task-dependent variations to two classes: D(f) = wTf+b , where f is the feature vector to be classied, and w and b are the normal vector to the hyperplane and the corresponding bias computed by: ˆµ=1 2(ˆµ1+ ˆµ2) , w=˜ Σ−1(ˆµ2−ˆµ1) , b=−wTˆµ ; where ˆµj is the sample mean of class j , ˆµ the sample global mean, and ˜ Σ the regularized sample covariance matrix (shared for both classes). Regularization aims to minimize the covariance estimation error E= |Σ−˜ Σ| , with Σ being the real covariance matrix, by penalizing very large and very small eigenvalues. The regularized covariance matrix was computed by: ˜ Σ = (1 −γ(n))ˆ Σ + γ(n)νI ; where n was the number of trials used for training, γ(n)∈[0,1] the regularization factor (whose value can be computed numerically [124]), ˆ Σ the sample covariance matrix, and ν=tr(ˆ Σ)/k the average eigenvalue of ˆ Σ , with k being the number of diagonal elements of ˆ Σ . To tackle the signal variations, the classier was adapted based on a sequential process in which labelled examples of the new task were used to modify the discriminant function of the LDA classier [72]. Namely, given a new example of class j at time t , fj(t) , the mean ˆµj was updated using an exponential moving average: ˆµj(t) = (1 −α)ˆµj(t−1) + αfj(t). (3.3) α∈[0,1] is the update parameter (xed to α= 0.05 [72]) and the initial values for ˆµj were obtained using data of another task or subtask. The corresponding discriminant function was recomputed using the new mean ˆµj to update the w and b parameters accordingly. The impact of signal variations was analyzed rstly without adaptation, and then with the supervised adaptation. Firstly, the classier was trained with examples of one task and tested with another one (inter-task) or trained with examples of one subtask and tested with the other subtasks (intratask). The results were then compared with the performance of the baseline classier computed using a ten-fold cross-validation scheme for each task (subtask) separately. Secondly, the adaptation was evaluated against the baseline classier performance as a function of the number of trials used to train/adapt the classier (i.e., calibration time). The adaptive classication results were obtained using the train-test sets, as follows: the classier was initially trained using the train dataset and then, the test dataset was split into two subsets ( D1 with 300 examples and D2 with the remaining). For each trial at time t , the classier was updated using equation 3.3 and tested on the D2 dataset. The results were then compared with the performance of baseline classiers built using trials [1 . . . t] of D1 and tested on D2 . This process was repeated 10 times to reduce variability in the results for adaptive and baseline classiers while shuing trial positions, and then averaging the 3.3. Results 35 obtained accuracies. 3.3 Results 3.3.1 Electrophysiology of potentials their signal variations The analysis comprises the raw EEG data along with the ltered EEG data (see subsection 3.2.3). The proposed lter eliminated 90% of the ICA components that were not estimated within ACC or pre-SMA. The majority of these components were ocular artefacts such as eye movements (estimated in frontal areas such as Brodmann areas 10, 11 or 38) and brain activity estimated either in parieto-occipital/occipital areas (Brodmann areas 17, 18 and 19) or the posterior cingulate cortex (PCC, Brodmann areas 23 and 31). Although these components contributed to the EEG, they were not originated in the main error-processing areas and thus were eliminated by the lter. Artefact correlation with the subtasks but not with the tasks is an eect worthy of mention, which might aect feature extraction and classication analysis of the generalization study. Figure 3.2 displays an example with the raw and ltered EEG for the subtask OT1.Right. Without ltering, a signed r2 discriminability test indicated that the most discriminant features were on frontal channels (originated by lateral eye movements); on the other hand, after ltering, the most discriminant features were due to fronto-central activations (originated by error-related potentials). Note that without ltering, the most discriminant features may greatly discriminate the potentials within the subtask OT1.Right (due to the lateral eye movements), but would not generalize for the other subtasks or task as they involved dierent eye movements. For both operational tasks and subtasks, the average dierence of the raw/ltered potentials for both conditions presented a small negative deection approximately at 250 ms (an N2 component) and prominent positive and negative peaks (P3 and N4 components) at approximately 300 ms and 500 ms (in agreement with the r 2 test, gure 3.3). The topographical scalp maps at these last two peaks showed fronto-central activations for the two operational tasks. When not ltering the data, source estimations for OT1 (at 500 ms of the error average) were in the paracentral lobule (Brodmann Area 5), whereas for OT2 (at 450 ms of the error average) were in the ACC (Brodmann area 24). On the contrary, when ltering the data the potentials from both operational tasks were estimated within the ACC, agreeing with previous studies on 36 Chapter 3. ErrP task-dependent variations Time (ms) Channels 0 100 200 300 400 500 600 700 800 900 1000 32.Oz 31.Pz 30.Cp2 29.Cpz 28.Cp1 27.Cz 26.Fc2 25.Fcz 24.Fc1 23.Fz 22.O2 21.O1 20.P8 19.P4 18.P3 17.P7 16.Cp6 15.Cp5 14.T8 13.C4 12.C3 11.T7 10.Fc6 9.Fc5 8.F8 7.F4 6.F3 5.F7 4.Af4 3.Af3 2.Fp2 1.Fp1 −0.5 −0.4 −0.3 −0.2 −0.1 0 0.1 0.2 0.3 0.4 0.5 Time (ms) Channels 0 100 200 300 400 500 600 700 800 900 1000 32.Oz 31.Pz 30.Cp2 29.Cpz 28.Cp1 27.Cz 26.Fc2 25.Fcz 24.Fc1 23.Fz 22.O2 21.O1 20.P8 19.P4 18.P3 17.P7 16.Cp6 15.Cp5 14.T8 13.C4 12.C3 11.T7 10.Fc6 9.Fc5 8.F8 7.F4 6.F3 5.F7 4.Af4 3.Af3 2.Fp2 1.Fp1 −0.5 −0.4 −0.3 −0.2 −0.1 0 0.1 0.2 0.3 0.4 0.5 RAW EEG SIGNAL FILTERED EEG SIGNAL Frontal Channels Fronto-Central Channels Frontal Channels Fronto-Central Channels Figure 3.2: Signed r2 discriminability test of non-error versus error, performed on the subtask OT1.Right, when not ltering the signal (Left) and when ltering the signal (Right). The x-axis represents the time (from 0 to 1000 ms) and the y-axis represents each recorded EEG channel. Topographic interpolation of the r2 is shown at 350 and 500 ms. The solid boxes mark the position of fronto-central channels, whereas the dashed boxes mark the position of frontal channels. When not ltering the signals, most of the discriminability comes from frontal channels with the sign reversed on the left and right hemispheres. When ltering the signals, most of the discriminability comes from fronto-central channels. 3.3. Results 37 error-related potentials [46,47,78,116]. These results indicated that the use of the ICA lter was advisable for the isolation of the error-processing brain activity. Regarding the inter-task analysis, visual inspection revealed that the shape of the averaged potentials diered between the two operational tasks (Figure 3.3). In OT1, the averaged error-related potentials presented two positive peaks at 200 and 300 ms, whereas in OT2 the rst positive peak was smaller, with a more prominent peak at 280 ms. This dierence was also appreciated in the sorted single-trial error-related potentials. Additionally, the ANOVA analyses reported statistical dierences on the latency of the three main components both for the raw ( F1,9= 9.574, p = 0.01 ; F1,9= 24.469, p = 0.001 ; and F1,9= 48.442, p = 1·10−4 for the N2, P3 and N4 components) and ltered data ( F1,9= 7.789, p = 0.02 ; F1,9= 24.970, p = 0.001 ; and F1,9= 28.809, p = 0.0005 ). For the amplitude of the components, statistical dierences were found only for the N4 component of the ltered data ( F1,9= 15.66, p = 0.003 ). These results indicated the existence of signal variations in the error-related potentials between operational tasks aecting mainly the latency of their main components. Regarding the intra-task analysis, visual inspection revealed that the shape of the averaged potentials was very similar among subtasks (Figure 3.3). The ANOVA analyses reported no statistical dierences ( p > 0.05 ) except for the N2 component of OT1 for the raw data ( F2,18 = 12.021, p = 0.0005 ). These results indicated that on average, the components did not change among subtasks of the same task. 3.3.2 Features analysis Regarding the inter-task analysis, visual inspection of the features showed that only the best feature ( f1 ) reected similar patterns between OT1 and OT2, whereas the other features presented dierent spatio-temporal combinations (see Figure 3.4 Top-Middle for representative examples). In the intra-task case, the features were very similar between subtasks. For instance, the best features ( f1 to f3 in the gure) were almost equal among subtasks while the worst feature f10 presented greater variations. For each class (non-error and error), feature distributions were signi- cantly dierent for interand intra-tasks (ANOVA test, p < 0.001 in all the cases). For the inter-task case, KL divergences were 0.64 ±0.34 and 1.71 ±0.46 for non-error and error respectively, while for the intra-tasks of OT1 and OT2 the divergences were 0.64±0.51 and 1.06±0.28 , and 0.44±0.19 and 1.92 ±0.60 . For all tasks and subtasks, the non-error KL divergences were signicantly lower than the error KL divergences (unpaired one-tailed 38 Chapter 3. ErrP task-dependent variations OPERATIONAL TASK 2 NON-ERROR POTENTIALS ERROR POTENTIALS OPERATIONAL TASK 1 −6 −4 −2 0 2 4 Amplitude (µV) Error Correct Difference Time (ms) −200 0 200 400 600 800 1000 0 0 . 05 0 .1 0 .1 5 −3 −2 −1 0 1 2 3 Best Match: Brodmann Area 5 −6 −4 −2 0 2 4 Amplitude (µV) Error Correct Difference Time (ms) −200 0 200 400 600 800 1000 0 0 . 05 0 .1 0 .1 5 −3 −2 −1 0 1 2 3 Best Match: Brodmann Area 24 −6 −4 −2 0 2 4 Amplitude (µV) Error Correct Difference Time (ms) −200 0 200 400 600 800 1000 0 0 . 05 0 .1 0 .1 5 OT1.UP −6 −4 −2 0 2 4 Amplitude (µV) Error Correct Difference Time (ms) −200 0 200 400 600 800 1000 0 0 . 05 0 .1 0 .1 5 OT1.LEFT −6 −4 −2 0 2 4 Amplitude (µV) Error Correct Difference Time (ms) −200 0 200 400 600 800 1000 0 0 . 05 0 .1 0 .1 5 OT1.RIGHT RAW EEG SIGNAL FILTERED EEG SIGNAL −6 −4 −2 0 2 4 Amplitude (µV) Error Correct Difference Time (ms) −200 0 200 400 600 800 1000 0 0 . 05 0 .1 0 .1 5 OT2.DOWN −6 −4 −2 0 2 4 Amplitude (µV) Error Correct Difference Time (ms) −200 0 200 400 600 800 1000 0 0 . 05 0 .1 0 .1 5 OT2.UP −6 −4 −2 0 2 4 Amplitude (µV) Error Correct Difference Time (ms) −200 0 200 400 600 800 1000 0 0 . 05 0 .1 0 .1 5 OT2.LEFT −6 −4 −2 0 2 4 Amplitude (µV) Error Correct Difference Time (ms) −200 0 200 400 600 800 1000 0 0 . 05 0 .1 0 .1 5 OT2.RIGHT −6 −4 −2 0 2 4 Amplitude (µV) Error Correct Difference Time (ms) −200 0 200 400 600 800 1000 0 0 . 05 0 .1 0 .1 5 −3 −2 −1 0 1 2 3 Best Match: Brodmann Area 24 −5 0 5 RAW EEG SIGNAL FILTERED EEG SIGNAL −6 −4 −2 0 2 4 Amplitude (µV) Error Correct Difference Time (ms) −200 0 200 400 600 800 1000 0 0 . 05 0 .1 0 .1 5 −3 −2 −1 0 1 2 3 Best Match: Brodmann Area 24 NON-ERROR POTENTIALS ERROR POTENTIALS −5 0 5 NON-ERROR POTENTIALS ERROR POTENTIALS −5 0 5 NON-ERROR POTENTIALS ERROR POTENTIALS −5 0 5 Figure 3.3: Time-locked grand averaged signals for the raw EEG data (left) and after data ltering (right) on channel FCz (averaged for all subjects) for OT1 (top), OT2 (down) and averages of each subtask for the raw data (centre). The time range is [−200,1000] ms with 0 being the onset of the action. Error and non-error potentials are in red and blue respectively, and the difference averages (error minus non-error averages) are in dashed lines. The r2 discriminability test [1] between error and non-error potentials is below each plot, where dark colours indicate high values (i.e., large dierences) between the potentials in both conditions. The spatial location of each peak of the dierence average is displayed as topographical scalp maps, as well as the source location of the error grand average at the most prominent negativity ( 500 ms and 450 ms for OT1 and OT2). The single-trial potentials sorted by the negative peak latencies are shown below the source localization as a colour encoded image (red and blue indicate amplitudes higher and lower than 0 µV respectively). 3.3. Results 39 −0.1 −0.05 0 0.05 0.1 Time (ms) Channels 0 100 200 300 400 500 600 700 800 900 1000 32.Oz 31.Pz 30.Cp2 29.Cpz 28.Cp1 27.Cz 26.Fc2 25.Fcz 24.Fc1 23.Fz 22.O2 21.O1 20.P8 19.P4 18.P3 17.P7 16.Cp6 15.Cp5 14.T8 13.C4 12.C3 11.T7 10.Fc6 9.Fc5 8.F8 7.F4 6.F3 5.F7 4.Af4 3.Af3 2.Fp2 1.Fp1 r2 . . . f1 f2 f3 f10 FEATURES OT1 . . .. . .. . . FEATURES OT1.UP FEATURES OT1.LEFT FEATURES OT1.RIGHT OPERATIONAL TASK 1 . . . f1 f2 f3 f10 . . .. . .. . . FEATURES OT2 FEATURES OT2.UP FEATURES OT2.DOWN FEATURES OT2.LEFT OPERATIONAL TASK 2 . . . FEATURES OT2.RIGHT −0.1 −0.05 0 0.05 0.1 Time (ms) Channels 0 100 200 300 400 500 600 700 800 900 1000 32.Oz 31.Pz 30.Cpz 29.Cz 28.Fcz 27.Fz 26.Cp2 25.Cp1 24.Cp6 23.Cp5 22.Fc2 21.Fc1 20.Fc6 19.Fc5 18.Af4 17.Af3 16.O2 15.O1 14.P4 13.P3 12.P8 11.P7 10.C4 9.C3 8.T8 7.T7 6.F4 5.F3 4.F8 3.F7 2.Fp2 1.Fp1 Up vs Left Up vs Right Left vs Right 0 0.5 1 1.5 2 2.5 Divergence (bits) INTRA-TASK OT1 Up vs Left Up vs Right Left vs Right Up vs Left Up vs Right Left vs Right 0 0.5 1 1.5 2 2.5 Divergence (bits) INTRA-TASK OT2 OT1 vs OT2 0 0.5 1 1.5 2 2.5 Divergence (bits) INTER-TASK KL DIVERGENCES Non−error Error r2 Figure 3.4: (Top-Middle) Representative examples of the feature extraction process for each operational task. The r2 metric (Left) was used to choose the time window of [200,800] ms in fronto-central channels and then extract the initial features, as in [79]. These features were the inputs to the spatio-temporal lter, whose outputs were the k -most ( k= 10 for the features analysis) discriminant features, each of them encoding combinations of time points and channels. The weights of some features for each task and subtask are shown as a colour encoded image (blue and red indicate negative and positive weights, respectively). (Lower part of gure) Bar plots of the KL divergences (mean ± SEM) between the features distributions for the inter-task and intra-task conditions (blue and red for the non-error and error distributions). 40 Chapter 3. ErrP task-dependent variations Ten Fold OT2 Ten Fold OT1 Ten Fold OT2.i Ten Fold OT1.i 30 40 50 60 70 80 90 100 Accuracy (%) Non−error Error Train OT2 − Test OT1 Train OT1 − Test OT2 Train OT2.i − Test OT2.j Train OT1.i − Test OT1.j 30 40 50 60 70 80 90 100 Accuracy (%) Non−error Error Figure 3.5: Mean ± std classication accuracies averaged for all the subjects, for the (Left) inter-task and their baseline and (Right) intra-task and their baseline. Blue and red bars indicate accuracies for non-error and error potentials respectively. t-test, p < 1·10−4 ). The inter and intra-task KL divergences of the error distributions for OT1 were signicantly dierent (unpaired two-tailed t-test, t38 = 5.36, p = 4 ·10−6 ), but the inter/intra-task divergences of the non-error distributions were not ( t38 = 0.04, p = 0.97 ). For OT2, inter/intratask KL divergences were signicantly dierent for the non-error distributions ( t68 = 2.69, p = 0.01 ), but no signicant dierences were found for the inter/intra-task divergences of the error distributions ( t68 =−1.05, p = 0.30 ). In summary, feature distributions changed signicantly between tasks and among subtasks, and the error distributions changed signicantly more than the non-error distributions. Furthermore, the features varied signicantly more when changing the task, than when changing the subtask. 3.3.3 Classication Analysis without adaptation The ten-fold accuracies of non-error and error potentials were, on average, 89.29% and 78.00% for OT1 and 86.64% and 73.00% for OT2; and 89.38% and 77.97% for subtasks of OT1, and 84.06% and 72.33% for subtasks of OT2 (see Figure3.5). All baseline classiers were above the chance level. Regarding the inter-task generalization results (Figure 3.5,Left), when training with OT1 and testing with OT2 there was signicant average increase of 6.76% (one-tailed paired t-test, t18 = 1.86, p = 0.04 ) in the detection of non-error potentials, but a signicant decrease of 21.54% ( t18 = −4.19, p = 0.0003 ) for error potentials. As can be seen, the standard deviation was also increased compared to the ten-fold accuracies. This indicated that the accuracy drops varied substantially from subject to subject, with subjects having large drops, and others having almost no accuracy decrease. When training with OT2 and testing with OT1, the accuracies pre- 4|Latency correction of eventrelated potentials between dierent experimental protocols 4.1 Introduction The previous chapter has shown that error-related potentials are aected by the task being performed, changing up to the point of preventing a classi- er to generalize or adapt among dierent tasks. This result was expected, as many neurophysiological studies have shown that event-related potentials (to whom error-related potentials belong) are aected by aspects related to the task performed [50]. In this chapter, we further study the reasons behind this drop in the classication performance. The chapter focuses on two widely studied ERPs: the error-related potentials and the P300 [84]. This requires characterizing the ERPs by the acquisition of enough trials to build a reliable model represented by their grand averages [50]. This is due to the poor signal-to-noise ratio of the EEG as well as several sources of variability that may aect the amplitude or the latency of the ERP components. For instance, the early ERP components (appearing within 200 ms from the stimulus presentation) are aected by application-specic factors such as the spatial attention [113] or the stimuli contrast [50]; as well as user-specic factors such as arousal or valence [49]. In turn, late ERP components (occurring later than 200 ms) are aected by application-specic factors such as the probability of occurrence of the expected stimulus [50] or the interstimulus interval [127]; user-specic factors such as the age and the cognitive capabilities [128,129]; or applicationand user-specic variability such as the stimulus evaluation time (i.e., the amount of time required to perceive and categorize a stimulus) [50,51]. Typically, experiments are designed in a well-controlled manner to reduce the ERP variability. In consequence, it is not clear whether the obtained 47 48 Chapter 4. Latency correction of event-related potentials Latency MODEL GENERALIZATION MODEL GENERALIZATION Latency compensation X EXPERIMENT 1 EXPERIMENT 2 EXPERIMENT 1 EXPERIMENT 2 Figure 4.1: (Left) Example of the latency between two grand averaged event-related potentials elicited from dierent experimental protocols. Such dierence prevents from having classiers that generalize among protocols. (Right) By estimating and removing the latency of the two ERPs, the classier would be able to work under dierent experimental protocols. model also reects the same neural phenomena under dierent performed tasks. This is of particular importance for practical BCI applications where decoding algorithms are expected to keep their performance level irrespective of external factors. Moreover, BCIs often exploit the same brain processes in dierent applications where changes are introduced in the used stimuli, the feedback modality or the controlled device (e.g., see [45,47,83,130] for dierent applications based on error-related processing). In the ideal case, these systems should be able to generalize across dierent operating BCIs independent of the controlled device. In practice, however, there is a need for training a model for each new experimental protocol or session, which is a time-consuming operation and a major issue when deploying BCIs out of the lab. To address this issue, previous researches have tried to reduce this calibration time either by using adaptive classiers [71,83], or by initializing the model with data from a pool of subjects [73,74]. Although previous studies have described the eect of variations in the ERP amplitudes [47] and latencies [131] within the same BCI experimental protocol, the eect of these variations among dierent protocols remains unclear. We hypothesize that it could be possible to build or re-adjust models that compensate for these variations by using information from previous experimental protocols, thus enabling generalization of existing BCI decoders to dierent protocols or applications. In fact, the previous chapter has shown that this is the case for a specic ERP, the error-related potentials (see section 3.3.1). The main idea is depicted in Figure 4.1 Left, where two experimen- 4.2. Experimental methods 49 tal protocols elicit the same ERP with similar waveforms and amplitudes but dierent latencies. If we could estimate the latency variations between the two experimental protocols, previous models could be used in the new protocol after compensating for the latency shift (see Figure 4.1 Right). In this chapter, we analyze the eect of ERP amplitude and latency variations among dierent experimental protocols based on the same cognitive process. We also present a method to analyze and compensate for the latency variations in BCI applications. Two widely used signals were analyzed: the P300 evoked potentials [9,50,84] and the observation error-related potentials (ErrP) [42,45,47]. For each kind of ERP, three dierent experimental protocols with dierent levels of diculty were designed. The latencies between protocols were studied from two points of view: the characteristics of the ERPs and the single-trial classication. The results illustrate (i) how the experimental protocols signicantly aect the latency of the recorded potentials but not the amplitudes, and (ii) how the use of latency-corrected data allows for the generalization of BCI decoders, reducing this way the calibration time when facing a new experimental protocol. This work extends our previous work [85] with a more robust technique to compensate the latencies and shows its application to ERPs of dierent nature. 4.2 Experimental methods We focus on two types of ERPs: the P300 evoked potentials and the observation error-related potentials (ErrP). For each of these signals, three types of experimental protocols were designed (i.e., three dierent ways of evoking the P300 and the ErrPs). 4.2.1 Data recording and experimental setup EEG was recorded with a gUSBAmp amplier ( gTec medical engineering , Schiedelberg, Austria) with 16 active electrodes, with the ground and reference placed on the forehead and the left earlobe. Dierent montages were made for the P300 and ErrP protocols (see details below). The EEG was digitized at 256 Hz, power-line notch ltered at 50 Hz, and zero-phase band-pass ltered at [1,10] Hz. Participants were seated on a comfortable chair facing the visual displays of the protocols approximately one meter away. During all experiments participants were asked to restrict eye movements and blinks to specic resting periods. 50 Chapter 4. Latency correction of event-related potentials Figure 4.2: Experiments performed for the (Top) P300 potentials and (Bottom) observation error potentials (from left to right: experiments 1 to 3). P300 experimental protocols For these protocols we recorded EEG signals with the BCI2000 framework [110] from 16 active electrodes located at Fp1, Fp2, Fz, FC1, FCz, FC2, Cz, CP1, CPz, CP2, P3, Pz, P4, O1, Oz and O2 according to the 10-10 system. Five participants (one female, mean age 27.80 ±2.49 years) took part in the study. We synchronized the onset of visual stimuli with the EEG by means of an optical trigger placed on the monitor [132]. This removed latencies introduced by the protocol implementation and thus the latency variations across experiments were restricted to the user side [51]. Three experimental protocols were used to evoke the P300 potentials (Figure 4.2, Top), with dierent types of stimuli. The stimulation process followed the oddball paradigm [84], where subsets of potential targets (e.g. an entire row or column) are sequentially highlighted in random order. The stimulus (row or column) remained highlighted for 125 ms on the screen, and the inter-stimulus interval was random within the range [1.7,3.0] s. The participants were instructed to observe the stimulation process xing their attention to a given target, and to count the times the target was highlighted while ignoring the other targets. All participants executed the experiments in the same order, each experiment lasting ≈1.5 hours and with a time between experiments of 1.10 ±0.81 days. Experiment 1, 2D Simulated Wheelchair (Figure 4.2 Left, Top) [9] The visual display showed a virtual environment with 20 possible targets 4.2. Experimental methods 51 to drive a wheelchair, located in 2D in a 4x5 matrix. For the stimulation process, the rows and columns were highlighted showing a blue dot over each possible target position. The probability of target appearance was 22% . For each subject, all possible target positions were recorded, obtaining 144 target (P300) and 720 non-target responses respectively. Experiment 2, 2D Speller (Figure 4.2 Middle, Top) [84] The visual display showed a matrix of 36 possible letters to spell represented in 2D as a 6x6 matrix. The stimulation was made by illuminating the corresponding row or column. The probability of target appearance was 17% . For each subject, all possible target positions were recorded, obtaining 200 and 700 target and non-target responses respectively. Experiment 3, 3D Augmented Reality Protocol (Figure 4.2 Right, Top) The display showed a gray background and 27 possible targets located in 3D in a 3x3x3 matrix. The stimulation was made by illuminating rows, columns, and depths. To facilitate the user's distinction among the three depths, each depth was illuminated with a dierent colour (green, blue or red). The probability of target appearance was 33% . For each subject, all possible target positions were recorded, obtaining 273 and 610 target and non-target responses respectively. Error potentials experimental protocols We recorded ErrPs with a custom C++ framework using 16 active electrodes located at Fz, FC3, FC1, FCz, FC2, FC4, C3, C1, Cz, C2, C4, CP3, CP1, CPz, CP2, and CP4 according to the 10-10 system. Six participants (one female, mean age 27.33 ±2.73 years) took part in the study. In these experiments, the use of an optical trigger was not possible since one experiment involved a real robotic device instead of visual stimuli on the screen (see Experiment 3). Thus, latency variations could be originated by both the subject and the implementation (i.e. the amount of time of receiving and executing the delivered command). The experimental protocols designed to elicit error potentials (Figure 4.2, Bottom) had dierent setups (and devices), where in all cases the goal of the device was to reach a target from dierent starting points. The device executed random movements with approximately 30% probability of performing an erroneous movement. The time between two movements was random within the range [1.7,4.0] s. The target position was randomly changed after 100 actions. The participants were instructed to observe the device movements and evaluate them as correct when there was progress towards the 52 Chapter 4. Latency correction of event-related potentials target position, and as incorrect otherwise. Each participant executed the experiments in the same order, each experiment lasting ≈2.5 hours and with a time between experiments of 17.58 ±10.09 days. Experiment 1, Virtual Moving Square (Figure 4.2 Left, Bottom) [47] The visual display showed a one-dimensional space with 9 possible positions (marked by a horizontal grid), a blue square (device) and a red square (target). The device could execute two discrete actions: move one position to the left or to the right. For each subject, the leftand rightmost positions were tested as targets, and around 250 and 600 error and non-error potentials were recorded. Experiment 2, Simulated Robotic Arm (Figure 4.2 Middle, Bottom) The display showed a simulated robotic arm (Barrett WAM) with 7 degrees of freedom (device) [133] moving within a two-dimensional space with 13 possible positions (marked in orange), and a target location (green square). The robot was situated behind the squares pointing at one position, and could perform four possible actions: moving one position to the left, right, up, or down. The robot actions were continuous, with each displacement lasting ≈500 ms. For each subject, the left-, right-, upand down-most positions were tested as targets, and around 300 and 700 error and non-error potentials were recorded. Experiment 3, Real Robotic Arm (Figure 4.2 Right, Bottom) This experiment followed the same conguration of Experiment 2 but using a real Barret WAM robotic arm ( Barret Technology Inc. ). The user was seated two meters away from the robot, and between them there was a transparent panel to mark the positions (the distance between two neighbor positions was 15 cm). For each subject, the left-, right-, upand down-most positions were tested as targets, and around 300 and 700 error and non-error potentials were recorded. 4.2.2 Analysis of Event-Related Potentials We assessed protocol-dependent variations in the latency and amplitude of the ERPs for each of the experimental protocols. First, the grand averaged signals were computed for each condition (target and non-target trials for the P300; error and correct trials for the ErrP), for the time window [−200,1000] ms, being 0 ms the stimulus/action onset. Following previous studies, we analyzed the activity over parietal areas from the target average [50] for the 4.2. Experimental methods 53 P300, and over fronto-central areas from the dierence average (error minus correct averages) for the ErrPs [47]. A one-way within-subjects ANOVA was performed separately for each type of signal (P300 or ErrP), where the factor was the experiment (three levels corresponding to each experiment), and two dependent variables were tested: the peak amplitudes and the peak latencies. For the P300 experiments, the peak amplitudes and latencies were measured from the P3 component (most prominent positive peak) of the target average from the parieto-occipital channels. For the ErrP experiments, the amplitudes and latencies were measured from the P3 and N4 components (most prominent positive and negative peaks) of the dierence average from the fronto-central channels. When needed, the Geisser-Greenhouse correction was applied to data to assure sphericity [50]. Pairwise post-hoc tests (t-tests with the Bonferroni correction) were computed to determine the dierences between pairs of experiments. 4.2.3 Estimation and evaluation of latencies among different protocols The rst goal is to estimate the temporal variations between two experimental protocols, which can be achieved using the cross-correlation. The crosscorrelation has been used in the past for the detection and analysis of brain signals with successful results [131,134,135]. In order to assure the best estimations, the input to the cross correlation (for each channel) were the grand averages of the condition of interest, with the time window narrowing to the event-related potential elicitation. For the P300 experimental protocols, the average ERP for target stimuli within the time window [50,400] ms was used; for the ErrP experimental protocols, the error average within the time window [0,500] ms was used. The outputs of the cross-correlation were the maximum correlation value of the two grand averages and the latency variation between them (i.e. the shift that yields the maximum cross-correlation). We then assessed whether the main ERP change was due to the latency variation and whether this variation could be compensated for. To do so, the latency variation across two protocols was estimated as described above using all the available data. Let Ei and Ej be the datasets from two experiments i and j , we compensate for the variation by shifting the trials in Ej by the estimated latency shift between them, dEiEj . Then, we computed the same ANOVA test for the peak latencies performed in subsection 4.2.2. We performed further analysis on how sensitive the latency estimation was with respect to (i) the number of trials used to compute the grand average for experiment j , and (ii) the channel used to perform the estimation. Assuming 54 Chapter 4. Latency correction of event-related potentials that data from a previous experiment i is available, we computed the latency variation using training data, ET r j , from the new experiment ( ET r j⊆Ej ). We assessed the estimation using dierent sizes of the training dataset (ranging from 10 to 200 trials with increments of 10). For each size, we perform 10 repetitions and report the average of the maximum cross-correlation value, max(CET r j Ei) , and the average latency variation, dEiET r j . In each repetition the training subset ET r j was randomly drawn from Ej , keeping the proportion of target/non-target and error/correct trial. The analysis was performed independently for each recorded channel. The latency variations were computed in a pair-wise manner among the three experiments for each of the signals of interest. The combinations of experiments tested were E1E2 , E1E3 , and E2E3 for both the P300 and the ErrPs. For each pair of experiments a within-subjects two-way ANOVA (factors: number of trials and repetitions) was performed on the latency estimations. The ANOVA results served to study the latency variations by the number-of-trials main eect, to determine whether the amount of trials used from Ej led to dierent latency estimations; and by the number of trials x repetitions interaction, to determine whether dierent data from a xed number of trials aected the latency estimations. As a sanity check, we also evaluated the method by computing the latency variation among datasets from the same experiment ( dE1 iE2 i ), with the two datasets E1 i and E2 i mutually exclusive. Therefore, this baseline latency computation should give correlations close to one for latencies near to 0 ms. 4.2.4 Single-trial classication of latency-corrected ERPs The objective of the single-trial classication study was to determine whether it is possible to reduce the calibration time of a new experiment by re-using latency-corrected data from a previous experiment. The same combination of experiments detailed before served as the basis of the study. The process of the latency estimations used the number of trials and repetitions tested in the previous subsection. All the classication parameters (including the latency estimations) were learned using only the training sets. Feature extraction and classication Feature extraction was based on a spatio-temporal lter [80]. The lter input was a dataset with labeled trials, and worked as follows: Firstly, the EEG data were common-average-reference (CAR) ltered and downsampled to 64 Hz. For each trial, the features were extracted using a combination of 4.2. Experimental methods 55 channels and time points. For the P300, eight centro-parietal and occipital channels (Cz, CPz, P3, Pz, P4, O1, Oz, and O2) were used within a time window of [100,700] ms. For the ErrP, eight fronto-central channels (Fz, FC1, FCz, FC2, C1, Cz, C2, and CPz) were used within a time window of [200,800] ms. For both cases, this resulted in a feature vector of 312 features per trial. Then, the features were normalized, and decorrelated using PCA retaining 95% of the explained variance, leading to an average of 45 ±10 features. Single-trial classication was carried out using a linear discriminant analysis (LDA) [124]. Analysis of the single-trial classication We compared the accuracies of three dierent classiers (see Figure 4.3). The rst one, denoted baseline classier, followed the standard calibration approach of current BCIs, where the classier for Ej was trained using only the data from ET r j . The size of the training data was increased as in the previous subsection to assess the accuracy of the classier for dierent calibration times. The second and third classiers were both trained using all the data from Ei and the data from the new experiment ET r j . The dierence between them is that in the latter one the latency of dataset Ei was corrected for each channel used for classication, c.f. subsection 4.2.4. These three classiers were tested on a left-apart xed subset of Ej (denoted ET e j ), composed of 400 trials (see Figure 4.3). Additionally, the performance of the three classiers was also compared to the result of performing ten-fold cross-validation (CV) using all the data from Ej . For each pair of experiments ( E1E2 , E1E3 , and E2E3 ), two metrics are reported: (i) the mean accuracy ( target _ acc+non _ target _ acc 2 and error _ acc+correct _ acc 2 for the P300 and ErrP experiments respectively); and (ii) the accuracy bias ( |target _ acc−non _ target _ acc| and |error _ acc−correct _ acc| for the P300 and ErrP experiments). To assess statistical dierences among the classi- cation results, two-tailed paired t-tests were computed, where the p-values were adjusted with the false discovery rate (FDR) procedure [136]. 4.2.5 Reducing calibration time during online sessions of experiments Finally, the latency estimation method was used for real online experiments. For this analysis, 6 new subjects (two females, mean age 27.17 ±4.07 years) performed the error-related potentials experiments where the calibration time was reduced using latency-corrected data from the previous experiment. For E1 , the standard calibration was followed. On the other hand, during E2 the 56 Chapter 4. Latency correction of event-related potentials Not used Training set Testing set Baseline Classifier EiEj Classifiers when not correcting or correcting the latency Variable size Tr Ej Te 400 trials Figure 4.3: Training and testing datasets used for each classier. For the baseline classier, the training dataset from the new experiment ET r j was used to train the classier. When information from a previous experiment Ei was available, it was used as part of the training set, either not correcting or correcting the latency dEiEj . To correct the latency, the training datasets Ei and ET r j were fed to the latency correction algorithm, and then the latency from Ei was corrected. To assess the accuracies for dierent calibration times, the number of trials of ET r j was variable (from 10 to 200 trials with increments of 10 trials), while the testing set ET e j remained xed to 400 trials. data from E1 was latency-corrected using a few trials from E2 . Similarly, during E3 the data from E2 was used. The latency between experiments was estimated using the dierence average of channel FCz within the window [0,500] ms. The groups of subjects were denoted control group (the subjects present in this work, where they followed a standard calibration approach) and experimental group (latency-corrected calibration). The calibration phase nished whenever a ten-fold mean accuracy of 75% was obtained in the training data. Then, both groups of subjects performed the online experiment presented in chapter 2, which was use to extract the mean online accuracy obtained. Mixed two-way ANOVAs (within factor: experiments; between factor: group of subjects) were performed to test whether (i) the number of calibration trials in the experimental group decreased from E1 to E2 and from E2 to E3 ; (ii) the number of calibration trials was signicantly dierent between groups; and (iii) the mean online accuracies of both groups were not dierent. In order to nd the underlying signicances, post-hoc one-tailed Bonferroni-corrected t-tests were performed. 4.3 Results 4.3.1 Analysis of Event-Related Potentials Figure 4.4 shows the ERP grand averages of all experiments. In the P300 experiments, as in previous studies [9,50,84], a clear sharp positive peak (P3) 4.3. Results 63 E E 0 50 100 150 200 0 50 100 Trials used from new experiment Bias (%) 0 50 100 150 200 50 60 70 80 Trials used from new experiment Accuracy (%) 0 50 100 150 200 50 60 70 80 Trials used from new experiment Accuracy (%) 0 50 100 150 200 0 50 100 Trials used from new experiment Bias (%) 0 50 100 150 200 50 60 70 80 Trials used from new experiment Accuracy (%) 0 50 100 150 200 0 50 100 Trials used from new experiment Bias (%) 1 2 E E 1 3 E E 2 3 ErrP EXPERIMENTS Baseline Not Correcting Latency Correcting Latency Figure 4.8: (Top) Mean accuracy and (Bottom) classier bias when correcting the latency from E1E2 , E1E3 and E2E3 for the ErrP experiments. allowed for an improvement both in the classier accuracy and calibration time. However, at least 90 trials were required for the latency correction method to perform as well as the no correction approach, seemingly due to errors in the latency estimation. The baseline classier showed large bias when few trials are available. This bias decreased as more trials were added reaching 23.85% after 200 trials. In contrast, re-using data from the previous experiment yielded a low bias since the beginning, both when correcting and not correcting the latency ( 15.70% and 13.70% after 200 trials, respectively). These values were similar to the ten-fold CV. As with the accuracy, the baseline classier bias was always signicantly worse than the results obtained with the other two classiers ( p < 0.0001 ). Compared to the previous case, going from E1 to E3 resulted in lower accuracies for all types of classiers (c.f. Figure 4.7 central column), always lower than the CV accuracies. After 200 trials the accuracies were of 59.33% , 58.89% and 56.48% for the baseline, latency-corrected and latency non-corrected classiers, respectively. When the latency was not corrected, signicantly lower accuracies were obtained with more than 100 trials ( p < 0.05 ). No signicant improvement was found when correcting the latency with respect to the baseline ( p > 0.05 ). This suggests a smaller eect of re-using data from the previous experiment. Furthermore, the baseline classier exhibited signicantly lower bias ( p < 0.0001 ) than the other classiers, reaching 12.67% after 200 trials, versus 26.60% and 28.59% when correcting and not correcting the latency, respectively. In the last case ( E2E3 , c.f. Figure 4.7 right), the latency correction mechanism yielded signicantly higher accuracies ( p < 0.05 ) than the baseline or no correction approaches ( 63.09% , 59.85% , and 60.31% after 200 trials, re- 64 Chapter 4. Latency correction of event-related potentials spectively), converging to the accuracy of the 10-fold CV. These dierences appeared even when a small number of trials were available. Thus, the use of latency-corrected data allowed for a signicant improvement in the accuracies. Moreover, although after 200 trials the bias level was similar for all approaches (around 12% ), not correcting the latency yielded signicantly higher bias than the other two classiers in almost all the cases ( p < 0.05 ). No signicant dierences were found between the latency-corrected and the baseline bias. Error potentials Results for the error potential protocols are shown in Figure 4.8. In the rst case, E1E2 , the accuracy of the baseline classier reached 72.05% after using 200 trials for training. The latency-corrected classier showed better performance of 74.27% , similar to the ten-fold CV accuracy. Notably, the latter always performed signicantly better than the baseline (two-tailed paired ttest, p < 0.0001 ). In contrast, the use of previous data without correcting the latency led to signicantly lower accuracies than the other classiers for less than 120 trials ( p < 0.05 ), reaching 62.81% after 200 trials. Furthermore, despite its high accuracy, the baseline classier had signicantly higher bias than the other classiers ( p < 0.0001 ), whereas the latency correction yielded to a signicantly lower bias than the classier when not correcting the latency ( p < 0.0001 ). In the second case, E1E3 , the latency-corrected classier signicantly outperformed the baseline when a small number of trials were used ( p < 0.05 , with less than 100 trials). These classiers performed similarly after 200 trials (with accuracies of around 73% ), without reaching the accuracy of the ten-fold CV. Again, the classier using not corrected data always performed signicantly worse than the others ( p < 0.0001 ), with an accuracy of 62.37% after 200 trials. Similarly, the bias of the latency-corrected classier was signicantly lower than the bias of the other two classiers when less than 160 trials are used ( p < 0.05 ). In the last case ( E2E3 , c.f. Figure 4.8 right), the baseline classier was always signicantly worse than the latency-corrected one ( p < 0.0001 ). After 200 trials, the baseline classier reached a 73.10% of mean accuracy versus a 76.60% when correcting the latency, comparable to the ten-fold CV classi- er. The latter classier was also signicantly better than the non-corrected classier with 20 trials or more ( p < 0.0001 ). Finally, the latency-corrected classier obtained signicantly lower biases than those obtained with the baseline classier for less than 130 trials ( p < 0.05 ). To summarize, apart from one case generalization from the P300 ex- 4.4. Discussion 65 periments E1 to E3  the latency-correction mechanism improved the classi- cation performance in all the analyzed cases. It allows signicantly higher accuracies and/or lower bias than the baseline condition when a small number of trials from the new experiment are available. Exp. 1 Exp. 2 Exp. 3 0 100 200 300 400 500 600 # Calibration trials Control Experimental (a) Exp. 1 Exp. 2 Exp. 3 0 20 40 60 80 100 Mean accuracy (%) Control Experimental (b) Figure 4.9: (a) Mean ± SEM calibration trials needed for the control group (red) and experimental group (blue), for each experiment. Signicant dierences ( p < 0.05 ) between pairs are marked with a star. (b) Mean classier accuracy during the online control phase for each group and experiment. 4.3.4 Online results Figure 4.9a shows the number of calibration trials needed for each group of subjects and experiment. The ANOVA test revealed a signicant interaction between the experiment and group ( F2,20 = 8.65, p = 0.002 ). Post-hoc tests showed that signicant dierences were found between groups in experiments 2 and 3 (one-tailed unpaired t-tests, p < 0.05 ), and also signi- cant dierences within the experimental group between experiments 1 and 2 (one-tailed paired t-test, p= 0.004 ), and between experiments 1 and 3 ( p= 0.004 ), see Figure 4.9a. On the other hand, the accuracies was significantly similar for all the experiments and subjects ( p > 0.85 ) (see Figure 4.9b). These results indicate that, given data from previous experiments, the latency correction algorithm allowed for a reduction in the calibration time of error-related potentials. 4.4 Discussion A practical issue in the study of event-related brain activity and its use for BCI applications is the time required to acquire sucient data to have a 66 Chapter 4. Latency correction of event-related potentials reliable model or a usable classier of the EEG signals. In general, each protocol is treated as a completely new experiment even if they are tapping into similar cognitive processes. Besides the increase in the required time and resources, this also provides little information about how similar responses are across experimental conditions. Using several protocols on two wellstudied signals, we showed that the experimental design mainly aected the ERP latencies. Moreover, we proposed a simple, yet powerful mechanism to compensate for these changes allowing to generalize BCI classiers across experiments using a reduced amount of new data. Although previous studies have reported modulations on the P300 and ErrP amplitude depending on the target or error probability [47,50], variations in our protocols did not result in statistically signicant amplitude dierences across experiments. In contrast, ERP latencies were found to be dierent across several experiments (c.f. Section 4.3.1). In the P300 experiments, signicant variations appeared when changing from experiments 1 or 2 to the third experiment. Interestingly, this experiment had the most complex visual stimuli (a three-dimensional grid), seemingly requiring the subject longer time to evaluate the stimulus. This is supported by neurophysiological studies suggesting that the P300 is related to the stimulus evaluation time [51,137]. Regarding the error potentials experiments, the latency changes were larger for both peaks (P3 and N4) than for P300 when changing the experimental conditions. The selected protocols were designed so as to have an increased level of complexity both in the number and type of possible actions: changing from two to four possible actions at each state (from E1 to E2 and E3 ); changing from 1D to 2D (from E1 to E2 and E3 ); and changing from a simulated to a real device with a wider eld of view needed (from E2 to E3 ). Accordingly, increasingly longer latencies were found from protocols E1 to E3 . However, it should be noticed that part of the measured latencies may be also due to implementation issues (e.g. the time it takes to the robot to actually start moving once the control command has been issued). Nonetheless, the use of a simple technique such as the cross-correlation allowed to signicantly remove these latency jitters among experiments, as presented by the ANOVA results after correcting the latencies. As expected, the latency estimation improves when more data is available. Nonetheless, 100 trials proved to be sucient to obtain a reliable estimation with correlations higher than 0.8 as seen on Figures 4.5 and 4.6. In this work, the latency was estimated and corrected separately for each channel. However, in principle more reliable estimations and higher correlation values should be obtained in those areas where the ERP is generated. This information could be combined to study protocol-dependent spatio-temporal 4.4. Discussion 67 variations of a given brain process. Focusing on applications of brain-computer interfacing, we propose a simple latency correction mechanism to re-use data from previous experiments when building classiers for new experiments on a related phenomenon. This yields a reduction in the calibration time as a smaller number of trials is required to achieve similar classication performance (in terms of accuracy and bias) than if a new classier is built from scratch (c.f. Figures 4.7 and 4.8). Moreover, in those cases where the experimental change has little eect on the ERP latencies, the eect of the latency correction mechanism does not have a negative eect on the classication performance (e.g. moving from P300 experiment E1 to E2 ). In a similar way, Thompson et al. [131] also found that the latency variations among trials but within the same experimental protocol were one of the main problems for the classication performance, and proposed the use of within-experiment latency variations as a predictor of online BCI accuracies. The authors also argued that a brute-force method (i.e. testing a classier for each possible latency and taking the classier with maximum accuracy) could be used to estimate these latencies. The latency variations have been assessed on two dierent ERPs (P300 and observation error potentials), showing their eect on the single-trial classication. In the future, the generalization of BCI decoders across protocols can be assessed for other ERPs, such as those generated during rapid visual processing [138], or the N2 evoked component [139]. Moreover, additional studies of event-related potentials in controlled and non-controlled applications may yield new ndings. The proposed method could be used to elucidate common patterns across conditions, not only in BCI applications but also in neurophysiological studies, e.g. comparing latency variations between error-related activity in choice reaction tasks [42], and in feedback tasks [32]. Furthermore, more sophisticated techniques could be tested to cope with the latency variations such as dynamic time warping [140,141]. Finally, one disadvantage of the proposed approach is that it relies on the assumption that there are only temporal changes in the ERPs, whereas the spatial contributions remain xed among experiments. However, this assumption may be wrong. Thus, a more complete approach could be designed by performing a spatio-temporal compensation of the ERP variations. 68 Chapter 4. Latency correction of event-related potentials 5|Shared-control BCI for a two dimensional reaching task using error-related potentials 5.1 Introduction Chapter 2 presented an alternative BMI paradigm where the user was able to teach a device how to solve a task. Besides the calibration phase needed prior to the control phase (addressed in chapters 3 and 4, and further studied in chapter 6), the bottleneck of the proposed paradigm is the scalability of the system. Due to the reinforcement learning framework, the device needed to explore the entire state space to learn the optimal behavior for each position. Furthermore, although ErrPs provide feedback about the device actions, the amount of information conveyed by them is limited. In particular, as shown in the previous chapters, the decoders do not contain any information about direction or magnitude, and have a non-negligible number of misdetections. As a result, the system presented in chapter 2 required more time to converge as the task space increased its state-action complexity (e.g. the system proposed needed more time to converge in experiments 2 and 3 than in experiment 1, see section 2.3). One way of solving these issues is to further increase the controller intelligence. In this sense, recent approaches started to explore shared-control strategies where the device does not only execute the decoded commands, but is also involved in executing the task [2] (e.g. taking into account the environment while reaching a target or avoiding an obstacle [9,10]). In this chapter, we propose the use of a shared-control scheme to greatly reduce the time needed to solve a task in the BMI paradigm proposed in chapter 2, by exploiting the intrinsic constraints of the task being performed. In this work, we present a 2D reaching task of a cursor over a discrete grid of possible targets (similar to the experiments performed in chapter 2 but 69 70 Chapter 5. Shared-control BCI using error-related potentials with greater complexity), using error-related potentials as supervision signals. Under this control scheme, the user has to evaluate whether the cursor is correctly moving towards the goal or is making wrong movements. Based on this signal, the device estimates which is the desired target while moving towards it based on the optimal motion policies for each of the possible goals. The results show that all users were able to reach predened and self-selected goal locations in around 23 actions (about 19 seconds of EEG signal). The use of this scheme has several interesting advantages with respect to the paradigm presented previously. First, it is possible to learn the entire motor behavior space using much less actions, and at the same time reach the desired target; secondly, it oers a natural way of not only learning motor behaviors, but also choosing among a set of previously stored ones; third, within a shared-control strategy it can be combined with other brain signals to correct decoding failures or recover from wrong or ambiguous decisions of the device. 5.2 Methods 5.2.1 Data Recording Electroencephalographic (EEG) and electrooculographic (EOG) activity were recorded using a gTec system (3 synchronized gUSBamp ampliers). For the EEG, 32 electrodes were recorded, distributed according to an extended 10/20 international system (FP1, FP2, F7, F8, F3, F4, T7, T8, C3, C4, P7, P8, P3, P4, O1, O2, AF3, AF4, FC5, FC6, FC1, FC2, CP5, CP6, CP1, CP2, Fz, FCz, Cz, CPz, Pz and Oz), with the ground on FPz and the reference on the left earlobe; for the EOG, 6 monopolar electrodes were recorded (placed above and below each eye, and from the outer canthi of the left and right eyes [142]), with the ground on FPz and the reference on the left mastoid. The EEG and EOG signals were digitized with a sampling frequency of 256 Hz, power-line notch ltered, and band-pass ltered at [0.5,10] Hz. The EEG was also common-average-reference (CAR) ltered. Additionally, the horizontal, vertical, and radial EOG were computed as in [142] to remove the EOG from the EEG using a regression algorithm [143]. The data acquisition and online processing was developed under a self-made BCI platform. 5.2.2 Experimental Protocol Four subjects (mean age 26 ±2 years) participated in the experiments. The participants were seated on a comfortable chair one meter away of a computer 5.2. Methods 71 (a) (b) (c) Figure 5.1: (a) Experimental protocol designed. The protocol showed a 5x5 grid with a virtual cursor (green circle) and a goal location (shadowed in red). The up-left-most position is the [1 1] (row and column). (b) The cursor could perform ve dierent actions (from top to bottom, move one position up, down, left or right, or performing a goal-reached action). (c) Correct actions (i.e. optimal policy) from each state for the goal exemplied on (a). screen that displayed all the information related to the experiment. The experimental protocol is shown in Fig. 5.1. The protocol consisted of a 5x5 squares grid, a virtual cursor (green circle), and a goal location shadowed in red. The cursor performed ve dierent instantaneous actions: move one position left, right, up or down; and a goal-reached action, represented as concentric blue circumferences (see Fig. 5.1b). The time between two actions (inter-action interval) was random within the range [3,3.5] s. The role of the subjects was to assess the cursor actions as correct when the cursor performed (i) a movement towards the goal position, or (ii) a goal-reached action over the goal position; and as incorrect otherwise (see Fig. 5.1c). The users were instructed not to move their eyes during the cursor actions, and to restrain blinks only to the resting periods. The experiment consisted of two phases: the training phase used to calibrate a classier able to distinguish between correct and error user's assessments; and the control phase, where the user controlled the cursor to a goal position. During the control phase, two dierent groups of goal locations were tested: (i) the rst group was xed for all the subjects, and consisted of ve predened goals and initial cursor positions (see table 5.1); (ii) for the second group, each user was asked to freely choose ve dierent initial cursor 72 Chapter 5. Shared-control BCI using error-related potentials 0 0.2 0.4 0.6 0.8 1 G (a) G (c) G (b) Figure 5.2: Likelihoods of each policy πi after performing dierent actions: (a) correct movement with p(ct= 1|xt)=0.2 (b) incorrect movement p(ct= 1|xt) = 0.8 (c) or a goal-reached action p(ct= 1|xt)=0.2 . The goal position is marked with a capital G. positions and goals to reach. During this group of goals, the goal position was not shadowed in red, since it was the user who chose it. Table 5.1: Initial and goal positions for the xed group of goals Run 1 Run 2 Run 3 Run 4 Run 5 Initial position [1 1] [2 3] [1 2] [1 4] [3 3] Goal position [5 5] [3 1] [3 3] [4 1] [3 3] 5.2.3 Calibration of error potentials For the calibration of the error potential detection, a training phase was rst executed to acquire sucient examples of potentials while the user assessed the device actions. During this phase, the virtual cursor performed random actions, with a 20% of probability of performing an erroneous action. This phase lasted for 30 minutes, acquiring around 80 correct and 320 erroneous trials. Once the training data was acquired, features were extracted from eight fronto-central channels (Fz, FC1, FCz, FC2, C1, Cz, C2, and CPz) within a time window of [200,800] ms downsampled to 64 Hz, forming a vector of 312 features. The feature vectors of all trials were normalized and decorrelated using PCA, retaining those that explained 95% of the variance. Finally, a regularized linear discriminant (LDA) [124] classier was trained using the retained features. The classier output has the form y( x ) = w 0 x +b , where y( x )≥0 is classied as a correct assessment (class 1), and y( x )<0 as 6|Zero-calibration BMIs for sequential tasks using error-related potentials 6.1 Introduction Previous chapters have shown how it is possible to teach a device how to solve a task following the proposed teaching paradigm. As for all non-invasive BMIs, this paradigm needed a decoder trained after a calibration phase that translated the error-related potentials into feedback for a device. Despite we have seen how this calibration phase can be reduced thanks to re-using data from previous experiments, several examples from the new task were always needed to train the decoder. Furthermore, these examples needed to be acquired in open loop, that is, where the user is not performing the task. This is a key problem for the paradigm developed throughout this thesis, since the time needed to learn the decoder adds to the time needed to learn the task as such, greatly delaying the nal device operation and hindering the deployment of the BMI out of the lab. This chapter addresses the problem of removing the calibration phase in the BMI paradigm proposed in chapter 2 and extended in chapter 5. Here, we show how it is possible to fuse both the calibration and the task learning phases in a single closed-loop phase that learns the classier and the task at the same time and in an unsupervised manner. The main idea is the same as the one presented in the previous chapter: to exploit the task constraints. The previous chapter showed how these constraints allowed to estimate all the possible optimal policies (one for each target) while a recursive lter estimated which target was the correct one. Similarly, it is possible to compute all the possible signal labels (i.e. whether the current EEG is generated from a wrong or correct assessment) associated to each optimal policy, and use a recursive lter that learns which labels best t the 79 80 Chapter 6. Zero-calibration BMI using error-related potentials EEG signals. Related work on this domain is very limited. This idea was addressed in a similar way by Kindermans et al. for P300 signals by exploiting the fact that these ERPs usually require of multiple visual stimulations [74,77]. Thus, a recursive lter was able to unsupervisedly estimate the optimal labeling for the signals. Furthermore, it was possible to further increase the convergence speed by using prior information from a pool of subjects. In invasive BCIs, Orsborn et al. [145] learned from scratch in closed loop a decoder for known targets using pre-dened policies to each target. However, the approach needed for a warm-up period of around 15 minutes. The main contribution of this chapter is a method that simultaneously calibrates the feedback decoder and executes in closed loop a task only known by the user. Our method exploit the described task constraints, namely optimal policies, to hallucinate virtual labels for the feedback signals. Using these labels, it is possible to compute the expected classication rate of the decoder of the task. Since the right task will assign the right labels to the EEG signals, the expected classication rate is also a good measure to identify the user's intended task. Furthermore, we show that it is possible to use model-based planning based on the uncertainty about the task and the feedback signals to explore the space eciently while learning. The method has been evaluated performing online experiments where four users directly controlled an agent on a 2D grid world to reach a target. The results show that the proposed method is able to learn good feedback models and solve the task eciently without any calibration. Oine experiments show that our unsupervised trained decoder has the same accuracy as a standard one and illustrate the benets of our active strategy. 6.2 BCI-feedback Based Control Without Explicit Calibration 6.2.1 Experimental protocol The experimental protocol followed the same reaching task as the one described in the previous chapter (see subsection 5.2.2), where a device was placed in a 5x5 grid world and whose objective was to reach a target position. The user goal was to assess the device actions (thus generating error-related potentials) and teach it to reach the desired position (see Figure 5.1). 6.2. BCI-feedback Based Control Without Explicit Calibration 81 TT Class 0 Class 1 T Figure 6.1: Task labels for a 1-D grid world. For the represented example, the arrows indicate for each state what action should elicit a positive feedback to reach the target position shadowed in blue (i.e., the optimal policies). 2D Gaussian distributions of binary feedback signals for three possible targets are shown below. While for the correct target the distributions shows a large separability (Left), the overlaps increases as the believed target moves away from the real one (Middle, Right). 6.2.2 Simultaneous Estimation of Task and Signal Model In a common BCI scenario, the EEG signals are usually trained in open loop, where the user has no control whatsoever over the device (see for instance subsection 5.2.3). Whenever the calibration is done, the user performs the closed loop experiment, where it controls (or teaches) the device. In this chapter, we faced the question of whether it is possible to estimate at the same time the task being performed (reach a target position) and the signal model (binary classier trained with ErrP signals). The main idea of this work is depicted in Fig. 6.1 for a toy 1D example with 7 possibles targets to reach (i.e. 7 possible tasks). The user wants the device to reach the right-most state. However, neither the target nor the signal model are known. The feedback signals (ErrPs) are generated as a response to the execution of an action a in state s according to the true unknown task the user wants to solve. The key point is that these signals are generated from an underlying model that for binary signals has two dierent classes. Given sucient feedback signals, it is possible to build the underlying distributions for each possible target. Only the right task will provide the right meanings (or labels) to each of the feedback signals (Fig. 6.1 Left), while the other tasks will gradually mix both classes as the task gradually diers more to the original task (Fig. 6.1 Middle-Right), up to the point of almost mirroring the labels when the target is mirrored. If we extrapolate this to the proposed protocol, it is clear to see that there are 25 possible tasks (one for each possible target). In the remainder of this section we show how this property can be exploited to estimate the task and the model generating 82 Chapter 6. Zero-calibration BMI using error-related potentials the feedback signals. Let ei∈Rn be the EEG measurements e obtained after the execution of action ai in state si . The meaning (or label) zi∈ {c, w} of each feedback signal belongs to one of two classes (correct or wrong). Following the literature [124] and the previous chapters, the EEG signals are modeled using independent multivariate normal distributions for each class, N(µc,Σc),N(µw,Σw) . Let θ be the set of parameters {µc,Σc, µw,Σw} . Regarding the tasks, the system has access to a set of task hypotheses ξ1, . . . , ξT which includes the task the user wants to solve 1 . We do not make any particular assumption on how the task is represented given that for each particular task ξ we are able to compute a policy π which represents the probability of choosing a given action a in state s , πξ(s, a) = p(a|s, ξ) . As mentioned above, these are the policies that, conditioned on the task, provide meanings to the feedback signals of an action-state pair. (e.g. in the proposed reaching task, progressing towards the goal will generate correct answers while moving apart from it will generate wrong ones). Our goal is to learn which task ˆ ξ the user wants to solve based on the assessment of the user extracted from EEG measurements collected while executing actions. Thus collected data are in the form {(ei, si, ai), i = 1, . . . , N} , this is, a sequence of states, actions and teaching signals triplets. Following the discussion of Fig. 6.1, a straightforward option to estimate the task ξ is to measure the coherence of the signal model for each possible task using the virtual meanings provided by the target policy. In other words, the best (ξ, θ) pair should provide the lowest predictive error (perr) on the observed signals p(e|s, a, ξ, θ) . One possible way of solving this problem is to maximize the expected predictive classication rate: ˆ ξ, ˆ θ=argmaxξ,θEe(δ(z(s, a, ξ), z(e, θ))) (6.1) where δ() is an indicator function, z(s, a, ξ) is the label (wrong or correct) corresponding to the execution of action a in state s under task ξ and z(e, θξ) is the label provided by the Gaussian classier with parameters θξ . The expected predictive error can be explicitly written dependent on the task and decoder model: Ee(δ(z(s, a, ξ), z(e, θ))) = X l=c,w p(z=l|s, a, ξ)p(z=l|e, θ) (6.2) where p(z=l|s, a, ξ) represents the probability of the user assigning label l when assessing task ξ . We add a noise term to cope with those situations 1 If this is not the case, the system will nd the most suitable task. 6.2. BCI-feedback Based Control Without Explicit Calibration 83 were the user assessment may be wrong. The model is then p(z=l|s, a, ξ) = (1−α, a = argmaxaπξ(s, a) α, otherwise (6.3) with α modeling error rate of the user. The term p(z=l|e, θξ) is just the probability of the meaning under the Gaussian model provided by θξ p(z=l|e, θ) = p(e|z=l, θ)p(z=l) Pk=c,w p(e|z=k, θ)p(z=k)=N(e|µl,Σl)p(z=l) Pk=c,w N(e|µk,Σk)p(z=k) (6.4) Note that the optimization process has be factored using the fact that given the task, the estimation of θ under the Gaussian model is trivial. It basically requires to compute the maximum-likelihood estimate θML ξ for each target ξ . Using the labels of target ξ , the estimation of θ under the Gaussian model described above simply accounts for computing the posterior mean µz and covariance Σz . In order to avoid numerical problems when estimating the covariance for a low number of examples, a regularization term was applied to penalize very large and very small eigenvalues [146]: Σz= (1 −λ)Σz+λtrace(Σz) n I n, (6.5) with n the feature dimension, I n the identity matrix of size n , and λ= 0.5 the regularization term. The second part of the optimization requires to estimate the expected classication rate to select the best target. Possible solutions are use crossvalidation or bootstrapping methods to estimate the classication rate using the available data up to time i . However, for small amounts of data, these methods result in estimates with high variance (see also the discussion in [147] about the variance of cross-validation estimators) and have a high computational cost. Alternatively, we propose another approximation of the expected estimation error: the Bhattacharyya coecient. This coecient has been related to the classication error of Gaussian models [148]. Although there is no analytical relation between the coecient and the classication rate, it is possible to derive bounds and good empirical approximations [149]. Under the Gaussian model, the Bhattacharyya coecient ρ between the distributions of the signals, associated to each meaning, N(µc,Σc),N(µw,Σw) , is simply ρ=e−DB where DB is the Bhattacharyya distance: DB(θ) = 1 8(µc−µw)t(Σc+ Σw)−1(µc−µw) + 1 2ln det(Σc+ Σw) √detΣcΣw. (6.6) 84 Chapter 6. Zero-calibration BMI using error-related potentials 6.2.3 Estimation of Task and Online Re-Estimation of Signal Model The use of the Bhattacharyya coecient provides a simple and ecient way to estimate the target when the number of examples is small and avoids the problems of cross-validation. However, when changing to a new target, the signal model does not change and does not have to be re-learned from scratch. Indeed, once the system has correctly identied a task, it has access to correctly labeled data and so it is possible to train a classier as in usual BCI calibration approaches. To achieve this, we factorize the joint distribution p(ξ, θ |Di)∝p(ξ|Di)p(θ|ξ, Di) =p(ei|ξ, (s, a)i, Di−1)p(ξ|Di−1)p(θ|ξ, Di), (6.7) where Di contains the triplets (ei, si, ai) up to time i associated to its labels. The factorization makes explicit that given the target the distribution p(θ| ξ, Di) can be easily evaluated using the labels of each target. We approximate this posterior using the maximum likelihood point estimate ˆ θML ξ per target. For the term p(ξ|Di) , we use a recursive Bayes lter p(ξ=t, |Di)∝p(ei|ξ=t, (s, a)i, Di−1)p(t|Di−1) ≈p(ek|ˆ θML ξ=t)p(ξ=t|Dk−1). (6.8) Notice that we are keeping a dierent symbol model ˆ θML Dξt i−1 for each possible target ξt , the maximum likelihood estimation needs to be done in relation to a dataset Dξt i−1 corresponding to the expected labels of target ξ up to time i−1 . Whenever a task is identied, its labels are transferred to all the triplets of the other tasks to correct the prior for the next step with the right labels. This scheme performs a long term adaptation of θ to accommodate slight variations of EEG such as non-stationarities or variations induced by the task. 6.2.4 Action planning The previous algorithm is able to acquire the needed knowledge but it is not explicit on the process it uses to acquire the data. The goal of the device is to fulll the task desired by the user but, as it has to simultaneously estimate which is the task, it has also to explore regions that allow to disambiguate among dierent tasks. There are several ecient model-based reinforcement learning exploration methods that besides using the task reward function to 6.2. BCI-feedback Based Control Without Explicit Calibration 85 plan for actions, add an additional exploration bonus for states that might provide more learning gains. Several theoretical results show that these approaches allow to learn tasks eciently [150,151]. We dene an uncertainty measure and use model-based planning to select sequences of actions that guide the agent to states that better identify the desired task. The proposed planning method is based on reducing the expected prediction error of the label corresponding to each brain signal. If we choose an action in a given state, given each possible task, we can predict to receive a correct or wrong label. Such label is linked to a signal generation model which dier from task to task, and to the optimal action at that particular state. A state-action where the optimal actions and the signal model is the same for all the hypothesis will be less informative that any other stateaction where either actions or models dier. Thus, our measure of global uncertainty U(s, a) will be higher when, for a given state-action there is an high incongruity between either actions or signal models. For this we compute a similarity matrix S where each element Sij(s, a) corresponds to the similarity of the distributions of the signals of the expected label according to tasks i and j if action a is performed in state s . In order to improve computation eciency we do not rely on a precise metric between distribution and only consider the similarity between the mean of the distribution (empirical tests did no show a signicant impact on the results). The nal uncertainty value U(s, a) can be computed as the sum of the upper diagonal of the similarity matrix. This measure is then used as a classical exploration bonus method by summing the task reward and the uncertainty measure. As the world dynamics is known we can plan actions that maximize the surrogate reward function R(s, a) = (1 −β)Rtask(s, a) + βU(s, a) . We can, for instance, use value iteration and then follow the optimal policy. The reward function R includes a parameter β that provides a gain schedule of exploration versus exploitation. We can either tune this parameter or switch to a pure exploitation approach after reaching the desired condence level on the task model. Interestingly the same approach can be used either in the case of known or unknown signal models. As the uncertainty function combines task and symbol uncertainty, when former is known, the latter becomes the sole source of ambiguity. 6.2.5 Methods for the online experiments The objective of the online evaluation was to determine whether a user was able to teach the virtual device how to reach a target while learning at the same time the classier. For this purpose, we followed the algorithm 86 Chapter 6. Zero-calibration BMI using error-related potentials presented in the previous subsections. For the executed experiments, EEG and EOG signals were recorded following the same conguration as in the previous chapter (see subsection 5.2.1). Four subjects (aged between 25 and 28 ) performed the online experiments. Each subject completed 5 runs of learning from scratch how to reach a randomly chosen target. A run was stopped whenever a task model had more than 99% of the probability of being the correct one, according to Eq. 6.8. During the online execution of the experiments, the device performed actions while the user assessed them. After each action, the device updated its knowledge about the environment and the task following the proposed algorithm (section 6.2.2), and using several features from the user's EEG. These features were extracted from two fronto-central channels (FCz and Cz) within a time window of [200,700] ms (being 0 ms the action onset of the agent) and downsampled to 32 Hz, leading to a vector of length 34 . 6.2.6 Methods for the oine experiments Additionally to the online experiments, further oine experiments were executed. For this analysis, a larger dataset was used to ensure that the results obtained were statistically signicant. Specically, we used the datasets acquired from ten subjects during the experiments described in chapter 3, namely OT1 and OT2. For each subject and each dataset, we simulated twenty runs of 400 steps following the control task by sampling a signal from the data. The objectives of the oine analysis were: (i) study the impact of the exploration approach proposed in section 6.2.4; (ii) evaluate whether the classier learned during the rst phase of the algorithm (section 6.2.2)) has similar decoding capabilities to classiers trained with standard calibration approaches; (iii) the advantages of using a classier rather than the Bhattacharyya coecient during the second phase of the algorithm; and (iv) the number of tasks (targets) that can be learned depending on the dataset tenfold accuracy. 6.3 Results 6.3.1 Online experiments The rst result obtained was that, for all the subjects, the device was able to reach the correct target while learning the signal model at the same time. Table 6.1 shows for each subject and run the number of steps needed to 6.3. Results 87 Table 6.1: Number of steps needed to reach the correct target. Run 1 Run 2 Run 3 Run 4 Run 5 mean ± std Subject 1 95 62 56 60 64 67 ± 16 Subject 2 89 77 98 60 62 77 ± 17 Subject 3 68 80 118 76 157 100 ± 37 Subject 4 98 142 57 142 47 97 ± 45 learn the correct task. Despite there was no sucient information to perform statistical tests, no substantial dierences were found among subjects or runs. On average, the number of steps needed to reach the task was of 85 ± 32. Notice that, in the previous chapter, the time needed to reach a target was 23 steps on average. However, the previous protocol required of a calibration phase of 400 actions on average, whereas in this approach the signal model was also learned during operation. Figure 6.2 shows the evolution of the probability of the correct task, averaged for all subjects and runs (Fig. 6.2a), and separated for each subject and run (Fig. 6.2b). Despite on average the probability function is increasing continuous, there were some runs were the probability had a more erratic behavior, especially on subjects 3 and 4. This indicated that, initially, there were other tasks which were increasing their probability. Thus, the algorithm needed more time to converge to the correct task. 0 50 100 150 200 0 0.2 0.4 0.6 0.8 1 Number of steps Probability of Correct Task 0 50 100 150 200 0 0.2 0.4 0.6 0.8 1 Number of steps Probability of Correct Task s1 s2 s3 s4 (a) (b) Figure 6.2: Evolution of the probability of the correct task (a) averaged across subjects and runs (the shadowed area indicates the standard deviation); and (b) for each separate subject and run. 88 Chapter 6. Zero-calibration BMI using error-related potentials 6.3.2 Oine experiments Planning Methods : We compared the average number of steps (with maximum values of 400 steps) needed to identify the rst task (section 6.2.2) when learning from scratch with dierent planning methods. Figure 6.3 shows the results averaged across subjects, runs and the two datasets. Notice that the large standard deviations are mainly due to the large variations in classication accuracy across subjects and datasets. Several planning methods were able to correctly estimate the tasks. On the other hand, a greedy exploration (i.e. always trying to follow the most probable task) does not allow the system to explore suciently, and at least some random exploration is necessary to allow a correct identication of the task (e.g. ε -greedy). The proposed planning method based on expected signal uncertainty (section 6.2.4), leads the system to regions that improve disambiguation among hypotheses in a faster way. Compared to assessing uncertainty only on the meaning space, the proposed algorithm performs better as such methods do not take into account the signal to meaning uncertainty inherent to the problem to solve. Online re-estimation of classier : After identifying the rst task, the second phase of the algorithm was training and adapting a classier (see section 6.2.3). The quality of the classier learned with the proposed method can be measured according to the percentage of labels correctly learned (according to the ground truth label), see Figure 6.4. Notice how this percentage varied with the ten-fold accuracy of a given subject and dataset. In general, having accuracies higher than 75% guaranteed that more than 90% of the labels were correctly learned. This result shows that our algorithm could also be used to collect training data for calibrating any other state-of-the-art error potentials classier, but has the important advantage of controlling the device at the same time. Figure 6.5a demonstrate the advantage of training a classier after the rst reached task instead of keeping the estimation given by the Bhattacharyya coecient. Indeed, the Bhattacharrya coecient works very well for small amounts of data because it directly compares model parameters. On the other hand, when there are sucient data, training a classier allows for a faster adaptation since the classier makes a much harder decision when evaluating a new EEG signal. Finally, gure 6.5b shows the number of tasks identied with respect to the accuracy of the dataset, and the number of tasks wrongly identied. Notice how the number of identied task is correlated to the quality of the dataset. Importantly, the algorithm was able to identify on average 20 tasks in 400 steps on average without the need for a calibration procedure. 7.1. Future work 95 process, its elicitation will depend on the subjective evaluation of each user. Thus, an important question is to evaluate whether there is any measurable assessment information in the user's brain during continuous movements performed by a device; and the feasibility of detecting these signals online. Our preliminary works with two subjects have shown that these potentials do exist during continuous motions of a virtual device, and we can detect them on single trial using frequency information [152]. Furthermore, these ndings have also been studied in a realistic task with a mobile robot, where preliminary results suggested that it is feasible to teach devices under continuous spaces following the proposed BMI paradigm [88]. Despite this thesis presented a deep analysis of the error potentials under constrained conditions, it is still unclear what are the mechanisms underlying the generation of these signals. In fact, there is still a debate on the eld of psychophysiology about whether the error-related negativity encodes a quantitative reward predictive error signal (RPE) or valence [41] (the dierence between expected and received outcome), or rather the absolute value of this RPE [35]. According to some authors, the latter would imply that the ERN actually encodes the surprisingness of an event. In the specic case of ErrPs, several works have shown that this signal is also present under situations where the percentage of errors is very high ( 40% −50% ) but with a decrease in amplitude [47,153]. In this thesis we showed that, during the beginning of any run of the teaching paradigm, there is a high percentage of errors (up to 70% ) due to the device exploring the task space. Even so, the ErrP was present and we were able to detect it online. After a certain number of actions, this percentage of errors decreased drastically to 10% when the device learns the policy. Nonetheless, what would happen under situations where there is always a large number of errors for a large amount of time? It could happen that the user nds the errors the expected outcome, being the correct task execution an unexpected outcome. Some of these situations are present when performing very demanding tasks, e.g. hitting a ball with a baseball bat. Thus, designing and performing experiments would be necessary to address this issue. Even after understanding the relationship between surprise, expectance and error potentials, there is still another open question linking error potentials and cognitive neuroscience: To what extent are the error-related potentials the same phenomenon that the error-related negativity? As stated in the introduction, the error-related potentials are characterized by three main peaks: N2, P3 and N4, whereas the ERN (and the FRN) is a single component appearing at around the same timing that the N2. There has been works comparing the N2, the FRN and ERN under similar experimental conditions [111], but what about the other ErrP components? In order 96 Chapter 7. Conclusions to answer this question, a reasonable option would be to perform experiments with simultaneous recordings of EEG and fMRI [154] to determine the underlying generation of the ErrPs. In fact, a recent work showed the eectiveness of this experiment to understand the brain mechanisms generating the feedback-related negativity [35]. The presented paradigm strongly relies on the fact that a feedback signal can be extracted from the brain. In the instances presented throughout the thesis, this feedback encoded binary information about the actions wellness. However, no study was performed to understand up to what point the feedback given by the users encode the wellness about long term goals. Recently, Knox et al. suggested that the feedback given by humans could actually be encoded as the expected return in the long term, rather than a binary reward of the current action performed by a device [155157]. In principle, this nding should be the case also for the error-related potentials. However, more studies would be needed to address this question. From the point of view of brain-machine interfaces, the future work is associated to further study the proposed paradigm. This framework has shown how the users are able to teach a device by giving feedback. Thus, the device can adapt to the user's preferences. But, at the same time, the user is implicitly adapting to the device since this device provides a visual feedback of the task to the user. This term is usually known as co-adaptation between user and machine, and it is well known in the non-invasive and invasive BMI community [72,100,145], demonstrating that co-adaptation improves the nal decoding performance that can be obtained with the BMI decoder. Despite a co-adaptation seemed to be present in the proposed paradigm, this characteristic has not be explicitly addressed. In which sense and how the user adapts to the machine (e.g. how the ErrPs vary in shape throughout time) is still an open question that needs to be studied. In the introduction, a broad denition of the proposed paradigm was given (i.e. use cognitive signals and plug them into intelligent device controllers). On the contrary, throughout this thesis only a specic instance of this paradigm (error potentials + reinforcement learning) has been presented. Thus, the proposed paradigm has room for many improvements as there are many other instances of this paradigm that could be either work independently or together with the proposed instance. In non-invasive BMIs, there are several cognitive signals that could help in improving the paradigm, ranging from anticipatory signals encoding goal direction or movement intention [22]; to user's (c)overt attention to increase or decrease the device learning rate [113, 158]. For instance, thanks to the combination of goaloriented cognitive signals (such as the target or hand direction), it would be possible to use the proposed system as a complementary paradigm tting the 7.1. Future work 97 approach into a two-step system: rst decide where to go, then how to go. On the other hand, in principle this paradigm could also work using invasive BMIs, where richer information can be extracted from the brain. Several works have already reported that it is possible to measure single-neuron ring rates associated to error processing, and encoding it as rewards [101]. The use of invasive approaches would allow to obtain more informative rewards, that could be use to nely adapt the device to the user's preferences. Finally, this paradigm was only tested in laboratory conditions. Thus, the natural next step of this thesis would be to apply the paradigm in realistic experiments so as to demonstrate its feasibility to work in out-of-the-lab scenarios. Indeed, the last two chapters solved two key problems of the paradigm: the calibration phase and the large amount of learning time. Thanks to these improvements, several experiments could be designed involving daily-life activities, such as drinking from a glass using a neuroprosthetic device. In the given example, this robotic device would learn from a set of optimal policies and adapt to the user's preferences without the need of any calibration phase. Once these experiments were designed, the proposed approach would be one step nearer to restore or replace the motor execution of limbs as a way of improving the daily life with real patients. Thus, performing experiments with patients would be the last and more demanding improvement of the current thesis. 98 Chapter 7. Conclusions Bibliography [1] Wolpaw JR, Birbaumer N, McFarland DJ, Pfurtscheller G, and Vaughan TM. Brain-computer interfaces for communication and control. Clin Neurophy , 113(6):76791, June 2002. [2] Millán JdR, Rupp R, Müller-Putz GR, Murray-Smith R, Giugliemma C, Tangermann M, Vidaurre C, Cincotti F, Kübler A, Leeb R, Neuper C, Müller KR, and Mattia D. Combining braincomputer interfaces and assistive technologies: State-of-the-art and challenges. Front Neurosci , 4, 2010. [3] Lebedev MA and Nicolelis MAL. Brain-machine interfaces: past, present and future. Trends Neurosci , 29(9):536546, 2006. [4] Carmena JM, Lebedev MA, Crist RE, O'Doherty JE, Santucci DM, Dimitrov DF, Patil PG, Henriquez CS, and Nicolelis MAL. Learning to control a brain-machine interface for reaching and grasping by primates. PLoS Biol , 1:193208, 2003. [5] Millán JdR, Renkens F, Mouri no J, and Gerstner W. Noninvasive brain-actuated control of a mobile robot by human EEG. IEEE Trans Biomed Eng , 51:10261033, 2004. [6] Musallam S, Corneil BD, Greger B, Scherberger H, and Andersen RA. Cognitive control signals for neural prosthetics. Science , 305:162163, 2004. [7] Müller-Putz GR, Scherer R, Pfurtscheller G, and Rupp R. EEG-based neuroprosthesis control: A step towards clinical practice. Neurosci Lett , 382:169174, 2005. [8] Velliste M, Perel S, Spalding MC, Whitford AS, and Schwartz AB. Cortical control of a prosthetic arm for self-feeding. Nature , 453:1098 1101, 2008. 99 100 Bibliography [9] Iturrate I, Antelis J, Kübler A, and Minguez J. Non-invasive brainactuated wheelchair based on a P300 neurophysiological protocol and automated navigation. IEEE Trans Robot , 25:614627, 2009. [10] Tonin L, Carlson T, Leeb R, and JdR Millán. Brain-controlled telepresence robot by motor-disabled people. In Conf Proc IEEE Eng Med Biol Soc , 2011. [11] Ethier C, Oby ER, Bauman MJ, and Miller LE. Restoration of grasp following paralysis through brain-controlled stimulation of muscles. Nature , 485:368371, 2012. [12] Hochberg LR, Bacher D, Jarosiewicz B, Masse NY, Simeral JD, Vogel J, Haddadin S, Liu J, Cash S, van der Smagt P, and Donoghue JP. Reach and grasp by people with tetraplegia using a neurally controlled robotic arm. Nature , 485:372375, 2012. [13] Carlson T and Millán JdR. Braincontrolled wheelchairs: A robotic architecture. IEEE Robot Autom Mag , 20:6573, 2013. [14] Collinger JL, Wodlinger B, Downey JE, Wang W, Tyler-Kabara EC, Weber DJ, McMorland AJC, Velliste M, Boninger ML, and Schwartz AB. High-performance neuroprosthetic control by an individual with tetraplegia. The Lancet , 381:557564, 2013. [15] Wolpaw JR and McFarland DJ. Control of a two-dimensional movement signal by a noninvasive brain-computer interface in humans. In Proc Natl Acad Sci USA , volume 101, pages 1784917854, 2004. [16] Santhanam G, Ryu SI, Yu BM, Afshar A, and Shenoy KV. A highperformance brain-computer interface. Nature , 442:195198, 2006. [17] Bradberry TJ, Gentili RJ, and Contreras-Vidal JL. Reconstructing three-dimensional hand movements from noninvasive electroencephalographic signals. J Neurosci , 30:34323437, 2010. [18] O'Doherty JE, Lebedev MA, It PJ, Zhuang KZ, Shokur S, Bleuler H, and Nicolelis MAL. Active tactile exploration using a brain-machinebrain interface. Nature , 479:228231, 2011. [19] Hauschild M, Mulliken GH, Fineman I, Loeb GE, and Andersen RA. Cognitive signals for brain-machine interfaces in posterior parietal cortex include continuous 3D trajectory commands. Proc Natl Acad Sci USA , 109:1707517080, 2012. Bibliography 101 [20] Ball T, Schulze-Bonhage A, Aertsen A, and Mehring C. Dierential representation of arm movement direction in relation to cortical anatomy and function. J Neural Eng , 6:016006, 2009. [21] Fried I, Mukamel R, and Kreiman G. Internally generated preactivation of single neurons in human medial frontal cortex predicts volition. Neuron , 69:548562, 2011. [22] Lew E, Chavarriaga R, Silvoni S, and Millán JdR. Detection of selfpaced reaching movement intention from EEG signals. Front Neuroeng , 5:13:doi: 10.3389/fneng.2012.00013, 2012. [23] Kim HK, Biggs J, Schloerb W, Carmena J, Lebedev MA, Nicolelis MAL, and Srinivasan MA. Continuous shared control for stabilizing reaching and grasping with brain-machine interfaces. IEEE Trans Biomed Eng , 53(6):11641173, 2006. [24] Scott SH. Optimal feedback control and the neural basis of volitional motor control. Nat Rev Neurosci , 5:534546, 2004. [25] Dimitrijevic MR, Gerasimenko Y, and Pinter MM. Evidence for a spinal central pattern generator in humans. Annals of the New York Academy of Sciences , 860:36076, November 1998. [26] Courtine G, Gerasimenko Y, van den Brand R, Yew A, Musienko P, Zhong H, Song B, Ao Y, Ichiyama RM, Lavrov I, Roy RR, Sofroniew MV, and Edgerton VR. Transformation of nonfunctional spinal circuits into functional states after the loss of brain input. Nat Neurosci , 12:13331342, 2009. [27] Takakusaki K, Saitoh K, Harada H, and Kashiwayanagi M. Role of basal ganglia-brainstem pathways in the control of motor behaviors. Neuroscience research , 50(2):13751, October 2004. [28] Kalaska JF, Scott SH, Cisek P, and Sergio LE. Cortical control of reaching movements. Current opinion in neurobiology , 7(6):84959, December 1997. [29] Rizzolatti G and Luppino G. The cortical motor system. Neuron , 31:889901, 2001. [30] Waldert S, Preissl H, Demandt E, Braun C, Birbaumer N, Aertsen A, and Mehring C. Hand movement direction decoded from MEG and EEG. J Neurosci , 28:1000108, 2008. 102 Bibliography [31] Holroyd CB and Coles MGH. The neural basis of human error processing: Reinforcement learning, dopamine, and the error-related negativity. Psychol Rev , 109:679709, 2002. [32] Nieuwenhuis S, Holroyd CB, Mola N, and Coles MGH. Reinforcementrelated brain potentials from medial frontal cortex: origins and functional signicance. Neuroscience and Biobehavioral Reviews , 28:441 448, 2004. [33] Frank MJ, Woroch BS, and Curran T. Error-related negativity predicts reinforcement learning and conict biases. Neuron , 47(4):495501, August 2005. [34] Holroyd CB, Krigolson OE, and Lee S. The impact of deliberative strategy dissociates ERP components related to conict processing vs. reinforcement learning. Frontiers in neuroscience , 6(April):43, January 2012. [35] Hauser TU, Iannaccone R, Stämpi P, Drechsler R, Brandeis D, WalitzaS, and Brem S. The feedback-related negativity (FRN) revisited: New insights into the localization, meaning and network organization. NeuroImage (to be published) , August 2013. [36] Allman JM, Hakeem A, Erwin JM, Nimchinsky E, and Hof P. The anterior cingulate cortex. Annals of the New York Academy of Sciences , 935(1):107117, 2001. [37] Schultz W. Predictive reward signal of dopamine neurons. Journal of Neurophysiology , 80(1):127, July 1998. [38] Bayer HM and Glimcher PW. Midbrain dopamine neurons encode a quantitative reward prediction error signal. Neuron , 47(1):12941, July 2005. [39] Niv Y, Du MO, and Dayan P. Dopamine, uncertainty and TD learning. Behavioral and Brain Functions , 1(6):19, 2005. [40] Glimcher PW. Understanding dopamine and reinforcement learning: The dopamine reward prediction error hypothesis. Proc Natl Acad Sci USA , 108(Supplement 3):1564715654, 2011. [41] Luo Q and Qu C. Comparison enhances size sensitivity: neural correlates of outcome magnitude processing. PloS one , 8(8):e71186, January 2013. Bibliography 103 [42] Falkenstein M, Hoormann J, Christ S, and Hohnsbein J. ERP components on reaction errors and their functional signicance: A tutorial. Biol Psychol , 51:87107, 2000. [43] Miltner WHR, Braun CH, and Coles MGH. Event-related brain potentials following incorrect feedback in a time-estimation task: Evidence for a generic neural system for error detection. Journal of Cognitive Neuroscience , 9(6):788798, 1997. [44] van Schie HT, Mars RB, Coles MGH, and Bekkering H. Modulation of activity in medial frontal and motor cortices during error observation. Neural Netw , 7:549554, 2004. [45] Schalk G, Wolpaw JR, McFarland DJ, and Pfurtscheller G. EEG-based communication: Presence of an error potential. Clin Neurophysiol , 111(12):21382144, 2000. [46] Ferrez PW and Millán JdR. Error-related EEG potentials generated during simulated brain-computer interaction. IEEE Trans Biomed Eng , 55:923929, 2008. [47] Chavarriaga R and Millán JdR. Learning from EEG error-related potentials in noninvasive brain-computer interfaces. IEEE Trans Neural Syst Rehabil Eng , 18(4):381388, 2010. [48] Duncan CC, Barry RJ, Connolly JF, Fischer C, Michie PT, Näätänen R, Polich J, Reinvang I, and van Petten C. Event-related potentials in clinical research: Guidelines for eliciting, recording, and quantifying mismatch negativity, P300, and N400. Clin Neurophysiol , 120:1883 1908, Sep 2009. [49] Olofsson JK, Nordin S, Sequeira H, and Polich J. Aective picture processing: An integrative review of ERP ndings. Biol Psychol , 77(3):24765, March 2008. [50] Luck SJ. An introduction to the event-related potential technique . The MIT Press, 2005. [51] Kutas M, McCarthy G, and Donchin E. Augmenting mental chronometry: The P300 as a measure of stimulus evaluation time. Science , 197(4305):792, 1977. [52] Butteld A, Ferrez PW, and Millan JdR. Towards a robust BCI: Error potentials and online learning. Neural Systems and Rehabilitation Engineering, IEEE Transactions on , 14(2):164168, 2006. 104 Bibliography [53] Blumberg J, Rickert J, Waldert S, Schulze-Bonhage A, Aertsen A, and Mehring C. Adaptive classication for brain computer interfaces. In Conf Proc IEEE Eng Med Biol Soc , pages 25369, January 2007. [54] Zeyl TJ and Chau T. A case study of linear classiers adapted using imperfect labels derived from human event-related potentials. Pattern Recognition Letters , May 2013. [55] Blankertz B, Dornhege G, Schafer C, Krepki R, Kohlmorgen J, Muller KR, Kunzmann V, Losch F, and Curio G. Boosting bit rates and error detection for the classication of fast-paced motor commands based on single-trial EEG analysis. IEEE Trans Neural Syst Rehabil Eng , 11(2):127131, 2003. [56] Ferrez PW and Millán JdR. Simultaneous real-time detection of motor imagery and error-related potentials for improved BCI accuracy. In Proc 4th Int BCI Workshop and Train Course , 2008. [57] Dal Seno B, Matteucci M, and Mainardi L. Online detection of P300 and error potentials in a BCI speller. Comput Intell Neurosci , 2010:307254, 2010. [58] Schmidt NM, Blankertz B, and Treder MS. Online detection of errorrelated potentials boosts the performance of mental typewriters. BMC Neurosci , 13(1):19, 2012. [59] Llera A, Gómez V, and Kappen HJ. Adaptive classication on braincomputer interfaces using reinforcement signals. Neural Computation , 24(11):29002923, 2012. [60] Spüler M, Bensch M, Kleih S, Rosenstiel W, Bogdan M, and Kübler A. Online use of error-related potentials in healthy users and people with severe motor impairment increases performance of a P300-BCI. Clin Neurophysiol , 123(7):132837, July 2012. [61] Sutton RS and Barto AG. Reinforcement learning: An introduction . MIT Press, 1998. [62] Sutton RS. Introduction: The challenge of reinforcement learning. In Reinforcement learning: An introduction , pages 13. Springer, 1992. [63] Sutton RS, McAllester D, Singh S, and Mansour Y. Policy gradient methods for reinforcement learning with function approximation. In Advances in Neural Information Processing Systems (NIPS) , volume 12, pages 10571063, 1999. Bibliography 111 [131] Thompson DE, Warschausky S, and Huggins JE. Classier-based latency estimation: A novel way to estimate and predict BCI accuracy. Journal of neural engineering , 10(1):016006, December 2012. [132] Wilson JA, Mellinger J, Schalk G, and Williams J. A procedure for measuring latencies in brain-computer interfaces. IEEE Transactions on Biomedical Engineering , 57(7):178597, July 2010. [133] Sauser E. Robottoolkit. Available online: lasa.ep.ch/RobotToolKit. [134] Woody CD and Nahvi MJ. Application of optimum linear lter theory to the detection of cortical signals preceding facial movement in cat. Exp Brain Res , 16:455465, 1973. [135] Levine SP, Huggins JE, BeMent SL, Kushwaha RK, Schuh LA, Rohde MM, Passaro EA, Ross DA, Elisevich KV, and Smith BJ. A direct brain interface based on event-related potentials. IEEE Trans Rehabil Eng , 8(2):180, 2000. [136] Benjamini Y and Yekutieli D. The control of the false discovery rate in multiple testing under dependency. Ann Stat , pages 11651188, 2001. [137] Jung TP, Makeig S, Westereld M, Townsend J, Courchesne E, and Sejnowski TJ. Analysis and visualization of single-trial event-related potentials. Hum Brain Mapp , 14(3):166185, 2001. [138] Gerson AD, Parra LC, and Sajda P. Cortically-coupled computer vision for rapid image search. IEEE Trans Neural Syst Rehabil Eng , 14(2):174179, Jun 2006. [139] Hong B, Guo F, Liu T, Gao X, and Gao S. N200-speller using motiononset visual response. Clin Neurophysiol , 120(9):16581666, 2009. [140] Berndt D and Cliord J. Using dynamic time warping to nd patterns in time series. In AAAI Workshop on Knowledge Discovery in Databases , volume 398, pages 229248, 1994. [141] Casarotto S, Bianchi AM, Cerutti S, and Chiarenza GA. Dynamic time warping in the analysis of event-related potentials. IEEE Eng Med Biol Mag , 24(1):6877, 2005. [142] Croft RJ and Barry RJ. EOG correction of blinks with saccade coecients: a test and revision of the aligned-artefact average solution. Clinical neurophysiology , 111(3):44451, March 2000. 112 Bibliography [143] Schlögl A, Keinrath C, Zimmermann D, Scherer R, Leeb R, and Pfurtscheller G. A fully automated correction method of EOG artifacts in EEG recordings. Clinical neurophysiology , 118(1):98104, January 2007. [144] Bishop CM. Pattern recognition and machine learning , chapter Linear models for classication. Springer, 2006. [145] Orsborn AL, Dangi S, Moorman HG, and Carmena J. Closed-loop decoder adaptation on intermediate time-scales facilitates rapid BMI performance improvements independent of decoder initialization conditions. IEEE Transactions on neural systems and rehabilitation engineering , 20(4), 2012. [146] Friedman JH. Regularized discriminant analysis. Journal of the American statistical association , 84(405):165175, 1989. [147] Bengio Y and Grandvalet Y. No unbiased estimator of the variance of k-fold cross-validation. The Journal of Machine Learning Research , 5:10891105, 2004. [148] Kailath T. The divergence and bhattacharyya distance measures in signal selection. IEEE Trans. Commun. Technol. , 15(3):5260, 1967. [149] Chulhee L and Euisun C. Bayes error evaluation of the gaussian ml classier. Geoscience and Remote Sensing, IEEE Transactions on , 38(3):14711475, 2000. [150] Brafman RI and Tennenholtz M. R-max-a general polynomial time algorithm for near-optimal reinforcement learning. Journal of Machine Learning Research , 3, 2003. [151] Kolter JZ and Ng AY. Near-bayesian exploration in polynomial time. In International Conference on Machine Learning . ACM, 2009. [152] Iturrate I, Omedes J, and Montesano L. Detection of event-less error related potentials. In IROS 2013 Workshop on Neuroscience and Robotics , 2013. [153] Ferrez PW. Error-Related EEG Potentials in Brain-Computer Interfaces . PhD thesis, École Polytechnique Fédérale de Laussane, 2007. [154] Rosa MJ, Daunizeau J, and Friston KJ. EEG-fMRI integration: a critical review of biophysical modeling and data analysis approaches. Journal of Integrative Neuroscience , 09(04):453476, December 2010. Bibliography 113 [155] Knox WB and Stone P. Interactively shaping agents via human reinforcement: The tamer framework. In Proceedings of the fth international conference on Knowledge capture , pages 916. ACM, 2009. [156] Knox WB and Stone P. Combining manual feedback with subsequent mdp reward signals for reinforcement learning. In Proceedings of the 9th International Conference on Autonomous Agents and Multiagent Systems (AAMAS'10) , pages 512, 2010. [157] Knox WB, Glass B, Love B, Maddox W, and Stone P. How humans teach agents: A new experimental perspective. International Journal of Social Robotics, Special Issue on Robot Learning from Demonstration , 2012. [158] Tonin L, Leeb R, and Millán JdR. Time-dependent approach for single trial classication of covert visuospatial attentions. Journal of Neural Engineering , 9(4):e045011, 2012.