Tiziana Cattai

Fabrizio De

Vico Fallani

Stefania Colonnese

Ville De

associate professor, HDR, rapporteur -Pinna Dimitri

associate professor, HDR, examinateur -Scarano Andrea

Gaetano

Leveraging brain connectivity networks to detect mental states during motor imagery

Keywords: Functional connectivity, signal processing, graph signal processing, brain-computer interface

The brain is a complex network and we know that inter-areal synchronization and de-synchronization mechanisms are crucial to perform motor and cognitive tasks.

Nowadays, brain functional interactions are studied in brain-computer interface (BCI) applications with more and more interest. This might have strong impact on BCI systems, typically based on univariate features which separately characterize brain regional activities. Indeed, brain connectivity features can be used to develop alternative BCIs in an effort to improve performance and to extend their real-life applicability. The ambition of this thesis is the investigation of brain functional connectivity networks during motor imagery (MI)-based BCI tasks. It aims to identify complex brain functioning, re-organization processes and timevarying dynamics, at both group and individual level. This thesis presents different developments that sequentially enrich an initially simple model in order to obtain a robust method for the study of functional connectivity networks. Experimental results on simulated and real EEG data recorded during BCI tasks prove that our proposed method well explains the variegate behaviour of brain EEG data. Specifically, it provides a characterization of brain functional mechanisms at group level, together with a measure of the separability of mental conditions at individual level. We also present a graph denoising procedure to filter data which simultaneously preserve the graph connectivity structure and enhance the signal-to-noise ratio. Since the use of a BCI system requires a dynamic interaction between user and machine, we finally propose a method to capture the evolution of time-varying data. In essence, this thesis presents a novel framework to grasp the complexity of graph functional connectivity during cognitive tasks.

List of Figures

2.1 Graphical representation of the most common functional connectivity estimators used in literature. They are organized in a coloured structure to identify different properties. We report undirected and directed estimators in green and red colours respectively. Then, they can be represented according to their ability to capture bivariate or multivariate interactions, in blue or green boxes. The characterization of possible non linearities is reported in grey boxes. . . . . . . . We plot the J-divergence over M = 9, 1s long, time intervals with 50% overlapping, versus the time interval index. The J-divergence is computed in β band and averaged across subjects. . . . . . . . . . we can reach an higher level. In fact, after 500ms, we have the first Laplacian estimation and the first classification can be performed.

General

After the second interval (i.e. at 1s) the first L1-PCs can be computed. We have to wait until 2s (i.e. the 4th interval) to collect all the samples to complete the first layer with k = 1 and to perform the second L1-PCA, reaching the second layer. At 4s, we have all the samples needed to compute the last L1-PCA. . . . . . . . . . . . 83

List of Tables

Introduction

Functional connectivity (FC), measuring statistical dependences between different brain signals, captures crucial information to investigate brain functioning during cognitive tasks. Indeed, we currently interpret the brain as a complex network [START_REF] Bassett | Network neuroscience[END_REF][START_REF] Sporns | Graph theory methods: applications in brain networks[END_REF], where neurons belonging to different regions collaborate to integrate information through anatomical and functional connections [START_REF] Bullmore | Complex brain networks: graph theoretical analysis of structural and functional systems[END_REF]. The possibility to explore brain interaction mechanisms in cognitive tasks, such as those involved during brain-computer interface (BCI) experiments, is gaining more and more interest [START_REF] Vico Fallani | Network neuroscience for optimizing brain-computer interfaces[END_REF][START_REF] Gonzalez-Astudillo | Network-based brain computer interfaces: principles and applications[END_REF].

This thesis aims to provide new insight by adopting a signal processing and graph signal perspective for the application of functional brain connectivity in the BCI context. In particular, we intend to study the importance of brain interactions in the characterization of brain states. We specifically focus on the study of motor imagery-based BCIs.

BCIs are communication systems in which the interaction between the subject and the external world is realized without the peripheral neuro-muscolar activity [START_REF] Wolpaw | Brain-computer interfaces: Definitions and principles[END_REF][START_REF] Wolpaw | Brain-computer interfaces for communication and control[END_REF][START_REF] Wolpaw | The wadsworth center brain-computer interface (bci) research and development program[END_REF]. BCI systems are based on the identification of mental states from brain signals, and because they depend on the subject's ability to voluntarily modulate the brain activity, they are being increasingly investigated for control, communication and rehabilitation [START_REF] Wolpaw | Brain-computer interfaces: principles and practice[END_REF]. Unfortunately, the real-life applicability of BCIs is still limited. It is estimated that a significant percentage of users (around 30%) cannot correctly interact with the interface (i.e. reaching the accuracy level needed to control BCIs) because BCIs are not simple systems and even trained subjects might not be able to voluntarily modulate their brain activity. This phenomenon is generally indicated as BCI-inefficiency [START_REF] Thompson | Critiquing the concept of bci illiteracy[END_REF][START_REF] Vidaurre | Towards a cure for bci illiteracy[END_REF] and it globally limits the usability of BCI systems in clinical applications [START_REF] Chaudhary | Neuropsychological and neurophysiological aspects of brain-computer-interface (bci) control in paralysis[END_REF][START_REF] Soekadar | Brain-machine interfaces in neurorehabilitation of stroke[END_REF].

In the last decade, great effort has been made to improve of BCI accuracy. One 2 possibility consists in finding the best mental strategy to identify the user's intent, looking for the most suitable feedback to efficiently interact with the subject [START_REF] Jeunet | Using recent bci literature to deepen our understanding of clinical neurofeedback: A short review[END_REF][START_REF] Roc | A review of user training methods in brain computer interfaces based on mental tasks[END_REF][START_REF] Schumacher | Towards explanatory feedback for user training in brain-computer interfaces[END_REF]. Another possibility concerns the development of advanced classification algorithms to improve the detection of the user intention, using refined signal processing tools to improve signals' characterization [START_REF] Lotte | A review of classification algorithms for eeg-based brain-computer interfaces[END_REF][START_REF] Lotte | A review of classification algorithms for eeg-based brain-computer interfaces: a 10 year update[END_REF]. These methods are able to enhance BCIs performances, but they are blind to the brain mechanisms underlying the classification and they do not allow a clear neuro-physiological interpretation. This aspect has a crucial importance in most clinical applications, where the brain's structure or functioning is damaged and other solutions able to take it into account are needed.

This thesis aims to present a possible solution. It consists in looking for alternative features, able to characterize brain functional mechanisms involved during complex tasks, such as motor imagery [START_REF] Vico Fallani | Network neuroscience for optimizing brain-computer interfaces[END_REF][START_REF] Gonzalez-Astudillo | Network-based brain computer interfaces: principles and applications[END_REF]. Functional connectivity (FC) describes interactions between spatially distributed areas by evaluating functional dependencies between activities localized in different regions [START_REF] Bastos | A tutorial review of functional connectivity analysis methods and their interpretational pitfalls[END_REF]. Those features appear more representative of brain processes than univariate ones, such as power spectral density [START_REF] Jin | Bispectrumbased channel selection for motor imagery based brain-computer interfacing[END_REF][START_REF] La Rocca | Human brain distinctiveness based on eeg spectral coherence connectivity[END_REF][START_REF] Park | Frequency-optimized local region common spatial pattern approach for motor imagery classification[END_REF].

Our first question is to investigate the actual possibility of applying functional connectivity and basic graph estimators to characterize mental states. To answer this question we consider two state-of-the-art FC estimators, i.e. the spectralcoherence and imaginary-coherence [START_REF] Bastos | A tutorial review of functional connectivity analysis methods and their interpretational pitfalls[END_REF][START_REF] Gonzalez-Astudillo | Network-based brain computer interfaces: principles and applications[END_REF][START_REF] Hamedi | Electroencephalographic motor imagery brain connectivity analysis for bci: a review[END_REF]. We find that integrating networkbased features enables a better characterization of brain states and we uncover the presence of a twofold mechanism of amplitude and phase synchronization between brain signals during motor imagery. At the same time, the inclusion of FC-based features in the classification, generally allows to increase the accuracy in the detection of the user's intent, with respect to standard univariate techniques only.

After testing the potential of graph connectivity features in the identification of mental states, we ask what improvement in states' separability can be obtained through robust functional connectivity estimation. To understand it, we use tools from different fields, from neuroscience to signal processing and graph signal processing [START_REF] Vico Fallani | Graph analysis of functional brain networks: practical issues in translational neuroscience[END_REF][START_REF] Ortega | Graph signal processing: Overview, challenges, and applications[END_REF][START_REF] Shuman | The emerging field of signal processing on graphs: Extending high-dimensional data analysis to networks and other irregular domains[END_REF]. We develop an original algorithm to denoise the graph Laplacian in order to obtain a more robust graph FC estimation. This procedure simultaneously preserves the signal-to-noise ratio and the graph connectivity properties. The validity of the denoising algorithm is verified by measuring the distance between Laplacians related to two states. The distance between those mental states is quantified with a novel formulation of the Jensen divergence [START_REF] Basseville | Divergence measures for statistical data processing-an annotated bibliography[END_REF].

From brain signals to functional connectivity and graph analysis

Modern neuroscience has started with the characterization of the human brain as a system of interacting neurons [START_REF] Ramon-Y Cajal | Histology of the nervous system[END_REF]. Indeed, research has demonstrated that the brain is not constituted by isolated systems, but neurons cooperate and interact [START_REF] Bassett | Network neuroscience[END_REF][START_REF] John | The neurophysics of consciousness[END_REF] to perform complex tasks. Today we interpret the brain as a complex network [START_REF] Sporns | Graph theory methods: applications in brain networks[END_REF] and new techniques have been developed to deeply investigate its functioning by merging different fields, such as network theory and graph signal processing [START_REF] Vico Fallani | Graph analysis of functional brain networks: practical issues in translational neuroscience[END_REF][START_REF] Huang | A graph signal processing perspective on functional brain imaging[END_REF][START_REF] Itani | A graph signal processing framework for the classification of temporal brain data[END_REF][START_REF] Medaglia | Brain and cognitive reserve: translation via network control theory[END_REF][START_REF] Pessoa | Understanding brain networks and brain organization[END_REF].

These methodological advances have provided novel insights into fundamental and clinical neuroscience as well as into several applications from neuro-rehabilitation, to disease modelling and therapy [START_REF] Davidson | Application of describing function analysis to a model of deep brain stimulation[END_REF][START_REF] Fornito | The connectomics of brain disorders[END_REF][START_REF] Guggisberg | Brain networks and their relevance for stroke rehabilitation[END_REF][START_REF] Padole | Early detection of alzheimer's disease using graph signal processing on neuroimaging data[END_REF]. This includes the study of BCIs, which have been mainly explored by using univariate brain features [START_REF] Jin | Bispectrumbased channel selection for motor imagery based brain-computer interfacing[END_REF][START_REF] Park | Frequency-optimized local region common spatial pattern approach for motor imagery classification[END_REF]. Indeed, functional connectivity and network features have been recently applied in BCI-applications [START_REF] Gonzalez-Astudillo | Network-based brain computer interfaces: principles and applications[END_REF][START_REF] Hamedi | Electroencephalographic motor imagery brain connectivity analysis for bci: a review[END_REF]. However, state-of-the-art results do not uniquely show the advantage of using those features instead of classical univariate ones [START_REF] Hamedi | Electroencephalographic motor imagery brain connectivity analysis for bci: a review[END_REF].

To facilitate the reader in the comprehension of the advanced methods presented here, we are going to propose in the following a concise overview of the state-of-the-art associated to FC estimation from time series data, graph signal processing theory and models to capture time-varying dynamics. Note, that a more specific bibliography is directly presented in each chapter.

How to estimate brain connectivity from time series

In the last decades, many methods have been proposed to measure functional interactions between brain signals, based on tools coming from several domains,

Chapter 2. From brain signals to functional connectivity and graph analysis 5

from signal processing to information theory. Although each estimator has its own characteristics, the general idea consists in quantifying the functional dependencies between each pair of signals.

The first distinction between FC estimators is between those able to capture mutual interactions and those able to quantify causal information flows [START_REF] Vico Fallani | Graph analysis of functional brain networks: practical issues in translational neuroscience[END_REF]. We commonly refer to them as undirected and directed estimators, respectively. Further distinctions can be done according to the possibility to capture linear or non linear interactions, or to measure bivariate or multivariate connections, or the potential to quantify FC in time or in frequency domain, as well as to describe timevarying dynamics [START_REF] Bastos | A tutorial review of functional connectivity analysis methods and their interpretational pitfalls[END_REF]. In Fig. 2.1, we graphically present some examples of FC estimators, organized in a subjective way. We remark that it is a non-exhaustive list and other distinctions are possible. We report undirected and directed estimators in green and red colours respectively. Then, they can be represented according to their ability to capture bivariate or multivariate interactions, in blue or green boxes. The characterization of possible non linearities is reported in grey boxes.

The choice of the FC estimator can be a tricky task [START_REF] Vico Fallani | Graph analysis of functional brain networks: practical issues in translational neuroscience[END_REF]. Comparative studies between several FC estimators have been proposed, but conclusions are not unambiguous since methods' performances depend both on the estimator itself and data characteristics [START_REF] Bakhshayesh | Detecting synchrony in eeg: A comparative study of functional connectivity measures[END_REF][START_REF] Krakovská | Comparison of six methods for the detection of causality in a bivariate time series[END_REF][START_REF] Wang | A systematic framework for functional connectivity measures[END_REF]. A first possible approach consists in applying all

How to estimate brain connectivity from time series

the methods to converge in the same estimation. This procedure completely lacks rationale because each estimator measures the signals' interactions in a different way and it captures different properties [START_REF] Bastos | A tutorial review of functional connectivity analysis methods and their interpretational pitfalls[END_REF].

The key solution to choose the connectivity estimator consists in a clear hypothesis according to the specific study. In practice, the scientific question underlying the research naturally determines the appropriate FC estimator. We consider as example the case of epileptic crises, which are episodes characterised by the increasing of brain activity causing convulsions or loss of awareness. A nonlinear estimator provides an appropriate description of this transient and rapid phenomenon. Since the epileptic crisis generally starts in the temporal lobe to propagate throughout the cortex, a directed method can be suitable to identify information sources [START_REF] Van Mierlo | Functional brain connectivity from eeg in epilepsy: Seizure prediction and epileptogenic focus localization[END_REF]. In this scenario, kernel Granger causality can be applied [START_REF] Marinazzo | Nonlinear connectivity by granger causality[END_REF]. Otherwise, if researcher are interested is the identification of synchronization mechanisms during the crisis, an undirected measure will be used, such as phase-locking value [START_REF] Aydore | A note on the phase locking value and its properties[END_REF], mutual information [START_REF] Vergara | An information theory framework for dynamic functional domain connectivity[END_REF], or wavelet coherence [START_REF] Chavez | Detecting dynamic spatial correlation patterns with generalized wavelet coherence and non-stationary surrogate data[END_REF].

Some other technical aspects guide the choice of the FC estimator [START_REF] He | Electrophysiological brain connectivity: theory and implementation[END_REF]. In fact, the study of rhythmic oscillatory phenomena requires a particular attention. Estimators defined in the frequency domain, such as the classical spectral coherence [START_REF] Carter | Coherence and time delay estimation[END_REF] and directed transfer function (DTF) [START_REF] Kaminski | A new method of the description of the information flow in the brain structures[END_REF], are suitable for the investigation of FC behaviours defined in frequency bands. Time signals can be transformed in the frequency domain through parametric (i.e. autoregressive models) or non parametric (i.e. Fourier or Hilbert transform). Since signals lie in the frequency domain, the frequency band of interest can be easily extracted. In the case of FC estimators defined in time-domain, like partial Granger causality [START_REF] Youssofzadeh | Partial granger causality analysis for brain connectivity based on event related potentials[END_REF] and transfer entropy [START_REF] Vergara | An information theory framework for dynamic functional domain connectivity[END_REF], band passing procedure can be applied.

Another important element to correctly choose the FC estimator is the available type of neuro-imaging data, because different techniques correspond to different temporal resolutions. EEG and MEG signals are characterized by high temporal resolution, in order of milliseconds, while fMRI typically have low time changes, in the order of seconds. Consequently, frequencies describing the EEG or MEG activities spans a broad interval, from portions of Hz to 100 Hz, according to the subject condition and the experiment [START_REF] Andreou | Theta and high-beta networks for feedback processing: a simultaneous eeg-fmri study in healthy male subjects[END_REF][START_REF] Bazanova | Interpreting eeg alpha activity[END_REF][START_REF] Yuan | Brain-computer interfaces using sensorimotor rhythms: current state and future perspectives[END_REF]. Conversely, fMRI data are characterized by a limited frequency range (less than 1 Hz) and FC estimators defined in time-domain are more appropriate [START_REF] Fox | Spontaneous fluctuations in brain activity observed with functional magnetic resonance imaging[END_REF].

In [START_REF] Gonzalez-Astudillo | Network-based brain computer interfaces: principles and applications[END_REF], we presented in detail state-of-the-art FC estimators with some associated critical aspects, such as the ability of some FC estimators to avoid the presence of spurious connections (e.g. imaginary coherence [START_REF] Nolte | Identifying true brain interaction from eeg data using the imaginary part of coherency[END_REF], partial directed coherence [START_REF] Sanei | Adaptive processing of brain signals[END_REF]) or to capture non-linear interactions (e.g. mutual information

Chapter 2. From brain signals to functional connectivity and graph analysis 7 [START_REF] Vergara | An information theory framework for dynamic functional domain connectivity[END_REF] and phase amplitude synchronization [START_REF] Aydore | A note on the phase locking value and its properties[END_REF]).

Graph representation of brain connectivity networks

To build a brain network, the first step consists in the identification of nodes and edges. The graph nodes are typically defined according to the specific neuroimaging technique used to extract the brain data. For voxel-based techniques, like fMRI and PET, nodes classically correspond to regions of interests (ROIs), identified using anatomical atlas [START_REF] Cohen | Defining functional areas in individual human brains using resting functional connectivity mri[END_REF][START_REF] Salvador | Neurophysiological architecture of functional magnetic resonance images of human brain[END_REF]. For sensor-based techniques, such as EEG and MEG, nodes usually correspond to recording sensors [START_REF] Lai | A comparison between scalp-and source-reconstructed eeg networks[END_REF]. In this latter case, brain nodes lies at the scalp level, but source-reconstruction techniques can be applied to re-define nodes on the cortex [START_REF] Edelman | Eeg source imaging enhances the decoding of complex right-hand motor imagery tasks[END_REF][START_REF] Jonmohamadi | Source-space ica for eeg source separation, localization, and time-course reconstruction[END_REF].

Neuroimaging techniques measure in specific modalities brain regional activities, without having access to their mutual relationships. Consequently, the edges' weights have to be inferred from signals through statistical procedures. Functional connectivity (FC) estimators, described in the previous section, are typically used to quantify links in a graph setting.

Considered together, nodes and links give rise to a new type of networked data, which cannot be analysed with standard techniques but requires appropriate tools coming from network science, a research domain merging graph theory, statistical mechanics and inferential modelling [START_REF] Alessandro | Twenty years of network science[END_REF]. Network science provides a novel perspective to analyse interacting data at several level [START_REF] Boccaletti | Complex networks: Structure and dynamics[END_REF][START_REF] Strogatz | Exploring complex networks[END_REF]. The networks of traffic movements, social relationships and molecular or protein interactions are only some examples [START_REF] Arias | Prediction of electric vehicle charging-power demand in realistic urban traffic networks[END_REF][START_REF] Hedayatifar | Us social fragmentation at multiple scales[END_REF][START_REF] Power | Collective ritual and social support networks in rural south india[END_REF][START_REF] Vinayagam | Controllability analysis of the directed human protein interaction network identifies disease genes and drug targets[END_REF] . In this scenario, a complex network is modelled as a graph [START_REF] Newman | Networks: an introduction[END_REF]. A graph G = {V, E, A} is defined as a set of nodes (or vertices) V with |V| = N, a set of links (or edges) E and an adjacency matrix A [START_REF] Newman | Networks: an introduction[END_REF]. The A i,j element differs from zero if there is a link e = (i, j) between the nodes i and j. At this stage, the graph is fully connected and weighted.

Widely used techniques in network domain involve thresholding strategies to reduce the number of links or to obtain a binary adjacency matrix. Reducing graph links provides advantages in improving the interpretabilty and mitigating false connections [START_REF] Vico Fallani | A topological criterion for filtering information in complex brain networks[END_REF][START_REF] Smith | Sift: Sphericaldeconvolution informed filtering of tractograms[END_REF]. The most intuitive way to proceed consists in defining a threshold on the number of edges or on the functional connectivity values.

This simple approach is parametric and the same analysis has to be often performed several times before obtaining stable results. Another possibility is to use nonparametric methods, based on the statistical or topological properties of the network [START_REF] Vico Fallani | A topological criterion for filtering information in complex brain networks[END_REF][START_REF] Tumminello | A tool for filtering information in complex systems[END_REF]. After thresholding, networks can be weighted or unweighted, depending on whether the links' weights are binarized. In Fig. 2.2, we report a

The Graph Signal Processing Framework

graphical representation of the main steps needed to estimate graph connectivity networks from time series.

Once the graph is completely defined through its nodes and edges, it is possible to extract some metrics, which concisely capture graph topological properties [START_REF] Gonzalez-Astudillo | Network-based brain computer interfaces: principles and applications[END_REF]. Network properties are defined at several scales. At a local scale, which refers to the node level, we typically have measures quantifying the node importance in the network. One example is the node strength which is computed as the weighted sum of the links including a node. At a meso-scale, which is denoted as the node group level, measures of nodes' tendency to group are defined.

One example is the study of modules and communities. At the global-level, we study properties related to the information transfer in the graph, such as globalor local-efficiency [START_REF] Fornito | Fundamentals of brain network analysis[END_REF]. From each node, brain signals are recorded. After choosing an appropriate FC estimator, functional connectivity is computed for each pair of nodes in order to obtain an adjacency matrix. Nodes and edges, corresponding to the elements in the adjacency matrix, constitute the final brain network.

The Graph Signal Processing Framework

There are many tools coming from different domains which are suitable to analyse graphs from different perspectives. This section presents graph signal processing (GSP) [START_REF] Sandryhaila | Discrete signal processing on graphs[END_REF] tools from a theoretical and practical point of view. The goal of GSP is to develop signal processing methods for data defined over graphs [START_REF] Ortega | Graph signal processing: Overview, challenges, and applications[END_REF]. Signals defined over graphs (SoGs) are a set of values lying on the graph nodes, linked by edges. Fig. 2.3 graphically shows SoGs with the underlying graph structure.

Graph signals have properties, as in signal processing, which have to be rede- fined in the specific scenario of graphs. For instance, one classical property is the smoothness, but the signal smoothness in this context needs to be defined. A smooth SoG can be a signal assuming similar values in close nodes. At the same time, smoothness can be referred to the frequency domain [START_REF] Ortega | Graph signal processing: Overview, challenges, and applications[END_REF].

Even basic and intuitive concepts in classical signal processing become challenging where signals are defined over a graph. One example is the translation operation, which does not have a direct meaning in the graph setting because the shift invariant property is not defined on irregular domains [START_REF] Sandryhaila | Discrete signal processing on graphs[END_REF]. Also the downsampling procedure, consisting on eliminating some samples in the classical signal processing, here is tricky because of the vertex dimension [START_REF] Ortega | Graph signal processing: Overview, challenges, and applications[END_REF][START_REF] Tanaka | Spectral domain sampling of graph signals[END_REF].

Another fundamental tool in signal processing is the Fourier Transform, which has been redefined as Graph Fourier transform (GFT) [START_REF] Shuman | The emerging field of signal processing on graphs: Extending high-dimensional data analysis to networks and other irregular domains[END_REF] in GSP. For the presentation of this transformation, a step back is needed. After defining the edges' weights obtaining the adjacency matrix A, we can compute the graph Laplacian, as follows:

L = D -A (2.1)
where D is the degree matrix, which is a diagonal matrix collecting in the i-diagonal element the sum of the edges including the node i.

We identify the set of orthonormal eigenvectors {u i } i=0,1,...,N-1 , where N is the number of graph vertices. {u i } relates to increasingly ordered eigenvalues 0

= λ 0 ≤ λ 1 ≤ λ 2 ... ≤ λ N-1 = λ max .
Laplacian eigenvectors {u i } play a fundamental role in GSP. In fact they are used as signals on graphs (s) and they represent a basis for the Graph Fourier Transform. The GFT for SoG s is computed as the projection of the SoG on the Laplacian eigenvectors:

ŝ(λ l ) = s H u l (2.2)

Representation of Time-Varying Networked Data

Graph Laplacian eigenvalues defined in the graph setting, have the similar notion of frequencies in the classical Fourier transform. Indeed, small eigenvalues correspond to eigenvectors which have small variation on connected vertices. Conversely, eigenvectors associated to large eigenvalues rapidly vary over connected nodes [START_REF] Ortega | Graph signal processing: Overview, challenges, and applications[END_REF][START_REF] Petrovic | Communityaware graph signal processing: Modularity defines new ways of processing graph signals[END_REF][START_REF] Shuman | The emerging field of signal processing on graphs: Extending high-dimensional data analysis to networks and other irregular domains[END_REF].

These considerations on graph eigenvectors will be used in Chapter 4 to define the denoising algorithm.

Among the other applications, GSP has been recently leveraged in brain network scenarios [START_REF] Huang | A graph signal processing perspective on functional brain imaging[END_REF]. Indeed, GSP can be applied to neuroimaging data and it provides tools to simultaneously investigate brain structure, contained in the graph, and the brain functioning, contained in the SoGs [START_REF] Mortaheb | A graph signal processing approach to study high density eeg signals in patients with disorders of consciousness[END_REF]. Medaglia et al. in [START_REF] Medaglia | Functional alignment with anatomical networks is associated with cognitive flexibility[END_REF] applied the GSP metrics of alignment and liberality to deepen the problem of attention switching with fMRI data. Findings show that subjects characterized by signals aligned with white matter structure can switch attention more rapidly. The same metrics (i.e. alignement and liberality, deriving from spectral graph filtering) have been applied by Bolton et al. in [START_REF] Bolton | Dynamics of brain activity captured by graph signal processing of neuroimaging data to predict human behaviour[END_REF] to predict behavioural variables, such as cognitive features, task skills and motor abilities. Specifically, authors studied the presence of temporal dynamics in the alignment and liberality. Results from the temporal investigation, have demonstrated the possibility of detecting network transitions through changes in GSP metrics.

Graph signal processing has recently been applied for brain-computer interface applications. Petrantonakis et al. [START_REF] Petrantonakis | Single-trial nirs data classification for brain-computer interfaces using graph signal processing[END_REF] have used GSP tools to extract novel features to improve the classification performances in NIRS-based BCIs systems.

Taken together, these first results show potential in the application of GSP to brain data analysis. In fact, GSP provides complete and flexible tools to analyse interacting systems, which find application in brain behaviour investigation at several levels and in a large set of applications, such as BCIs.

Representation of Time-Varying Networked Data

Nowadays networked data are studied by different perspectives and with several innovative tools. One intrinsic property characterizing almost all the data is the time-varying evolution. Social, biological, traffic networks change their descriptors in time and investigation tools must be able to capture this dynamic behaviour [START_REF] Kim | Inference of dynamic networks using time-course data[END_REF][START_REF] Preti | The dynamic functional connectome: State-of-the-art and perspectives[END_REF]. Concerning the brain investigation, many studies limit the analysis to static functional connectivity networks, assuming that the interactions between brain signals are stationary. Today we know the brain coupling mechanisms are timevarying processes [START_REF] Boveroux | Breakdown of within-and between-network resting state functional magnetic resonance imaging connectivity during propofol-induced loss of consciousness[END_REF][START_REF] Chang | Time-frequency dynamics of resting-state brain connectivity measured with fmri[END_REF][START_REF] Schrouff | Brain functional integration decreases during propofol-induced loss of consciousness[END_REF] with a dynamic alternation of synchronizations and desynchronizations that constitutes the basis for many neuro-physiological tasks, as perception and memory [START_REF] Freyer | Bistability and non-gaussian fluctuations in spontaneous cortical activity[END_REF][START_REF] Hutchison | Dynamic functional connectivity: promise, issues, and interpretations[END_REF].

In order to investigate the time-varying brain functioning, the simplest strategy is to reduce the length of the time interval, considering sliding windows with possible overlap. While reducing the time windows assures the signal (quasi) stationarity, the reliability of the FC estimates can rapidly decrease. This situation becomes even more critical when multivariate or non-linear FC estimators are adopted, since they need even longer time signals to compute reliable estimates [START_REF] Netoff | 11 detecting coupling in the presence of noise and nonlinearity[END_REF][START_REF] Pereda | Nonlinear multivariate analysis of neurophysiological signals[END_REF].

To overcome this limitation, a possible strategy for time-varying FC estimation consists in using techniques specifically designed for non-stationary signals, such as detrended fluctuation analysis [START_REF] Gö Güş | Identification of apnea-hypopnea index subgroups based on multifractal detrended fluctuation analysis and nasal cannula airflow signals identification of apnea-hypopnea index subgroups based on multifractal detrended fluctuation analysis and nasal cannula airflow signals[END_REF] and wavelet transformation [START_REF] Chavez | Detecting dynamic spatial correlation patterns with generalized wavelet coherence and non-stationary surrogate data[END_REF]. Among the other FC measures [START_REF] Gonzalez-Astudillo | Network-based brain computer interfaces: principles and applications[END_REF], some have been specifically developed for time-varying scenarios, such as wavelet coherence [START_REF] Tewarie | Tracking dynamic brain networks using high temporal resolution meg measures of functional connectivity[END_REF] and adaptive partial directed coherence [START_REF] Sanei | Adaptive processing of brain signals[END_REF].

This is a first-level analysis that deals only with short-time FC estimation procedures. Another possibility is to go further and model graph temporal evolution. One strategy is presented in [START_REF] Ozdemir | Recursive tensor subspace tracking for dynamic brain network analysis[END_REF], where factor decomposition [START_REF] Kolda | Tensor decompositions and applications[END_REF].

State-of-the-art graph learning methods use either simplified graph models for SoGs or models that are more expensive in terms of computational cost. In the case of time-varying graphs, these aspects are even exacerbated. Ghoroghchian et al. [START_REF] Ghoroghchian | Node-centric graph learning from data for brain state identification[END_REF] proposed a graph learning method based on representation learning on graphs. Specifically, the representation learning algorithm produces embeddings for graph vertices, by considering neighbouring nodes. This procedure was tested on intra-cranial electroencephalographic (iEEG) data recorded from epileptic patients.

In the context of graph learning, many recent studies address the problem of time-varying graph characterization by assuming the smoothness on the temporal domain [START_REF] Ghoroghchian | Node-centric graph learning from data for brain state identification[END_REF][START_REF] Natali | Online time-varying topology identification via prediction-correction algorithms[END_REF]. In this direction, Kalofolias et al. [START_REF] Kalofolias | Learning time varying graphs[END_REF] and Yamada et al. [START_REF] Yamada | Time-varying graph learning with constraints on graph temporal variation[END_REF] proposed advanced optimization algorithms. Specifically, in [START_REF] Yamada | Time-varying graph learning with constraints on graph temporal variation[END_REF], authors introduced a method to combine several regularizations for graph learning applications even when few observations are available. Jiang et al. [START_REF] Jiang | Dynamic graph learning based on graph laplacian[END_REF] developed a model to characterize graph temporal evolution in human connectome. Authors mathematically describe the problem as a quadratic objective function on graph vertices across short-time intervals. A regularization procedure is applied to reflect the smoothness and other properties, such as the graph Laplacian evolution.

Another strategy consists in modelling the spatio-temporal structure of a 3 -D graph. To this aim, Romero et al. [START_REF] Romero | Kernel-based reconstruction of space-time functions on dynamic graphs[END_REF] developed a kernel-based decomposition.

The advantage of this method consists in simultaneously estimating the graph structure on time and space domain, with a reduction of the number of vertices to be analysed. Liu et al. [START_REF] Liu | Graph learning based on spatiotemporal smoothness for time-varying graph signal[END_REF] proposed an alternative way to estimate 3 -D graph, by introducing a different algorithm based on smoothness prior to learn the graph and simultaneously identify the time-correlation pattern.

Ortiz et al. [START_REF] Ortiz-Jiménez | Sampling and reconstruction of signals on product graphs[END_REF] presented a different strategy based on the structure of product graphs, that, whenever obtained by applying Kronecker and Cartesian products, can be used for sampling and reconstruction applications. This method has broad applicability since product graphs can be used to model sensors' values acquired at different time points.

In [START_REF] Isufi | Observing and tracking bandlimited graph processes from sampled measurements[END_REF], Isufi et al. developed a sampling procedure for time-varying SoGs, to observe and track signals described by a linear state-space model. A mathematical study to identify the role of graph, graph signals and sample location is provided.

Chapter 2. From brain signals to functional connectivity and graph analysis 13

In conclusion, novel tools coming from different research domains pave the way for the study of brain organizational processes and time varying dynamics from an original perspective, constituted by merging signal processing, network theory and graph signal processing.

Chapter 3

Graph connectivity estimation to detect mental states

Functional connectivity (FC) is widely explored for its ability to capture functional interactions between brain signals. However, the role of FC in the context of brain-computer interface applications is still poorly understood.

To address this gap in knowledge, we consider a group of 20 healthy subjects performing a motor imagery (MI) task. We study two state-of-the-art FC estimators, i.e. spectral-and imaginary-coherence, and we describe how they change during MI tasks. We characterize the resulting FC networks by computing the node strength of each EEG sensor and we compare its discriminant ability with respect to standard univariate features. At the group level, we show that while spectral-coherence based features increase in the sensorimotor areas, those based on imaginary-coherence significantly decrease. We demonstrate that this opposite, but complementary, behaviour is

Introduction

Based on the detection of cognitive states from brain signals, brain-computer interfaces (BCIs) are explored for control, communication and rehabilitation, via the ability of subjects to voluntary modulate their brain activity through mental imagery [START_REF] Wolpaw | Brain-computer interfaces: Definitions and principles[END_REF][START_REF] Wolpaw | Brain-computer interfaces for communication and control[END_REF][START_REF] Wolpaw | The wadsworth center brain-computer interface (bci) research and development program[END_REF]. The ability to correctly identify the user's intent is therefore a key point to design BCI systems [START_REF] Allison | Could anyone use a bci? In Brain-computer interfaces[END_REF][START_REF] Blankertz | Predicting bci performance to study bci illiteracy[END_REF][START_REF] Carlson | Brain-controlled wheelchairs: a robotic architecture[END_REF][START_REF] Edlinger | How many people can use a bci system? In Clinical systems neuroscience[END_REF][START_REF] Thompson | Critiquing the concept of bci illiteracy[END_REF][START_REF] Vidaurre | Towards a cure for bci illiteracy[END_REF].

In this direction, investigators have explored different approaches based on various theoretical and experimental perspectives. One possibility is to look for the best mental strategy to detect the user response or to identify the adapted feedback to convey the most relevant information to the user [START_REF] Ahn | Performance variation in motor imagery braincomputer interface: a brief review[END_REF][START_REF] Müller-Putz | Towards noninvasive hybrid brain-computer interfaces: framework, practice, clinical application, and beyond[END_REF][START_REF] Vidaurre | Co-adaptive calibration to improve bci efficiency[END_REF].

Another possibility is to develop advanced signal processing methods and sophisticated classification algorithms to improve the signal-to-noise ratio and to correctly identify the user's intent [START_REF] Lotte | A review of classification algorithms for eeg-based brain-computer interfaces[END_REF][START_REF] Lotte | A review of classification algorithms for eeg-based brain-computer interfaces: a 10 year update[END_REF]. Although these methods provide performance increments, they are intrinsically blind to the neural mechanisms that enable to identify the user's cognitive state and may not have a direct physical or physiological interpretation [START_REF] Van Delden | Comparing unilateral and bilateral upper limb training: the ultra-stroke program design[END_REF]. However, this is crucial especially in clinical applications where brain functioning can be compromised and other solutions must be identified.

A different strategy consists in looking for alternative -potentially more informative -features reflecting the human brain organization processes. To this end, functional connectivity (FC) can be adopted to estimate the interaction between spatially distributed brain areas by quantifying the dependences between the regional activities [START_REF] Vico Fallani | Graph analysis of functional brain networks: practical issues in translational neuroscience[END_REF]. In contrast to univariate features such as frequency band power, FC appears more appropriate to grasp the oscillatory network processes involved in brain (re)organization during cognitive tasks [START_REF] Vico Fallani | Network neuroscience for optimizing brain-computer interfaces[END_REF]. Recent studies have demonstrated the potential of FC features in BCI, albeit the results are variable and difficult to compare because of the different FC estimators, tasks and limited number of subjects [START_REF] Brunner | Phase relationships between different subdural electrode recordings in man[END_REF][START_REF] Hamedi | Electroencephalographic motor imagery brain connectivity analysis for bci: a review[END_REF][START_REF] Song | Phase synchrony rate for the recognition of motor imagery in brain-computer interface[END_REF][START_REF] Wei | Amplitude and phase coupling measures for feature extraction in an eeg-based brain-computer interface[END_REF]. More importantly, the neurophysiological interpretation of FC features is still poorly understood in motor imagery tasks but this is crucial to eventually design alternative FC-based BCIs.

In order to investigate the actual possibility of applying FC features in BCI context, we consider two state-of-the-art FC estimators, i.e. the spectral coherence and imaginary coherence [START_REF] Carter | Coherence and time delay estimation[END_REF][START_REF] Nolte | Identifying true brain interaction from eeg data using the imaginary part of coherency[END_REF]. From a theoretical perspective, these estimators bring complementary information because the first quantifies the synchronization between the signal amplitudes while the latter also depends on their phase difference [START_REF] Nolte | Identifying true brain interaction from eeg data using the imaginary part of coherency[END_REF][START_REF] Rosenblum | Phase synchronization of chaotic oscillators[END_REF]. Our hypothesis is that the integration of these complementary features will allow a better characterization of the BCI-related cognitive states and that including them in the feature extraction block will improve the BCI accuracy as compared to standard approaches solely based on univariate features.

Material and methods

Notation

To verify our predictions, we explore brain FC networks extracted from EEG data recorded in a group of 20 healthy subjects performing the motor imagery (MI) of the right hand. In order to compare our approach with power spectrum features (P), we compute for each sensor the node strength, an intuitive graph metric which describes for each node its overall connectivity within the network. At the group level, we statistically compare the spatial patterns obtained from graph-related features associated to motor imagery and resting state. At the individual level, we evaluate the discriminant potential of network metrics by means of an off-line classification simulation.

To facilitate the reader, we list in In Table 3.1 the main notation used in this chapter.

Material and methods

Experimental protocol and preprocessing

Twenty healthy subjects (aged 27.60 The BCI experiment consisted in a standard 1D, two-target box task [START_REF] Wolpaw | The wadsworth center brain-computer interface (bci) research and development program[END_REF].

The subject was in front of a screen with a distance of 90 cm. When the target was up, the subject was instructed to imagine moving his/her the right hand (i.e.

grasping); when the target was down, the subject had to remain at rest. EEG A pre-preprocessing step preceded the analysis/ Specifically, we performed on the entire dataset an independent component analysis (ICA) to eliminate ocular and cardiac artefacts, via the Infomax algorithm [START_REF] Bell | An information-maximization approach to blind separation and blind deconvolution[END_REF] available in the Fieldtrip toolbox [START_REF] Oostenveld | Fieldtrip: open source software for advanced analysis of meg, eeg, and invasive electrophysiological data[END_REF]. The ICA was operated by the visual inspection of both time signals and their associated topographies. We removed no more than two independent components in average. In Appendix A we report one example of preprocessing on the first subject.

Functional connectivity and brain network features

We consider two state-of-the-art functional connectivity estimators [START_REF] Bastos | A tutorial review of functional connectivity analysis methods and their interpretational pitfalls[END_REF], i.e. spectral coherence (C) [START_REF] Carter | Coherence and time delay estimation[END_REF] and imaginary coherence (IC) [START_REF] Nolte | Identifying true brain interaction from eeg data using the imaginary part of coherency[END_REF].

While other FC estimators, directed and undirected, have been already applied in the BCI context [START_REF] Grosse-Wentrup | Understanding brain connectivity patterns during motor imagery for brain-computer interfacing[END_REF][START_REF] Gysels | Phase synchronization for the recognition of mental tasks in a brain-computer interface[END_REF][START_REF] Hamedi | Electroencephalographic motor imagery brain connectivity analysis for bci: a review[END_REF][START_REF] Krusienski | Value of amplitude, phase, and coherence features for a sensorimotor rhythm-based brain-computer interface[END_REF][START_REF] Pichiorri | Brain-computer interface boosts motor imagery practice during stroke recovery[END_REF], here we explore C and IC because of their relatively simplicity and intuitiveness.

Given two EEG time series y i and y j in a time interval T s , the computation of C ij and IC ij at the frequency ω k can be respectively obtained as:

C ij (ω k ) = P ij (ω k ) P i (ω k ) • P j (ω k ) 1/2 (3.1)

Material and methods

IC ij [ω k ] = P ij (ω k ) P i (ω k ) • P j (ω k ) 1/2 (3.2)
where P i (ω k ) contains the samples of the power spectral density P ii (e iω ); and P ij [ω k ] are samples of the cross-spectrum P ij (e iω k ) between y i and y j .

These quantities are evaluated by means of Welch's method with Hanning time windows of 1s and an overlap of 50% [START_REF] Welch | The use of fast fourier transform for the estimation of power spectra: a method based on time averaging over short, modified periodograms[END_REF]. While C has an intuitive interpretation because it captures linear correlations in the frequency domain, IC, by neglecting zero-lag contributions, is more robust to spurious connectivity due to volume condition [START_REF] Nolte | Identifying true brain interaction from eeg data using the imaginary part of coherency[END_REF].

As evident from Fig. 3.1, spectral coherence is generally higher and it more sensitive to short-distance interactions than imaginary coherence [START_REF] Cattai | Characterization of mental states through node connectivity between brain signals[END_REF].

∆ ij (ω k ) = φ i (ω k ) -φ j (ω k ) (3.3)
where φ i (ω k ) , φ j (ω k ) are the phase terms of the discrete Fourier transforms (DFTs) computed from y i and y j .

After estimating C, IC and ∆ for each pair of EEG channel, we obtain symmetric N × N matrices where N = 74 is the number of EEG electrodes. These matrices correspond to fully connected and weighted graphs of N nodes or edges, possibly studied via graph theoretic tools [START_REF] Vico Fallani | Graph analysis of functional brain networks: practical issues in translational neuroscience[END_REF]. In this thesis, we focus on an intuitive local centrality measure, i.e. the node strength S, which is computed as the sum of the weights of all links coming into each node. This metric describes in a simple and intuitive way how much one brain node, or EEG channel, is connected to all the others at a certain frequency ω k . thereby, node strengths are derived for each connectivity estimator used to extract the network, as follows :

S C i (ω k ) = N ∑ j=1 C ij (ω k ), (3.4)
S IC i (ω k ) = N ∑ j=1 IC ij (ω k ), (3.5)
S ∆ i (ω k ) = N ∑ j=1 ∆ ij (ω k ) (3.6)

Statistical Analysis and Classification

At group level, we average for each subject the associated connectivity matrices across trials and within frequency bands, namely: theta(θ) = 4 -7Hz, al pha(α) = 8 -13Hz, beta(β) = 14 -29Hz and gamma(γ) = 30 -40Hz. Node strength features are computed from each of these resulting networks. The same procedure is performed for power spectrum-based features.

We statistically compared connectivity and node strength values between the two mental states of MI and resting conditions. More specifically, for each condition we consider the distributions of the connectivity-based features obtained from the entire population of 20 subjects. We used non parametric permutation statistical tests (2000 permutations) with a statistical threshold of 0.05 [START_REF] Zar | Biostatistical analysis[END_REF] corrected for multiple comparisons with false discovery rate [START_REF] Benjamini | Controlling the false discovery rate: a practical and powerful approach to multiple testing[END_REF]. .

At individual level, we did not average the features across trials or within frequency bands. We let the classification procedure automatically select the set of best discriminant features for MI and resting for each subject. We only impose some constraints to limit the research complexity. First, we consider frequency bins from 4 to 40 Hz, because state-of-the-art results prove their involvement in similar motor tasks [START_REF] Neuper | Event-related dynamics of cortical rhythms: frequency-specific features and functional correlates[END_REF]. Second, we limit the research among a subset of electrodes spatially covering the sensorimotor areas [START_REF] Pfurtscheller | Spatiotemporal erd/ers patterns during voluntary movement and motor imagery[END_REF].

In order to investigate the contribution of the three different types of features (S C ,S IC and P) to the classification we considered all their possible combinations, i.e. seven in total. To normalize the values in each combination, we apply a z-

Results

score transformation to original features, i.e. channels × frequency bins. Then, we perform a 100 repeated ten-fold cross-validation classification with linear discriminant analysis (LDA) [START_REF] Lotte | A review of classification algorithms for eeg-based brain-computer interfaces[END_REF][START_REF] Lotte | A review of classification algorithms for eeg-based brain-computer interfaces: a 10 year update[END_REF]. In addition, we perform a sequential feature selection procedure [START_REF] Ferri | Comparative study of techniques for large-scale feature selection[END_REF] within a nested cross-validation framework to automatically identify the best discriminant features for each subject.

Results

EEG network connectivity changes during motor imagery

In order to verify that subjects were actually performing a proper motor imagery task, we show sensorimotor power decreases in Fig. 3.2 mainly in frequency bins within alpha and, more pronouncedly, beta band (p < 0.05, FDR-corrected). Network connectivity changes between MI and rest conditions are evident in all the frequency bands (Appendix B, Fig. 3.11, Fig. 3.12). Interestingly, in the beta band, results tend to occur in motor-related areas contralateral to the imagined movement. It is evident both for single connection and node strength values (Fig. 1A). These changes were even stronger when considering node strength values at individual frequency bins (p < 0.05, FDR-corrected) (Fig. 3.3B). Notably, the direction of the change is opposite depending on whether we estimate EEG networks with spectral-coherence (C) or with imaginary-coherence (IC). From our results, we report significant MI-related increases when we consider C estimators, while we observe significant decreases when using IC. In terms of spatial locations these differences involve both intra-hemispheric and inter-hemispheric interactions, while the largest changes in node strength tend to concentrate around the brain areas corresponding to the EEG electrode C3.

We observe a similar behaviour for gamma frequency band, although the most involved regions are more heterogeneously distributed. On the contrary, we cannot recognize similar trends for theta and alpha bands, where connectivity changes are in the same direction (Appendix B, Fig. 3.11 ,3.12).

We remark that the intensity of network changes is generally larger compared to power spectral P features (Fig. 3.3B, 3.2). For sake of clarity, we preliminary verified that results obtained with node strength are not correlated with those obtained by using P values (Pearson's correlation < 0.1). These results indicate that the motor imagery of the hand grasping produces significant brain network changes that might be useful to characterize MI-based

Results

BCI tasks. These changes uncover the existence of two parallel connectivity behaviours (i.e. increase for C and decrease for IC) that primarily involved sensorimotor areas within beta frequencies.

Modulation of amplitude and phase synchronization between brain signals

To deepen the nature of such dichotomy, we investigated more in detail the be- 3.2).

Table 3.2:

Spearman correlation coefficient between IC and ∆ for links including C3 for each subject.

Results

Taken together, our findings highlight a twofold mechanism which occurs during MI tasks. It consists in a simultaneous amplitude and phase synchronization among sensorimotor brain regions.

Mental state detection in single individuals

Finally, we test the ability of these brain connectivity features to discriminate MI and resting states at single subject level. To increase specificity, we consider a finer frequency resolution of 1 Hz -from 4 to 40 Hz -and we reduce the features to the EEG electrodes in the sensorimotor areas.

For each subject and each condition (i.e. MI and rest) we extract three type of features: power spectrum P, coherence-based node strength S C and imaginary coherence-based node strength S IC . We perform an automatic sequential forward feature selection within a cross-validation linear discriminant analysis (LDA). We consider the overall accuracy to quantify the average classification performance across cross-validation folds.

We report classification accuracies in Table 3.3. Results show that the best classification accuracy was in general moderate regardless of the feature combination.

Table 3.3: Average accuracy across cross-validation fold and repetitions is reported the best feature combination for each subject. In the top rows, we have classification results obtained when selected electrodes are located in both in contralateral and ipsilateral sensorimotor areas, ie. FC5, FC3, FC1, C5, C3, C1, CP5, CP3, CP1, FCZ, CZ, CPZ, FC2, FC4, FC6, C2, C4, C6, CP2, CP4, CP6. In the bottom rows, results come from classification framework including only controlateral electrodes, ie.FC5, FC3, FC1, C5, C3, C1, CP5, CP3, CP1.

After collecting accuracy values, for each subject we compute the relative dif-ference between the accuracy obtained with the best feature combination and the accuracy referred to P features. Let us remark that we consider the relative difference equal to zero when the best feature combination is P. We report results in Fig. 3.15. Interestingly, we observe that the inclusion of node strength features leads to improve performances in terms of relative difference with respect to P features alone. The improvement of classification performances when node strength features are included is particularly evident when features are selected from the controlateral hemisphere only (16 subjects over 20) as compared to when we consider both hemispheres (12/20). For features in controlateral areas, the performance improvement associated to the best features combination was 13% in average across

Results

subjects. To identify the spatial and spectral characteristics of the selected features, we show their cumulative occurrence in a frequency-sensor plot (Fig. 3.7).

In general, we observe a concentration of features in the 10 -14 Hz range within the C-CP zone. For both P and S C , the occurrences at higher frequencies tends to fade out (Fig. 3.7A,B), while the situation appears more heterogeneous for S IC features (Fig. 3.7C). For comparison's sake, we perform an alternative classification by substituting power spectrum features with those obtained with another state-of-the-art method, namely the filter bank common spatial pattern (FBCSP) [START_REF] Ang | Filter bank common spatial pattern algorithm on bci competition iv datasets 2a and 2b[END_REF], widely used in MI -based BCI applications. Results show that using FBCSP features does not significantly improve the overall accuracy and that node strength features still contribute to the relative performance improvement (Appendix D).

Taken together, these results show the potential of brain network features, de-rived from both amplitude and phase synchronization, to play an important role in the characterization of motor imagery states in healthy subjects.

Discussion

Brain activity changes during motor tasks have been largely documented by means of invasive and non-invasive neuroimaging techniques in non-human and human primates, as well as in animal models [START_REF] Murthy | Coherent 25-to 35-hz oscillations in the sensorimotor cortex of awake behaving monkeys[END_REF][START_REF] Pfurtscheller | Eeg event-related desynchronization (erd) and synchronization (ers)[END_REF][START_REF] Pfurtscheller | Event-related eeg/meg synchronization and desynchronization: basic principles[END_REF]. These changes are not limited to specific brain areas, but also occur in a synchronized manner across larger spatial scales reflecting the need for coordination of information exchanges to perform the task [START_REF] Crone | High-frequency gamma oscillations and human brain mapping with electrocorticography[END_REF][START_REF] Denker | Lfp beta amplitude is linked to mesoscopic spatio-temporal phase patterns[END_REF][START_REF] Jiang | Modulation of functional connectivity during the resting state and the motor task[END_REF][START_REF] Meirovitch | Alpha and beta band event-related desynchronization reflects kinematic regularities[END_REF].

Functional connectivity methods, measuring dependences between spatially remote brain areas, represent a unique opportunity to investigate large-scale brain network changes during motor imagery tasks from EEG recordings. Previous works reported FC changes in both healthy and diseased subjects [START_REF] Vico Fallani | Evaluation of the brain network organization from eeg signals: a preliminary evidence in stroke patient[END_REF][START_REF] Lynall | Functional connectivity and brain networks in schizophrenia[END_REF], as well as in BCI contexts [START_REF] Grosse-Wentrup | Understanding brain connectivity patterns during motor imagery for brain-computer interfacing[END_REF][START_REF] Gysels | Phase synchronization for the recognition of mental tasks in a brain-computer interface[END_REF][START_REF] Hamedi | Electroencephalographic motor imagery brain connectivity analysis for bci: a review[END_REF][START_REF] Krusienski | Value of amplitude, phase, and coherence features for a sensorimotor rhythm-based brain-computer interface[END_REF][START_REF] Pichiorri | Brain-computer interface boosts motor imagery practice during stroke recovery[END_REF]. However, different FC estimators have been applied in those studies and a deeper interpretation of obtained FC mechanisms was in general overlooked. As a result, a common direction and key interpretation of the FC processes observed during BCI motor-related tasks is still lacking [START_REF] Krusienski | Value of amplitude, phase, and coherence features for a sensorimotor rhythm-based brain-computer interface[END_REF].

To deepen this aspect, we study the intrinsic nature of two popular FC estimators, spectral coherence and imaginary coherence, extracted during motor imagery tasks in a group of healthy subjects. Our results indicate that motor imagery elicits two major parallel oscillatory phenomena in the beta frequency band: i) the increase of synchronization between the EEG signals' amplitudes, ii) a decrease of phase difference, which means an increase of synchronization between signals'

phases. Both amplitude and phase synchronization increases have been respectively reported in previous studies. The former typically codes for a basic substrate of neural communication [START_REF] Van Wijk | Neural synchrony within the motor system: what have we learned so far[END_REF], while the latter favours a further information binding [START_REF] Neuper | Event-related dynamics of brain oscillations[END_REF]. These FC modulations were region-specific and more localized in the sensorimotor areas. They particularly emerge at the node strength level and are not correlated with other regional measures, such as standard power spectral densities. We observe the presence of complementary mechanisms of amplitude and phase synchronization in the higher frequency bands (i.e. beta and gamma),

but not in theta and alpha bands. While this is in line with previous results reporting high-frequency FC changes in motor imagery-based BCI, the role the lowest frequency bands needs to be further explored [START_REF] Grosse-Wentrup | Understanding brain connectivity patterns during motor imagery for brain-computer interfacing[END_REF]. Future research should assess whether these network changes only reflect direct motor-related response or, also include indirect effects due to mirror-neuron activity as well as attentional efforts

Discussion

associated to the task complexity [START_REF] Murthy | Coherent 25-to 35-hz oscillations in the sensorimotor cortex of awake behaving monkeys[END_REF].

The ability to detect cognitive states from non-invasive neuroimaging recordings has concrete consequences in our daily-life, from the early detection of brain diseases to the design of brain-computer interface systems [START_REF] Müller | Machine learning for real-time single-trial eeg-analysis: from brain-computer interfacing to mental state monitoring[END_REF]. In the BCI context, much of the efforts has focused on the improvement of the classification algorithms, such as the recent developments in Riemannian geometry-based approaches [START_REF] Barachant | Riemannian geometry applied to bci classification[END_REF][START_REF] Gaur | A multi-class eeg-based bci classification using multivariate empirical mode decomposition based Bibliography filtering and riemannian geometry[END_REF]. Although these methods improve the overall classification performances, they generally lack of physical and physiological interpretations [START_REF] Vico Fallani | Network neuroscience for optimizing brain-computer interfaces[END_REF][START_REF] Lotte | A review of classification algorithms for eeg-based brain-computer interfaces[END_REF][START_REF] Lotte | A review of classification algorithms for eeg-based brain-computer interfaces: a 10 year update[END_REF]. Looking for alternative features, beyond the characterization of regional responses, is therefore a fruitful research field. [START_REF] Brunner | Online control of a brain-computer interface using phase synchronization[END_REF][START_REF] Daly | Brain computer interface control via functional connectivity dynamics[END_REF][START_REF] Hamedi | Motor imagery brain functional connectivity analysis via coherence[END_REF][START_REF] Krusienski | Value of amplitude, phase, and coherence features for a sensorimotor rhythm-based brain-computer interface[END_REF][START_REF] La Rocca | Human brain distinctiveness based on eeg spectral coherence connectivity[END_REF][START_REF] Wei | Amplitude and phase coupling measures for feature extraction in an eeg-based brain-computer interface[END_REF].

Whether the use FC and network-based features enable to significantly improve BCI classification improvements, is still under debate [START_REF] Caramia | Optimizing spatial filter pairs for eeg classification based on phase-synchronization[END_REF][START_REF] Daly | Brain computer interface control via functional connectivity dynamics[END_REF][START_REF] Gouy-Pailler | Topographical dynamics of brain connections for the design of asynchronous brain-computer interfaces[END_REF][START_REF] Gysels | Phase synchronization for the recognition of mental tasks in a brain-computer interface[END_REF][START_REF] Krusienski | Value of amplitude, phase, and coherence features for a sensorimotor rhythm-based brain-computer interface[END_REF][START_REF] Wei | Amplitude and phase coupling measures for feature extraction in an eeg-based brain-computer interface[END_REF].

Our results demonstrate that FC-based features bring complementary information with respect to power spectrum and have the potential to improve cognitive states' identification. This might have implications in the design of future BCIs, where the features selection is typically performed in the training phase. However, we remark the performance improvements measured here significantly vary across individuals and overall accuracies are generally moderated. Although our contribution highlights the potential of brain network features to address BCI inefficiency problem [START_REF] Thompson | Critiquing the concept of bci illiteracy[END_REF][START_REF] Vidaurre | Towards a cure for bci illiteracy[END_REF], further studies are crucial to assess the actual ability to discriminate between multiple mental states [START_REF] Salazar-Varas | Feature extraction for multi-class bci using eeg coherence[END_REF].

Methodological considerations

The estimation of spectral-and imaginary coherence assumes the stationarity of time series within the time window of computation. [START_REF] Nolte | Identifying true brain interaction from eeg data using the imaginary part of coherency[END_REF]. In our study, we considered time windows lasting 5 s, which could be too long the stationarity hypothesis [START_REF] Kaplan | Nonstationary nature of the brain activity as revealed by eeg/meg: methodological, practical and conceptual challenges[END_REF]. More in general, for real-time BCI applications the use of shorter time windows and FC estimators not requiring stationarity assumptions (e.g. wavelets [START_REF] Van Quyen | Comparison of hilbert transform and wavelet methods for the analysis of neuronal synchrony[END_REF], tracking algorithms [START_REF] Ozdemir | Recursive tensor subspace tracking for dynamic brain network analysis[END_REF]), would naturally allow to circumvent this issue.

We perform the study of FC networks on the EEG sensor space. Coherencebased FC estimates are affected by volume conduction distortions which cause spurious signal interactions [START_REF] Nolte | Identifying true brain interaction from eeg data using the imaginary part of coherency[END_REF][START_REF] Pascual-Marqui | Coherence and phase synchronization: generalization to pairs of multivariate time series, and removal of zero-lag contributions[END_REF]. Although source-reconstruction techniques can be used to mitigate such effect [START_REF] Jatoi | A survey of methods used for source localization using eeg signals[END_REF], we choose to work on the sensor space for two main reasons. The first motivation is a practical one. Indeed, we do not have access to the individual magnetic resonance images (MRIs) needed to have a realistic model of the head and its structures [START_REF] Baillet | Electromagnetic brain mapping[END_REF][START_REF] Edelman | Eeg source imaging enhances the decoding of complex right-hand motor imagery tasks[END_REF][START_REF] Michel | Eeg source imaging[END_REF]. The second motivation is more methodological because FC estimators can be really sensitive to signal transformations and results can strongly depend to the selected reconstruction Chapter 3. Graph connectivity estimation to detect mental states 29 algorithm [START_REF] Mahjoory | Consistency of eeg source localization and connectivity estimates[END_REF]. A detailed study on the effects of source-reconstruction is beyond the scope of our analysis. Future research is necessary to better investigate the stability of our results at source space level.

In order to combine intrinsically different estimators of brain activity (i.e.

power, node strengths), we perform a fusion at the feature level [START_REF] Ruta | An overview of classifier fusion methods[END_REF]. An alternative possibility is to perform the fusion at the classifier level, by combining the posterior probabilities of each separate classification [START_REF] Corsi | Integrating eeg and meg signals to improve motor imagery classification in brain-computer interface[END_REF]. The disadvantage of the latter approach is that it forces the research of significant features in each modality despite their absolute discriminant power. We prefer let the classifier to automatically select the best absolute combination of features.

Conclusions

Consistent with our hypothesis, we demonstrate the contribution of brain network connectivity features to detect cognitive states during typical MI-based BCI tasks.

More importantly, we find that hand MI is characterized by a twofold connectivity

Appendix C: Coherence and imaginary coherence estimation for sine waves

Let y i and y j be two signals of length T such that they only differ of a time shift t 0 .

In this simple situation, their cross-spectrum in the continuous domain reads as

P ij (e iω ) = P j (e iω )e iωt 0 (3.7)
where ω is the angular frequency of a signal [START_REF] Welch | The use of fast fourier transform for the estimation of power spectra: a method based on time averaging over short, modified periodograms[END_REF]. In the discrete domain, this can be rewritten as

P ij (ω k ) = P j (ω k )e i2πω k t 0 /T (3.8)
where the time shift t 0 becomes a linear phase term.

In the mathematical formulation of coherence

C C ij (ω k ) = P ij (ω k ) P j (ω k ) • P k (ω k ) 1/2 (3.9)
the numerator is the real part of the cross-spectrum and the exponential term in Eq. 2 is cancelled out. This indicates that C values do not depend on the amount of time shift between the signals.

Instead, in imaginary coherence IC,

IC ij (ω k ) = P ij (ω k ) P i (ω k ) • P j (ω k ) 1/2 (3.10)
it is trivial to show that there is a remaining term related to t 0 in the numerator.

Indeed, by rewriting the cross-spectrum via trigonometric functions:

P ij (ω k ) = P j (ω k )(cos (2πω k t 0 /T) + i sin (2πω k t 0 /T)) (3.11)
Hence, by taking the imaginary part one obtains

(P ij (ω k )) = P j (ω k ) sin (2πω k t 0 /T) (3.12)
This indicates that IC values do depend on the relative delay between the signals in a very specific way. More in general, in it has been emphasized that the estimated imaginary coherency between two time series can be expressed as a function of the instantaneous phase difference of their analytic signals [START_REF] Stam | Phase lag index: assessment of functional connectivity from multi channel eeg and meg with diminished bias from common sources[END_REF].

Appendix D: Alternative state-of-the-art classification based on FBCSP

In order to compare network-based feature with another state-of-the-art method, we perform a shrinkage LDA classification with FBCSP (filter bank common spatial pattern) features [START_REF] Ang | Filter bank common spatial pattern algorithm on bci competition iv datasets 2a and 2b[END_REF], alone and in combination with network features (S C and S IC ). The parameter µ f corresponds to a weight for each type of feature f . We observe from accuracy results that the classification with CSP did not perform better than classification with P, in our specific application.

Improving graph connectivity estimation with graph signal processing

Functional connectivity (FC) can be adopted to uncover brain (re)organization processes and it is applied to detect cognitive tasks, such as motor imagery (MI). The detection performances are strongly affected by the errors in FC estimation. In this chapter, we address the problem of robust graph FC estimation to improve the separability of cognitive states. Specifically, we propose a denoising algorithm based on graph signal processing (GSP) tools which acts on graph Laplacian.

In order to quantify the separation between different states, we derive a novel formulation of the Jensen divergence. We firstly apply the denoising procedure to synthetic data, demonstrating its potential in the improvement of the Jensen divergence between two simulated conditions. In order to address the problem of robust connectivity estimation together with the improvement of separability between mental states for BCI control, the combination of tools from different fields is fundamental. For instance graph signal processing is suitable for this scenario [START_REF] Ortega | Graph signal processing: Overview, challenges, and applications[END_REF][START_REF] Sandryhaila | Discrete signal processing on graphs[END_REF][START_REF] Shuman | The emerging field of signal processing on graphs: Extending high-dimensional data analysis to networks and other irregular domains[END_REF]. It has already been applied in context of biological data, and specifically with brain data [START_REF] Huang | A graph signal processing perspective on functional brain imaging[END_REF][START_REF] Medaglia | Functional alignment with anatomical networks is associated with cognitive flexibility[END_REF].

GSP has the potential to simultaneously integrate properties of the brain structure, represented by the graph itself, and brain functioning, represented by signals on graph.

Other fundamental tools to deal with the functional connectivity estimation problem come from signal detection theory. Detection techniques can be applied to measure statistical differences between FC features associated to two brain states, corresponding to motor imagery and resting state in our application. There are several measures that can be adopted in this context, such as the Likelihood Ratio (LR) of the features [START_REF] Kullback | Information theory and statistics[END_REF][START_REF] Pezeshki | Canonical coordinates are the right coordinates for low-rank gauss-gauss detection and estimation[END_REF][START_REF] Poor | An introduction to signal detection and estimation[END_REF][START_REF] Scharf | Low rank detectors for gaussian random vectors[END_REF][START_REF] Scharf | Statistical signal processing[END_REF][START_REF] Schürmann | Pattern classification: a unified view of statistical and neural approaches[END_REF][START_REF] Theodoridis | Pattern recognition[END_REF][START_REF] Trees | Detection, Estimation, and Modulation Theory-Part l-Detection, Estimation, and Linear Modulation Theory[END_REF] and the linear detector which maximizes the so-called deflection [START_REF] Chevalier | Complex linear-quadratic systems for detection and array processing[END_REF][START_REF] Picinbono | Detection and contrast[END_REF][START_REF] Picinbono | On deflection as a performance criterion in detection[END_REF][START_REF] Picinbono | Optimal linear-quadratic systems for detection and estimation[END_REF].

Obtaining a measure of the distance between two states is simplified under the assumption of normally distributed observations. Indeed, the maximum deflection test coincides with the LR test for normally distributed observations with equal conditional variance and different conditional means. Under this hypoth-

Related Work

esis, LR test can be extended to a linear quadratic detector so as to deal with observations with different conditional variances [START_REF] Picinbono | Detection and contrast[END_REF][START_REF] Picinbono | Optimal linear-quadratic systems for detection and estimation[END_REF]. To evaluate a measure of separability between features under the two conditions, the Jensen divergence can be explored, corresponding the maximum deflection test performance [START_REF] Pezeshki | Canonical coordinates are the right coordinates for low-rank gauss-gauss detection and estimation[END_REF][START_REF] Scharf | Low rank detectors for gaussian random vectors[END_REF][START_REF] Schürmann | Pattern classification: a unified view of statistical and neural approaches[END_REF].

In this chapter we propose a novel graph Laplacian denoising algorithm, able to improve the brain connectivity estimates. We propose a subspace-based Laplacian denoising which simultaneously preserves graph connectivity structure and rejects noise components. This algorithm preserves i) the subspaces more directly associated to the graph topology, summarized by the eigenvectors corresponding to the smallest Laplacian eigenvalues, ii) the subspaces associated to favourable signal-to-noise ratio, summarized by the eigenvectors corresponding to the largest Laplacian eigenvalues. We demonstrate that this twofold approach enhances the separability between two states.

In order to quantify the separability between two conditions and to measure the improvement achieved by the Laplacian denoising in the FC estimates, we present an original formulation of the Jensen divergence. The improvement of the J-divergence of the graph Laplacian coefficients under different states is proved by numerical simulations on synthetic data.

Finally, we test our method on real EEG data recorded during motor imagerybased BCI experiments, and we prove that our denoising algorithm increases the J-divergence evaluated between motor imagery and resting states, even when graph Laplacians are computed from short time-windows. As result of the Jdivergence analysis, we attribute a score to each Laplacian coefficient representing its marginal contribution to the total J-divergence. The score admits a relevant biological interpretation confirming the efficacy of the approach. These results can be assessed by further studies on the brain connectivity features.

To facilitate the reader, we list in Table 4.1 the main notation used in the chapter.

Related Work

The problem of graph connectivity estimation has been widely explored in literature in several fields, from neuroscience to signal processing and graph theory [START_REF] Boccaletti | Complex networks: Structure and dynamics[END_REF][START_REF] Bullmore | Complex brain networks: graph theoretical analysis of structural and functional systems[END_REF][START_REF] Segarra | Statistical graph signal processing: Stationarity and spectral estimation[END_REF]. Many state-of-the-art graph learning approach present oversimplified models to signals defined over graphs to overcome computational and memory issues. Some recent works propose different methods to address graph learning problems. For example, in [START_REF] Ghoroghchian | Node-centric graph learning from data for brain state identification[END_REF] [START_REF] Ortega | Graph signal processing: Overview, challenges, and applications[END_REF][START_REF] Sandryhaila | Discrete signal processing on graphs[END_REF][START_REF] Shuman | The emerging field of signal processing on graphs: Extending high-dimensional data analysis to networks and other irregular domains[END_REF]. GSP has already demonstrated the potential to characterize brain functioning in [START_REF] Huang | Graph frequency analysis of brain signals[END_REF] and [START_REF] Huang | A graph signal processing perspective on functional brain imaging[END_REF]. Indeed, GSP representation naturally reflects the human brain, where the structure is described by the graph itself while brain functioning corresponds to graph signals. An interesting application is represented by graph filtering [START_REF] Ortega | Graph signal processing: Overview, challenges, and applications[END_REF][START_REF] Rabbat | Graph laplacian distributed particle filtering[END_REF], that is useful to extract brain properties [START_REF] Rui | Dimensionality reduction of brain imaging data using graph signal processing[END_REF]. In [START_REF] Wang | Topographic regularity for tract filtering in brain connectivity[END_REF], authors present a mathematical method to track brain fibers in order to describe neurophysiological processes. The model, based on GSP techniques, selects a subset of graph eigenvectors representing a basis for filtering fiber tracts from brain imaging data.

GSP has already been applied in the context of brain-computer interfaces with

NIRS signals [START_REF] Petrantonakis | Single-trial nirs data classification for brain-computer interfaces using graph signal processing[END_REF]. Specifically, GSP analysis is applied in [START_REF] Petrantonakis | Single-trial nirs data classification for brain-computer interfaces using graph signal processing[END_REF] for feature extraction to obtain spatial information from the NIRS signals and it is proved to increase classification performances.

Classical signal processing and eigenvector-based filtering have already been 4.3. Signal Model adopted with brain data [START_REF] Spencer | Adaptive filters for monitoring localized brain activity from surface potential time series[END_REF][START_REF] Strobach | Event-synchronous cancellation of the heart interference in biomedical signals[END_REF]. In [START_REF] Chen | Eigenvector based spatial filtering of fetal biomagnetic signals[END_REF] and [START_REF] Zhang | Eigenvector-based spatial filtering for reduction of physiological interference in diffuse optical imaging[END_REF], eigenvector-based filtering has been applied to fetal magnetic signals and diffuse optic imaging data to have more localized activities and reduce artefacts and noise. In particular, in [START_REF] Zhang | Eigenvector-based spatial filtering for reduction of physiological interference in diffuse optical imaging[END_REF], classical eigenvector-filtering, which is typically based on eigenvectors associated to larger eigenvalues, is used in diffuse optical imaging in order to improve connectivity estimation.

In the following sections, we propose a novel Laplacian denoising algorithm, and we demonstrate that it improves the detectability of two states. To this aim,

we provide an analysis of the J-divergence, which naturally quantifies the distance between two distributions. Recently, the J-divergence has been applied in [START_REF] Li | Time series irreversibility analysis using jensen-shannon divergence calculated by permutation pattern[END_REF] to investigate the time series' irreversibility . Another recent application of the J-divergence is described in [START_REF] Nielsen | On a generalization of the jensen-shannon divergence[END_REF], where authors propose a novel framework to vector-skew the J-divergence. This method allows to preserve J-divergence properties and simultaneously to fine tune parameters for specific applications.

J-divergence has been also adopted in BCI design, to address one of the most challenging problem of EEG-based BCIs, that is the long calibration time. Indeed, the number of data necessary to calibrate the model is generally high, because of the presence of noise and artifacts the non-stationarity of EEG data. In [START_REF] Giles | A subject-to-subject transfer learning framework based on jensen-shannon divergence for improving brain-computer interface[END_REF],

authors propose a subject-to-subject transfer learning to improve the classification performance when limited training data are available. J-divergence is used in a transfer learning framework to test the method by comparing the data of the target subject with the data from previous subjects.

In the following sections, we study the J-divergence under a different points of view, specifically i) we assess the performance of the denoising algorithm in distinguishing brain states and ii) we provide a score definition for Laplacian coefficients based on their contribution to the total J-divergence.

Signal Model

We analyse signals defined on an undirected, connected, weighted graph G = {V, E, A}, which consists in a finite set of vertices (or nodes) V with |V| = N, a set of edges (or links) E and a weighted adjacency matrix A. If there is an link e = (i, j) connecting nodes i and j, the element A i,j represents the weight of the link; otherwise, A i,j = 0. The graph Laplacian, is a real symmetric matrix defined * as:

L = D -A (4.1)
* We refer here to the non-normalized graph Laplacian, also called the combinatorial Laplacian.

Chapter 4. Improving graph connectivity estimation with graph signal processing 45

where D is the degree matrix, which is a diagonal matrix whose i th diagonal element d i is equal to the sum of the weights of all the edges incident to vertex i.

We define the set of orthonormal eigenvectors {u i } i=0,1,...,N-1 , corresponding to increasingly ordered eigenvalues 0

= λ 0 ≤ λ 1 ≤ λ 2 ... ≤ λ N-1 = λ max .
In GSP, the Laplacian eigenvectors are represented as SoGs and constitute a basis for the Graph Fourier Transform (GST). The GFT is defined for a SoG s as the projection of s on the l-th eigenvector of the graph Laplacian:

ŝ(λ l ) = s H u l (4.2)
The graph Laplacian eigenvalues λ l , l = 0, • • • N -1 have an similar interpretation to Fourier transform frequencies. Indeed, eigenvectors associated to smaller eigenvalues present smoother variations over connected nodes.

In many applications, such as brain functional connectivity estimation, SoGs are represented by discrete sequences, obtained by sampling time signals acquired at each graph node.

We denote the discrete sequences collected during an observation period T oss with sampling pace T s as y n [νT s ], n = 0, . . . N -1, k = 0, . . . N s , N s = T oss /T s , or in vector form as y[νT s ] = [y 0 [νT s ] . . . y N-1 [νT s ]] The vector sequence y[νT s ], k = 0, . . . N s is used to compute the adjacency matrix A by estimating a similarity metric on each pair of nodes. There are many state-of-the-art procedures to estimate A i,j , i, j = 0, • • • N-1, which quantify links' weights according to specific interaction property [START_REF] Bastos | A tutorial review of functional connectivity analysis methods and their interpretational pitfalls[END_REF][START_REF] Carter | Coherence and time delay estimation[END_REF][START_REF] Friston | Functional and effective connectivity: a review[END_REF][START_REF] Nolte | Identifying true brain interaction from eeg data using the imaginary part of coherency[END_REF]]. Since we do not have access to the real adjacency matrix A, its estimated version  is actually used. It contains the connectivity values Âi,j estimated for each graph node pair (i, j), i, j = 0, • • • N -1.The estimated degree matrix D is derived, so as to compute the estimated Laplacian L as in Eq. Le us remark that any estimation error on the adjacency matrix affects the Laplacian estimate, and it results into less distinguishable connectivity states. In the following section we propose a denoising algorithm of the graph Laplacian in order to improve the separation of connectivity states. We generally refer to functional connectivity pattern associated to one condition as connectivity state.

Graph Connectivity Denoising

Graph Connectivity Denoising

In order to describe the graph Laplacian denoising algorithm, we consider the eigenvalue decomposition of the estimated Laplacian L as follows:

L = N-1 ∑ i=0 λi ûi ûH i (4.4)
Perturbations affect graph Laplacian estimation in terms of eigenvalues and/or eigenvectors. To elaborate on the effect of perturbations, we explicit the first, second and third order error contributions to L as:

L = N-1 ∑ i=0 (λ i + λ i )(u i + u i )(u i + u i ) H = N-1 ∑ i=0 λ i u i u H i L + λ i u i H u i + λ i u i u H i + λ i u i u H i first order error + λ i u i H u i + λ i u i H u i + λ i u i u H i second order error + λ i u i H u i third order error (4.5)
We can approximate the estimated graph Laplacian at the first order as the sum of N terms:

L ≈ N-1 ∑ i=0 (λ i + λ i )u i u H i + λ i (u i H u i + u i u H i ) (4.6)
Eq.(4.6) shows that the first order error is due to relative perturbation of the Laplacian eigenvalues as well as of the eigenvectors' direction. We are interested in the Laplacian components whose perturbation is low because either the relative eigenvalue perturbation λ i /λ i or the eigenvector perturbation u i is small. We consider the set of orthonormal eigenvectors U ALL = ûl , l = 0, 1, ..., N -1 with increasingly eigenvalues 0 = λ0 ≤ λ1 ≤ λ2 ... ≤ λN-1 := λmax , and we consider three subsets of eigenvalues and associated eigenvectors: 1) the subset U L with the N L smallest eigenvalues; 2) the subset U H with the N H largest eigenvalues; and 3) the subset U M with the remaining

N M = N -N L -N H central eigenvalues, with U L ∪ U M ∪ U H = U ALL .
The first consideration is that the N H largest eigenvalues are more robust to eigenvalue perturbations; this hypothesis is widely adopted in classical signal processing, where the subspace U H is used for the estimation of the covariance matrix because of its favourable signal-to-noise ratio [START_REF] Ollila | Complex elliptically symmetric distributions: Survey, new results and applications[END_REF].

Chapter 4. Improving graph connectivity estimation with graph signal processing 47

The second consideration is based on recent results [START_REF] Ceci | Graph signal processing in the presence of topology uncertainties[END_REF], according to which the subspace U L is partially robust in terms of eigenvector perturbations. In fact, Ceci and Barbarossa [START_REF] Ceci | Graph signal processing in the presence of topology uncertainties[END_REF] demonstrated that a connectivity estimation error on the A ij adjacency matrix element, i.e. on the weight of the link between the ith and the j-th nodes, causes a perturbation u i of the i-th eigenvector related to the difference between the i-th and the j-th coefficients of u i . As consequence, eigenvectors smoothly varying across the i-th and the j-th nodes are less affected by errors on A ij .

On the other hand, it is well known in GSP theory that U L eigenvectors corresponds to low frequency elements in the Graph Fourier Transform [START_REF] Luxburg | A tutorial on spectral clustering[END_REF], [START_REF] Shuman | The emerging field of signal processing on graphs: Extending high-dimensional data analysis to networks and other irregular domains[END_REF] since they are characterized by the smallest changes over connected graph nodes.

Thereby, the eigenvectors in U L are characterized by a natural resilience to connectivity estimation error between connected nodes. To sum up, the eigenvectors in U L are strongly related to the network topology, and therefore they need to be involved in the proposed denoising method.

Stemming on these observations, we introduce a novel denoising method which preserves the contribution to the Laplacian associated to the subspaces U L , U H while eliminating those associated to the subspace U M . In formulas, given the estimated Laplacian

L = ∑ i∈U L ∪U M ∪U H λi ûi ûH i (4.7)
we write the denoised Laplacian L as follows:

L = ∑ i∈U L λi ûi ûH i + ∑ i∈U H λi ûi ûH i (4.8)
The proposed graph Laplacian denoising algorithm presents a twofold approach which simultaneously preserves

• the subspace U L , which is directly related to the graph functional connectivity;

• the subspace U H , which is estimated with a favourable signal-to-noise ratio.

The proposed Laplacian denoising method, is synthetically presented in Algorithm 1. In essence, our method preserves the information relevant for graph connectivity characterization, and it rejects noisy components. In order to measure the improvement achieved in terms of connectivity states' separability, we use to the J-divergence as a measure of the distance between two states. In the following section we derive a novel formulation of the J-divergence.

Jensen divergence of connectivity states

Many state-of-the-art measures can be adopted to quantify the separability of two connectivity states [START_REF] Basseville | Divergence measures for statistical data processing-an annotated bibliography[END_REF], as represented by the Laplacian matrices L. Here, we apply the J-divergence to describe the separability of two connectivity states, and we propose a novel formulation to adapt it to our problem. Thus, J-divergence is applied to identify the Laplacian coefficients that are mostly relevant for detection and to measure the improvement achieved by the denoising algorithm.

We will assume in this analysis that the Laplacian coefficients obtained after the denoising algorithm are normally distributed. We remark that the Gaussian assumption stands in many applications * , including the case of FC estimation on real brain signals, and thereby it is often assumed in the literature, e.g. for Laplacian computation procedures [START_REF] Ortega | Graph signal processing: Overview, challenges, and applications[END_REF]. In our case, we assume that the vector l = Vec( L) is distributed according to a multidimensional Gaussian probability whose mean vector and covariance matrix differ under two connectivity states, referred to as the null and the alternative hypotheses H 0 , H 1 in the following: †

H 0 : l ∼ N (η 0 , K 0 ) H 1 : l ∼ N (η 1 , K 1 ) (4.9)
We compute the Jensen divergence in order to obtain an information theoretic measure of distance between l under H 0 and H 1 .The J-divergence is defined as the expected value of the difference of the Log Likelihood Ratio under the two hypothesis H 0 and H 1 [START_REF] Pezeshki | Canonical coordinates are the right coordinates for low-rank gauss-gauss detection and estimation[END_REF]. The J-divergence formulation enables us to evaluate to which extent the connectivity states represented by the Laplacian elements are distinguishable from each other.

Let us first assume that the moments of the Laplacians η 0 ,η 1 ,K 0 ,K 1 are known.

Detection procedure can be applied on a linear transformation of original observations:

x = T l-η 0
where η def = T (η 1 -η 0 ) and T = T (K 0 , K 1 ) is an affine transformation that simultaneously ‡ whitens the observations in the H 0 hypothesis and produces uncorre-* The reason why this occurs is that the Gaussian assumption tightly models laplacian diagonal elements, computed in each row as the sum of extradiagonal elements in that column, as well as extradiagonal elements which are often computed as the result of correlation estimates.

† The notation l ∼ N η j , K j , with j ∈ {0, 1} indicates that the random vector l is Gaussian distributed with mean vector η j and covariance matrix K j . ‡ The matrix T and the diagonal matrix Σ 2 def = diag(σ In real systems, the Laplacian moments η 0 , η 1 , K 0 , K 1 can be either estimated in a training set, e.g. during a BCI training, or tracking procedures can be applied [START_REF] Chang | Semiparametric estimation and selection for nonstationary spatial covariance functions[END_REF][START_REF] Greco | Maximum likelihood covariance matrix estimation for complex elliptically symmetric distributions under mismatched conditions[END_REF][START_REF] Guttorp | Covariance structure of spatial and spatiotemporal processes[END_REF]. Besides, data in the transformed domain x can be obtained even avoiding computation of moments and of T, by applying the Laplacian coefficients l to a trained network [START_REF] Ali | Biometricnet: deep unconstrained face verification through learning of metrics regularized onto gaussian distributions[END_REF], in order to enforce the afore-mentioned statistical constraints.

The observation model becomes:

H 0 : x ∼ N (0, I) versus H 1 : x ∼ N η, Σ 2 (4.10)
The J-divergence is defined as:

J def = E (R(x)|H 1 ) -E (R(x)|H 0 ) (4.11)
eralized eigenvectors and the generalized eigenvalues matrices of the pencil (K 1 , K 0 ), respectively. Given any square root

Q 0 of K -1 0 , i.e. such that Q H 0 • K 0 • Q 0 = I, we may conveniently employ the unitary transformation V 1 obtained from the eigenanalysis Q H 0 • K 1 • Q 0 = V 1 • Λ 1 • V H 1 ; in fact, it is easily proved that the matrix T = V H 1 • Q H 0 verifies T • K 0 • T H = I ; T • K 1 • T H = Σ 2 with Λ 1 = Σ 2 .
R(x) = x H I -Σ -2 x + 2η H Σ -2 x (4.12)
Here, we associate the variables x n whose variance σ 2 n = 1 to the first P indexes and the remaining ones to the indexes n = P + 1, . . . , N † so as to rewrite the LLR as follows: The P variates x n , n = 1, . . . , P, have different conditional variances under the hypotheses H 0 , H 1 and they contribute to the LLR by the P terms quadratic terms.

R(x) = N ∑ n=1 1 σ 2 n σ 2 n -1 x 2 n + 2η n • x n = P ∑ n=1 σ 2 n -1 |x n | 2 + 2η n • x n σ 2 n + N ∑ n=P+1 2η * n • x n
The N -P variates x n , n = P + 1, . . . , N , with equal unitary conditional variances under the hypotheses H 0 , H 1 , contribute to the LLR by the N -P linear terms.

To deepen the J-divergence, we present the following theorem, whose demonstration is in Appendix E.

Theorem 1 Let ξ be a vector formed by the N statistically independent random variables:

ξ n = x n + η n σ 2 n -1 2 , n = 1, . . . , P ξ n = x n , n = P + 1, . . . , N (4.15)
* The Log-Likelihood Ratio R(x) is widely adopted classical detection problems: R(x)

H 1 ↑ ≷ ↓ H 0 θ,
where θ is selected according to the desired detection versus missing probability trade-off. † We might have P = N or P = 0.

Chapter 4. Improving graph connectivity estimation with graph signal processing 51

The LLR is written as R(x) = a H LLR • ξ being a LLR constant coefficients defined as in Eq.(4.26) and the J-divergence in Eq.(4.11) is computed as follows:

J = P ∑ n=1 σ n -σ -1 n 2   1 + |η 2 n | σ n σ n + σ -1 n σ n -σ -1 n 2   + N ∑ n=P+1 2|η 2 n | = P ∑ n=1 J (σ,η) n + N ∑ n=P+1 J (η) n (4.16)
Theorem 1 generalizes the result in [START_REF] Pezeshki | Canonical coordinates are the right coordinates for low-rank gauss-gauss detection and estimation[END_REF][START_REF] Scharf | Statistical signal processing[END_REF] where only the case of variables having equal conditional means and different covariances (i.e. η 1 = η 0 , K 1 = K 0 ) has been considered.

The J-divergence as formulated in Eq.(4. 16) is a measure of the statistical distance of the Laplacian coefficients under two conditions, and it is suitable to quantify the separability of connectivity states achieved by the denoising algorithm described in section 4.4.

In addition, the J-divergence analysis provides an information on the variables that mostly contribute to the states separability. From Eq.(4. [START_REF] Bastos | A tutorial review of functional connectivity analysis methods and their interpretational pitfalls[END_REF], we see a one-toone correspondence between the transformed variables x n and the terms of the J-divergence J; besides, this term can be:

J (σ,η) n = σ n -σ -1 n 2   1 + |η 2 n | σ n σ n + σ -1 n σ n -σ -1 n 2    J (η) n = 2|η 2 n | (4.17)
The nature of J-divergence terms depends on whether the variable changes both in conditional mean and standard deviation, or in conditional mean only. The functions

J (σ,η) n , J (η)
n are plotted in Fig. 4.2 for η between 0 and 1 and σ -1 from 0 to 10. Fig. 4.2 shows that a conditional variance change gives a higher contribution to J than an equal conditional mean change.

In conclusion, the J-divergence analysis highlights the contribution of each transformed variable to the separability of the connectivity states. Consequently, it suggests a ranking procedure for x n according to their contribute to J n , in order to identify variables which mostly differ under the two hypothesis. This opens new ways to score Laplacian coefficients.

J-Divergence based Laplacian coefficients scoring

After the J-divergence analysis, we are able to identify Laplacian elements (i.e. links weights or nodes degrees), mostly contributing to the connectivity states separability. This is obtained by associating a score to each Laplacian element according to its contribution to the J-divergence.

By definition, the Laplacian elements ln , n = 0, • • • N -1, are combined in the transformed domain and they generate the variable x n . We define a score S n computed by means of a backpropagation procedure of the J n terms on each contributing Laplacian coefficient. Specifically, the n-th coefficient score is computed as

S n = ∑ n J n • t nn ∑ k t nk (4.18)
where we recognize that the weight t nn that brings the contribution of the n-th

Laplacian coefficient to the n-th transformed variable is normalized with respect to the sum ∑ k t nk of the weights of all the contributing coefficients.

A representation of the score rational is provided in Fig. To sum up, the score procedure enables to numerically quantify the relevance of the Laplacian elements ln to separate connectivity states. The Algorithm 2 presents the main steps of the J-divergence computation and scoring procedure.

Results on synthetic data

In this section, we investigate the ability of graph Laplacian denoising described in section 4.4 to increase the J-divergence between two connectivity states derived from synthetic SoGs. We firstly define a graph and a model for SoGs under two connectivity states, identified to simulate an over-simplified model of brain EEG signals functional connectivity.

For comparison's sake, we also consider the case of Laplacian without filtering U ALL and and other eigenvector-based filters ( i.e. U L , U H ). Here, we explain in detail the Laplacian generation procedure and the analyses we perform.

Signal on Graph generation and connectivity estimation

In order to test our procedure on synthetic data, we define signals at graph nodes under two hypothesis H 0 and H 1 to obtain graphs related to two connectivity

= N L -1 ∑ i=0 λi ûi ûH i + N-1 ∑ i=N-N H λi ûi ûH i states.
Under H 1 , we model the brain activity by considering H generator signals

s (h) [νT s ], h = 0, • • • H -1.
Each generator signal contributes to the signals measured over a subset G (h) , h = 0, • • • H -1 of nodes identified by the non-zero components of the N × 1 binary vector g (h) , h = 0, • • • H -1. Perturbations are also present on the model, both in term of a noise component w[νT s ] as well as a common component across all the nodes b[νT s ] • 1. Under H 0 , only these perturbation-related components are observed. Taken together, we define the vector of the observed signals y[νT s ] under the two hypotheses H 1 and H 0 as follows:

H 1 : y[νT s ] = H-1 ∑ h=0 s (h) [νT s ] • g (h) + w[νT s ] + b[kT s ] • 1 H 0 : y[νT s ] = w[νT s ] + b[νT s ] • 1 (4.19)
The noise w[νT s ] is a realization of a discrete, stationary, white Gaussian process, with

E{w[ν]} = 0, E{w[ν]w[ν] T } = σ 2
w I ∀ν; the samples of discrete sequences b[νT s ] are extracted from a zero mean Gaussian random distribution with variance σ 2 b ; and s (h) [νT s ], h = 0, • • • H -1 are extracted from a zero mean unit variance Gaussian random variable.

Once SoG samples y[νT s ] are defined, we estimate the adjacency matrix. There are many state-of-the-art methods to estimate the adjacency matrix such as spectral coherence [START_REF] Carter | Coherence and time delay estimation[END_REF], wavelet coherence [START_REF] Chavez | Detecting dynamic spatial correlation patterns with generalized wavelet coherence and non-stationary surrogate data[END_REF], Granger causality [START_REF] Youssofzadeh | Partial granger causality analysis for brain connectivity based on event related potentials[END_REF] which describe brain interactions between two nodes i, j in different ways. Here, we use spectral coherence:

Output: J n , S n , n = 0, • • • N -1 1:
Step 1: Transform computation a: Compute the square root matrix

Q 0 ← K -1/2 0 and the eigenvectors V 1 and the eigenvalues σ 2 0 , • • • σ 2 N-1 of the eigen- decomposition Q H 0 K 1 Q 0 = V 1 diag σ 2 0 , • • • σ 2 N-1 V H 1 b: Compute T as ← V H 1 Q H 0 and Σ ← (eig(Q H 0 K 1 Q 0 ) ) 2: Step 2: J-divergence computation a: Define a threshold θ b: Compute J n , n = 0, • • • N -1 as J n ←          2|η 2 n |, ⇐⇒ σ 2 n > θ ∪ |σ 2 n -1| > θ σ n -σ -1 n 2   1 + |η 2 n | σ n σ n + σ -1 n σ n -σ -1 n 2    , otherwise
S n = ∑ n J n • t nn ∑ k t nk C ij (ω k ) = | Pij (ω k )| Pi (ω k ) • Pj (ω k ) (4.20)
In Eq. (5.11), Pi (ω k ), Pj (ω k ) and Pij (ω k ) are the the estimated auto-spectra and cross-spectrum of the signals y i [νT s ], y j [νT s ] at the nodes i and j, computed at the frequency bin * ω k = 2π N s k. Given C ij (ω k ) as in Eq. (5.11), the adjacency matrix Â, is estimated by averaging across the N s frequency bins as follows:

Âij = N s -1 ∑ k=0 C ij (ω k ) (4.21) 56

Results on synthetic data

The proposed signal model for synthetic data generates a simple graph connectivity under the two hypotheses H 1 and H 0 . This model successfully simulates graphs characterized by distinct connectivity states in presence of controlled perturbations. Once we have computed the adjacency matrix under H 1 and H 0 , we derive the estimated Laplacians as in Eq. (5.9) and then, we decompose it with its eigenvalues and eigenvectors as in Eq. (4.4). In order to visualize the eigenvectors' behaviour on graph, we represent in Fig. 4 the Laplacian denoising framework. In this direction, we produce 20 random repetitions (or trials) of synthetic SoGs for each statistical hypothesis, as i.i.d. realizations of our model with a fixed set of parameters.

Subspace robustness on synthetic data

In this subsection, we investigate the robustness of Laplacian denoising based on U L ∪ U H , shortly denoted as U L∪H , and we compare it with other subspaces, notably U ALL , U L and U H .

In order to quantify the subspace robustness, for each subspace (ie. U L∪H , U ALL , U L and U H ) we take into consideration two cases, namely the absence and the presence of perturbation, which we indicate as the ground truth (GT) and the noisy cases, respectively. Each GT subspace is compared to several noisy configurations, corresponding to σ w = 0, 1.2 for noise and σ b = 0, 2 for polarization. To measure subspace robustness on synthetic data, we compute the Frobenius subspace distance F [START_REF] Baksalary | On subspace distances determined by the frobenius norm[END_REF] between the GT case and the noisy configurations, changing the perturbation levels * . We report results of Frobenius distance analysis in Fig. 4.6. We plot F as function of trials for the different subspaces, specifically U H (red), the subspace U L (green), and the subspace U L∪H (blue) in several perturbation scenarios. In Fig. 4.6a) the GT case is represented, in which, not surprisingly, F = 0 for every subspace and every trial. If we gradually increase perturbations (i.e. only noise in Fig. 4.6b) or only polarization in Fig. 4.6c)), the most favourable case is U H for almost every trial. When perturbations dramatically increase Fig. 4.6d), performances decrease in particular for U H configuration. In this figure, we do not have results for U ALL case because F = 0 for all the trials and independently from perturbations.

It is then clear that the eigenvectors in U H are more robust compared to the others. This result is not surprising, since in classical signal processing U H eigenvectors are widely used because of their advantages in terms of signal-to-noise-ratio (SNR). Besides, the subspace U L∪H maintains the robustness, while being relevant to describe the network topology In the following analyses. In the folloowing, we verify that the Laplacian denoising based on the subspace U L∪H leads to more separate connectivity states in presence of perturbations.

J-divergence computation on synthetic data

Finally, we test the ability to separate graph Laplacians under H 1 and H 0 . The J-divergence analysis in Section 4.5 ends with a measure of the statistical distance J between two states. In the following subsection, we apply this analysis to graph Laplacians derived from synthetic data in order to compare the discriminant ability of the proposed denoising method with respect to the other subspace configurations, i.e. U ALL , U L and U H .

In several perturbation configurations, we compute the total J as a measure of statistical distance between the two conditions i.e. related to the hypothesis that the Laplacian matrix comes from H 1 or H 0 ) and we estimate the marginal J n as measure of the contribution of each n-variable to the total separability. Table 4.2 contains J-divergence values obtained for several perturbation levels and for the different subspaces. Results show that in absence of perturbations the most favourable subspace is U ALL . This result is quite intuitive because without perturbations there is no reason why reduced subspaces should better discriminate. If perturbations (i.e. noise and polarization) increase, the most favourable case becomes U L∪H , which gives the highest J. It means that graph Laplacian denoising based on U L∪H offers the highest separability between the two connectivity states, even in presence of high perturbation. representations show the contribution of n variables to the total J and we find that increasing perturbations, variables in U L∪H generally give higher J n contributions compared to U ALL .

Table 4.2: J-divergence values on synthetic data. We report in bold characters the highest J-divergence value for each perturbation configuration.

To conclude, our results with synthetic data demonstrate that in presence of perturbations the Laplacian denoising algorithm succeeds in distinguishing graphs under two conditions. This consideration remains true if the system is perturbed by noise but also if there is an artefact of a different nature , i.e. a common artefact that we indicated as polarization.

Real BCI measurements

In this section, we apply the Laplacian denoising procedure on real data, recorded during motor-imagery BCI experiments. In this case the H 1 and H 0 hypotheses directly correspond to the hypotheses that subject performs motor imagery (H 1 )

or he/she is in resting state (H 0 ).

Experimental Protocol and Preprocessing

Twenty healthy subjects (aged 27.60 ± 4.01 years, 8 women), all right-handed, participated to the study. They received financial compensation for their participation and signed a written informed consent. The ethical committee CPP-IDF-VI of Paris approved the experimental protocol. During the experiments, subjects were seated in front of a screen with a target. They were instructed to perform a right hand-based motor imagery task ( i.e. grasping movement ) when the target was up, and to rest when the target was down [START_REF] Wolpaw | The wadsworth center brain-computer interface (bci) research and development program[END_REF]. EEG data were recorded with 74 channels in a standard 10-10 configuration.

The reference for EEG recordings was set to mastoids and the ground electrode was located on the left scalpula; the impedences were lower than 20 kOhms. The sampling frequency was originally 1 kHz, and then downsampled to 250 Hz. For each subject, EEG data have been segmented to obtain N T trials for motor imagery and N T trials for resting state. The length of each trial was 5s.

A preliminary preprocessing analysis was performed. More precisely, an Independent Component Analysis (ICA) with the Infomax Algorithm [START_REF] Bell | An information-maximization approach to blind separation and blind deconvolution[END_REF] was applied to original data to eliminate artefacts, such as ocular and cardiac signals [START_REF] Delorme | Enhanced detection of artifacts in eeg data using higher-order statistics and independent component analysis[END_REF].

J-divergence of brain connectivity states

We perform the J-divergence analysis on real EEG data. To this aim, we consider EEG signals from one subject and N T trials for H 1 and N T trials for H 0 , with N T = 20. We use spectral coherence to estimate the connectivity matrix, as in Eq.

(5.11). Then, we derive the estimated adjacency matrix  as in Eq. (4.21), and through the Eq. (5.9), we can compute L. As for synthetic data, we derive the filtered graph Laplacian L with U L∪H . In order to compare results of denoising procedure applied to real data, we also consider the U ALL , U H , U L subspaces.

In each case, we estimate the J-divergence J as in Eq. (4.11) and the marginal contribution J n associated to the n-th variable as in Eq.(4.16).

Table 4.3:

J-divergence values obtained with real EEG data. We highlight with bold characters the highest J-divergence value.

In Table 4.3, we report J-divergence values related for each subspace configuration. Comparing all the cases, the highest J-divergence value is obtained with U L∪H case. This result is very important because it means that the subspace U L∪H is appropriate to separate real EEG data and it is useful to correctly identify the subject's mental state.

In order to understand the contributions of different variables, we firstly compute the J n marginal contributions to associate a weight to each and every variable in the transformed domain. Then, we compute the cumulative sum of the first n variables. Specifically, once the J n vector is sorted, we estimate the cumulative J-divergence CJ n to quantify the impact of the variables to the J-divergence:

CJ n = n ∑ k=1 J k (4.22)
Results in Fig. 4.8 demonstrate that the cumulative sum CJ n of the first 20 variables is generally higher for U L∪H than all the other subspace cases. Express differently, if a given number of variables are retained, the overall achieved J-divergence is always larger using the proposed denoising algorithm. This result confirms the improvement to the discrimination of the two cognitive states achieved by the proposed Laplacian denoising.

Scoring of Laplacian coefficients in beta band

Given that the denoising with U L∪H allows a better discrimination between motor imagery and resting state, we now investigate the score performances to identify which Laplacian coefficients mostly contribute to separate the two mental states.

To proceed, it is important to remark that the brain response to cognitive tasks is in general not uniform across the frequencies, but it is mostly evident in α (8-13

Conclusion and further work

conditional (H 1 , H 0 ) estimated adjacency matrix  as in Eq.(5.12), the estimated graph laplacian L as in Eq. (4.4), and its denoised version L as in Eq.(4.8). Then, we evaluate the J-divergence between the two hypotheses as in Eq. (4.16). Finally, we average the J obtained on the m-th window m = 0 • • • M -1 in each timeinterval across subjects. For comparison's sake, we repeat the above computations in absence of denoising (i.e. with U ALL subspace). We plot the J-divergence over M = 9, 1s long, time intervals with 50% overlapping, versus the time interval index. The J-divergence is computed in β band and averaged across subjects.

Our findings on real EEG data show that the proposed Laplacian denoising applies also on short time-windows and improves the detectability of motor imagery states.

Conclusion and further work

In this chapter, we proposed a Laplacian denoising algorithm for the purpose of connectivity states detection. A novel formulation of the Jensen divergence has been derived to measure the performance of the denoising algorithm, as well as to attribute a score to the Laplacian coefficients. The Laplacian denoising performances are assessed by numerical simulations on synthetic data. Furthermore, the Chapter 4. Improving graph connectivity estimation with graph signal processing 65

Laplacian denoising algorithm has been applied to real EEG data recorded during motor imagery BCI experiments. Our results show that the proposed denoising strategy improves the separation of the two cognitive states of motor imagery and resting, even under short time intervals. In addition, the J-divergence based scoring highlights the contribution of Laplacian coefficients to the separability between two cognitive states. A critical aspect in the actual implementation of our framework will be the choice of parameters, such as the number of eigenvectors in each subspace or the frequency band of interest. One possibility to identify adapted parameters, consists in defining them for each specific subjects during the BCI training. Taken together, the proposed approach is promising for the robust detection of connectivity states while being appealing for application in real-time BCI systems.

Appendix E: Theorem 1

Let us consider the problem of binary classification of Gaussian variables H 0 :

x ∼ N (0, I), H 1 : x ∼ N η, Σ 2 , corresponding to the uncommon mean, uncommon covariance case, by means of the LLRT formulation in Eq.(4.12). By simple algebraic manipulation, we recognize that the test R(x)

H 1 ↑ ≷ ↓ H 0 t corresponds to: R (x) = P ∑ n=1 σ 2 n -1 σ 2 n x n + η n σ 2 n -1 2 P quadratic terms + N ∑ n=P+1 2η n • x n N -P linear terms H 1 ↑ ≷ ↓ H 0 t (4.24) with t = t + ∑ P n=1 |η 2 n | σ 2 n σ 2 n -1 -1 .
Let us consider the linear-quadratic observation space Ξ of the N-dimensional random vector ξ def = [ξ 1 . . . ξ N ] T defined as (see Eq. (4.15))

ξ n = x n + η n σ 2 n -1 2 = x 2 n + 2 x n η n σ 2 n -1 + η n σ 2 n - 1 2
; n = 1, . . . , P

ξ n = x n , n = P + 1, . . . , N (4.25) 66 4

.8. Conclusion and further work

In the space Ξ the LLRT R (x)

H 1 ↑ ≷ ↓ H 0
t rewrites as follows:

N-1 ∑ 0 a LLR,n ξ n = a H LLR • ξ H 1 ↑ ≷ ↓ H 0 t (4.26)
where the elements of a LLR def = [a LLR,1 , . . . , a LLR,N ] T are:

a LLR,n def =      σ n -σ -1 n σ n for n = 1, P 2η n for n = P + 1, N (4.27)
With these positions,

J def = E (R(x)|H 1 ) -E (R(x)|H 0 ) = E R (x)|H 1 -E R (x)|H 0 = N-1 ∑ n=0 a LLR,n (E (ξ n |H 1 ) -E (ξ n |H 0 )) (4.28)
By computing the above expectations it can be straightforwardly shown that the n-th term a LLR,n (E (ξ n |H 1 ) -E (ξ n |H 0 )) of the above sum equals to their moments change over time [START_REF] Boveroux | Breakdown of within-and between-network resting state functional magnetic resonance imaging connectivity during propofol-induced loss of consciousness[END_REF][START_REF] Chang | Time-frequency dynamics of resting-state brain connectivity measured with fmri[END_REF][START_REF] Schrouff | Brain functional integration decreases during propofol-induced loss of consciousness[END_REF]. The development of a robust tool to estimate graph functional connectivity is crucial for BCI applications, where a short-time feature extraction is needed to realize a real-time interaction between user and machine [START_REF] Wolpaw | Brain-computer interfaces for communication and control[END_REF]. L1-Principal Component Analysis (L1-PCA) is a stateof-the-art computational tool able to extract relevant components from data with noise and outliers [START_REF] Markopoulos | Adaptive l1-norm principalcomponent analysis with online outlier rejection[END_REF]. The idea of applying L1-PCA to capture the evolution of time-varying data has already been explored, for instance to track time-varying L1-PCA components [START_REF] Markopoulos | L1-norm principal-component analysis via bit flipping[END_REF]. A moving window is usually applied for tracking and it gains in computational cost because only the innovation components are computed at each new time interval. The disadvantage is that this method requires a preliminary selection of the window length, which is defined according to the time duration of the dynamic evolution to be analysed. Another possibility consists in applying a novel hierarchical approach, referred to as Deep L1-PCA. A deep computational architecture has already been investigated to extract Euclidean norm PCA in the context of face recognition in [START_REF] Liong | Face recognition using deep pca[END_REF], even though the impact of layering parameters was not discussed. Our Deep L1-PCA algorithm is composed by the sequence of two recursive steps: data partitioning and L1-PCA analysis with first rank component extraction.

J n (σ n ,η n ) = σ n -σ -1 n σ n σ 2 n + η 2 n + 2 η 2 n σ 2 n -1 -1 = σ n -σ -1 n σ n -σ -1 n + η 2 n σ n σ 2 n + 1 σ 2 n -1 = σ n -σ -1 n 2   1 + η 2 n σ n σ 2 n + 1 σ n -σ -1 n (σ 2 n -1)   = σ n -σ -1 n 2   1 + η 2 n σ n σ n + σ -1 n σ n -σ -1 n 2    , n = 1, • • • P J n (η n ) = 2η 2 n , n = P, • • • N -1.
In this chapter, we firstly describe the mathematical formulation of the Deep L1-PCA method and its advantages in minimizing the effect of perturbations with decreasing computational cost. This approach is applied here to obtain an estimation of graph connectivity in short-time windows and to define a classification scheme to discriminate two connectivity states.

For sake of clarity, we report in Table 5.1 the list of the main notation used in this chapter. procedure can be repeated as many times as necessary to end up, eventually, with only one principal component, which we refer to as global L1-PC. For sake of clarity, we remark that from the second layer on, the input data are the unit-norm W j . We assume that iterative partitions of the data lead to a single global component, i.e. M = ∏ K k=1 W k . This hypothesis is necessary because in some cases, the size M of the original dataset cannot be divided by the chosen W k obtaining equal batches. As consequence, a number of elements W k -1 have to be separately analysed in some manner.

Deep L1-PCA computational framework

At the m-th batch and k-th layer, the L1-PC is computed using the fast estimator introduced by Markopoulos et al. [START_REF] Markopoulos | Efficient l1norm principal-component analysis via bit flipping[END_REF]:

r k (m) = R k (m)b k (m) • ∆ k (m) (5.1)
b k (m) is the optimal binary vector, computed as:

arg max b∈{±1} W k R k (m)b 2 (5.2)
and ∆ k (m) comes from the equation that follows:

∆ k (m) = 1 R k (m)b k (m) 2 (5.3)
As from the equations (5.1) and ( 5.

3) , the r k (m) is a normalized version of R k (m) b k (m).
These operations provide a set of M k principal components, where each of them is the first L1-PC related to one batch m k . The iterative repetition of this procedure finally reaches the last layer K.

The Deep L1-PCA can be written in a compact manner using the Khatri-Rao product * To this aim, we define the matrices:

B k = [b k-1 (1), b k-1 (2), . . . , b k-1 (m), . . . , b k-1 (M k-1 )] ∈ {±1} W k-1 ×M k-1 which contains the optimal binary vectors b k-1 (m) ∈ R W k -1 R k = [r k-1 (1), r k-1 (2), . . . , r k-1 (m), . . . , r k-1 (M k-1 )] ∈ R D×M k-1 which collects the M k-1 L1-PCs computed in k -1 for k > 1 and for k = 1 it contains the original data R 1 = Ξ ∆ k = diag([∆ k-1 (1), ∆ k-1 (2), . . . , ∆ k-1 (M k-1 )])
Applying the Khatri-Rao product to B k ∈ R W k-1 ×M k-1 , we derive

I M k-1 B k =                      b k-1 (1) 0 W k-1 . . . 0 W k-1 0 W k-1 b k-1 (2) . . . 0 W k-1 . . . . . . . . . . . . 0 W k-1 0 W k-1 . . . b k-1 (M k-1 )                     
We write the whole set of projectors at the first layer as follows:

R 2 = [r 1 (1), . . . , r 1 (M 1 )] = [R 1 (1), . . . , R 1 (M 1 )] • (I M 1 ×M 1 B 2 ) • ∆ 2 = Ξ • (I M 1 ×M 1 B 2 ) • ∆ 2 R 2 D×M 1 = Ξ D×M ( I M 1 ×M 1 B 2 W 1 ×M 1 M 1 • W 1 × M 1 = M × M 1 ) • ∆ 2 M 1 ×M 1 (5.4)
and, for the other k = 2, 3, . . . , K layers:

R k = R k-1 • (I N k-1 B k ) • ∆ k ∈ R D×W k-1 (5.5)
We can write R k as:

R k = Ξ • (I M 1 ×M 1 B 2 ) • ∆ 2 • . . . • (I M k-1 ×M k-1 B k ) • ∆ k = Ξ • k ∏ j=2 [(I M j-1 ×M j-1 B j ) • ∆ j ]
(5.6)

Application to graph synthetic data

At the final layer K, the global principal component can be computed:

r K = R K • (I M K ×M K B K+1 ) • ∆ K+1 (5.7)
where

M K = 1, B K+1 = b K (1) ∈ {±1} M K-1 ×1 and ∆ K+1 = ∆ K (1) ∈ R + .
Some considerations are in order: i) L1-PCA differs from the overall L1-PCA computed on the whole original data; ii) Deep L1-PCA intrinsically has an additional information, which is the relative distance of L1-PCA solutions computed at the intermediate layers.

Application to graph synthetic data

In this section, we test the ability of the Deep L1-PCA to characterize two estimated graph Laplacian matrices, related to two connectivity states. We test our approach measuring the mean-square error (MSE) -to quantify the robustness of the procedure-and implementing a classification scheme -to verify the ability of the L1-PCs to discriminate two connectivity states.

Connectivity matrices simulation

In order to validate our method on synthetic data, we define two adjacency matrices under the two hypotheses H 0 and H 1 .

Under H 1 , we define an adjacency matrix characterized by all zero connections, except for H c connections set at 1, to simulate strong links. A noise component is added to the adjacency matrix in order to simulate perturbations. Under H 0 , H c = 0 (i.e. there are not strong connections) and there is only a noise contribution. We can write the estimated adjacency matrix Â(m) for each observation m as the sum of the real adjacency matrix A(m) and the noise contribution w(m):

H 1 : Â(m) = A (1) (m) + w(m) H 0 : Â(m) = A (0) (m) + w(m) (5.8)
where w is the noise, which is a realization of a discrete, stationary, white Gaussian process, with E{w} = 0, E{ww T } = σ 2 w I ∀k; A (0) is a zero matrix; A (1) is a zero-matrix containing ones for h C = 0, • • • H C -1 elements. Once obtained the synthetic adjacency matrices, we can derive the estimated Laplacians:

L(m) = D(m) -Â(m) (5.9)
where D is the estimated degree matrix, whose generic diagonal element d i is equal to the sum of the weights of all the edges incident to the node i. For each observation, and under each hypothesis, we derive the Laplacian matrix L(m).

For the following analysis on those synthetic data, we consider M = 8 observations and σ w = 0.6.

Robustness analysis via MSE

In the following, we will verify the robustness of the Deep L1-PCA. In order to do so, we will compute for each layer k and each observation m, the mean squared error between PCs obtained with synthetic data generated under two conditions :

(1), the presence of noise; (1)REF, the reference case without perturbations. Specifically, we consider the data matrices at each m observation ξ (1) . We collect M= 8 degree matrices for each condition.

Following the procedure described in section 5.2, we perform two Deep L1-PC analyses with W K = 2 and M k = 1, for k = 1, and we obtain: i) R

(1)REF k from Ξ (1) and ii) R

(1)REF k from Ξ (1)REF .
This procedure provides two Deep L1-PCA structures, one computed from data corrupted by noise and the other obtained from data without perturbations.

For each set of original data, a general scheme of the Deep L1-PCA procedure is presented in Fig. 5.4.

After a preliminary normalization step, for each observation m and each layer k, we compute the mean squared error (MSE) as follows: Let us remark that the definition of the MSE is expressed as in (5.10) (i.e. with the alternative on the sign of the PC estimates) because the direction of principal components is not uniquely defined through the procedure in [START_REF] Markopoulos | Adaptive l1-norm principalcomponent analysis with online outlier rejection[END_REF].

MSE k (m) = min(||R (1)REF k (m) -R (1) k (m)|| 2 , ||R (1)REF k (m) + R (1) k (m)|| 2 ) (5.10)
Results of the MSE computation are in Fig. 5.5, where MSE values are reported for each layer k and each observation m in a color scale. Our findings show that MSE values generally decrease by considering higher layers 19 to 9. This results is perfectly in line with a previous analysis in [START_REF] Orrú | Deep ll-pca of time-variant data with application to brain connectivity measurements[END_REF].

Since MSE values are computed between PCs obtained on noisy data with respect to those in absence of perturbations, results demonstrate that the impact of noise generally decreases in higher layers.

Classification framework

Here, we investigate the possibility to apply the Deep L1-PCA to obtain robust features to detect two connectivity states corresponding to H 1 and H 0 .

In order to perform a classification procedure, we synthetically produce N T repetitions, or trials, to obtain N T adjacency matrices for m = 1, .., M under each condition. In this manner, for each condition, we collect N T trials for the M adjacency matrices to allows a classification analysis for each layer k and each observation m. After matrices' generation, we follow the procedure presented in section 5.2

for each trial, in order to obtain for each condition a N T -dimensional set of R (1) k and N T -dimensional set of R (0) k . R (1) k is the Deep L1-PCA estimation when the original data matrix corresponds to ξ In this study we perform one classification test for each observation (or batch) m and each layer k in order to quantify the discrimination power associated to each specific case. More in detail, we implement a classification using a 10 repeated 10-fold LDA classifier for each batch and each layer, with N T = 100 for each connectivity state.

Results are presented in Fig. 5.6A), where we report the classification accuracy values. Interestingly, the accuracy, which randomly varies across batches of the same layer, generally increases considering higher layers (from k = 1 to k = 3). Specifically, results on the first layer oscillate around 65% while these values increase until reaching 78% in the final layer with k = 3.

In this first classification scenario, the feature vectors consist of the degree matrices under the two conditions. Another possibility is to build the feature vectors by means of the entire Laplacian matrices, taking into consideration both node strength and link values.

More formally, for each of the N T trials, the R (1) k under H 1 and R 0 k H 0 are obtained from ξ Our results on synthetic data show that in presence of perturbations Deep L1-PCA allows to separate graphs under two conditions. We reach this conclusion through two steps: i) a preliminary analysis of the mean squared error between a noisy graph and its ground truth counterpart to verify the robustness to noise of Deep L1-PCA; ii) a classification test to demonstrate the potential of Deep L1-PCA to separate two connectivity states. Table 5.2: Classification accuracy for synthetic data related to the original M = 8 observations. In the first row, we have accuracies when only the Laplacian diagonal elements are taken as classification features. In the second row, accuracy values relates to when all the Laplacian elements are considered.

Results on real BCI data

In this section, we present experimental results on the application of the Deep L1-PCA framework to EEG data recorded during motor imagery experiments.

Deep L1-PCA offer new tools in BCI data analysis in order to better identify mental states. For this analysis, H 1 and H 0 directly correspond to the hypotheses that subject performs motor imagery (H 1 ) or resting state (H 0 ).

Experimental Protocol and Preprocessing

The experimental protocol is constituted by twenty healthy subjects (aged 27.60 ± 4.01 years, 8 women), all right-handed. The study was approved by the ethical committee CPP-IDF-VI of Paris. The subjects, who did not have any physical or psychological disorder, received a financial compensation for their participation and they signed a informed consent.

During the experiments, the subject was seated in front of a screen, with a target. Specifically, when the target was up, subject has to perform a right hand motor imagery task and when it was down, he/she had to remain at rest [START_REF] Wolpaw | The wadsworth center brain-computer interface (bci) research and development program[END_REF] .

EEG data were recorded with a 74-channel system in a standard 10 -10 configuration. The sampling frequency was originally 1 kHz, and then downsampled to 250 Hz. Each subject performed several time the task, obtaining N T trials of motor imagery and N T trials of resting state. The trial length was 5s.

Before the following analyses, data were preliminarily preprocessed. An Independent Component Analysis (ICA) , through Infomax Algorithm [START_REF] Bell | An information-maximization approach to blind separation and blind deconvolution[END_REF] was performed with Fieldtrip toolbox [START_REF] Oostenveld | Fieldtrip: open source software for advanced analysis of meg, eeg, and invasive electrophysiological data[END_REF] to eliminate artifacts, such as ocular and cardiac signals [START_REF] Delorme | Enhanced detection of artifacts in eeg data using higher-order statistics and independent component analysis[END_REF].

As a proof of concept, we only show results associated to one subject and N T = 74.

Functional connectivity estimation procedure

In this subsection, we aim to demonstrate that the Deep L1-PCA is suitable to identify time-varying data related to two cognitive states.

The first step necessary to apply the Deep L1-PCA computation to brain EEG data consists in building the graph to finally estimate a Laplacian matrix. In this specific application, EEG data recorded at each trial directly constitute Signals on Graphs (SoGs).

In order to obtain M Laplacians for each condition, we split the total length of each trial into M time windows of length 500 ms each, obtaining 10 batches. Without loss of generality, we consider only M k = 8 for k = 1 for sake of simplicity, but alternative strategies can be identified to manage the last 2 windows. For instance in real applications, the number of the original batches (i.e. time windows) can be tuned to not loose data and even an overlapping factor can be applied to increase the number of windows.

Once SoG samples belonging to each m window are identified, the adjacency matrices can be estimated. Among the state-of-the-art methods able to derive the estimation of the adjacency matrix, we use spectral coherence [START_REF] Carter | Coherence and time delay estimation[END_REF], which has already exhibited advantages in motor-imagery applications [START_REF] Cattai | Characterization of mental states through node connectivity between brain signals[END_REF][START_REF] Cattai | Combination of connectivity and spectral features for motor-imagery bci[END_REF]. It is computed for each pair of nodes i, j as follows, obtaining one estimate for each trial and each batch:

C ij (ω k ) = | Pij (ω k )| Pi (ω k ) • Pj (ω k ) (5.11)
In Eq. (5.11), Pi (ω k ), Pj (ω k ) and Pij (ω k ) are the estimated auto-spectra and cross-spectrum of the signals at the nodes i and j, computed using Welch method at the frequency bin ω k = 2π N s k.

In order to derive the adjacency matrix it is possible to average C ij (ω k ) across all the N S frequency bins or to consider only some bins of interest depending on the application. Our strategy is to perform the analysis in precise frequency bands.

Regarding the choice of the appropriate frequency band of interest, we underline that the brain activity in response to motor tasks varies in the frequencies and it is generally more evident in in α (8-13 Hz) or β (14-29 Hz) bands depending on the subject and on the task [START_REF] Meirovitch | Alpha and beta band event-related desynchronization reflects kinematic regularities[END_REF].

Here, we perform the analysis in β band, knowing that in a real BCI scenario,

Results on real BCI data

the frequency band or the frequency bin of interest can be tuned during the calibration phase according to the subject characteristics. We compute the adjacency matrix by filtering the connectivity matrix in the selected frequency band:

Âij (m) = ∑ ω k /T s ∈β C ij (ω k ) (5.12)
After the computation of the adjacency matrix for each trial and each time window, we estimate the Laplacian matrices L as in (5.9) for motor imagery (H 1 )

and for resting state (H 1 ). The Deep L1-PCA structure can be computed as in sect.5.2, obtaining N T repetitions of R

k and N T repetitions of R

k .

Classification analysis on real EEG data

In this subsection, we verify the ability of the Deep L1-PCA procedure to identify the two mental states of motor imagery and resting in an off-line classification scenario.

Similarly to synthetic data, we perform a a 10 repeated 10-fold LDA classification for each batch and each layer. In Fig. 5.7 we report the associated results. In presents findings obtained when the feature vector is build with all the Laplacian matrix. In both the cases, we can recognize a general accuracy improvement considering progressively higher layers from k = 1 to k = 3. Classification results generally improve when also extra-diagonal Laplacian elements are considered with respect to those obtained with only diagonal elements (maximum value of 82% compared to 75%).

For sake of completeness, we report in Table 5. Table 5.3: Classification accuracy for real EEG data related to the original M = 8 observations. In the first row, we have accuracies when only the Laplacian diagonal elements are taken as classification features. In the second row, accuracy values relate to when all the Laplacian elements are considered as features.

Deep L1-PCA applicability in BCI systems

Deep L1-PCA applicability in BCI systems

BCIs require real-time interaction between user and interface [START_REF] Wolpaw | Brain-computer interfaces: Definitions and principles[END_REF], which implies time-varying computation of the BCI features. Deep L1-PCA provides a useful tool to have robust short-time estimates of graph Laplacians with low computational cost. Novel techniques have been recently developed to capture temporal evolution of graphs. Tracking algorithms has been proposed in [START_REF] Shen | Tensor decompositions for identifying directed graph topologies and tracking dynamic networks[END_REF] and [START_REF] Ozdemir | Recursive tensor subspace tracking for dynamic brain network analysis[END_REF] under the hypothesis of low rank temporal dynamics of FC. Other state-of-the-art strategies require smoothness on time dimension to solve advanced regularization problems [START_REF] Kalofolias | Learning time varying graphs[END_REF][START_REF] Yamada | Time-varying graph learning based on sparseness of temporal variation[END_REF]. Besides, techniques based on the estimation of spatio-temporal graphs describe the connectivity patterns across time but often the computation cost is really high [START_REF] Ortiz-Jiménez | Sampling and reconstruction of signals on product graphs[END_REF][START_REF] Romero | Kernel-based reconstruction of space-time functions on dynamic graphs[END_REF]. The advantage of Deep L1-PCA with respect to the other methods is that it is naturally suitable to deal with noisy data, it does not require any specific assumption and it has reduced computational cost.

We propose a scheme to actually apply the Deep L1-PCA to BCI scenarios.

In order to understand the potential of the Deep L1-PCA in BCIs, let us imagine to perform an online experiment, where time samples are progressively available and we suppose to have already trained the classifier. Once we have received all the samples in the first m time-interval, we can estimate functional connectivity and derive the Laplacian matrix. After first computation, the first classification can be done to give a feedback to the user. As the time passes, we collect other batches and we start to build the Deep structure. For instance, once the second Laplacian at m = 2 is estimated, we can compute the PC to obtain r 1 [START_REF] Ahn | Performance variation in motor imagery braincomputer interface: a brief review[END_REF]. We can continue this procedure until the end of the trial, when we have the whole Deep structure. Since we verified that accuracy increases with higher layers, it is more convenient to perform classification at the highest available layer. The possibility of considering high layers depends on the available data already collected. In This procedure enables to simultaneously i) perform classification at each new short-time interval to immediately interact with the subject with a possible feedback ii) build the Deep structure to collect all the layers to obtain higher layer estimations, which can be progressively used to detect subject's mental state.

Let us remark that for this analysis we have one Laplacian computation each 500 ms but modifications can be done to have more frequent estimations. In real BCI applications, parameters can be tuned to change the length of the timewindows or to apply an overlap. These results pave the way for the application of Deep L1-PCA framework to EEG data to online control BCI systems.

Conclusions

We proposed a novel formulation of the L1-PCA, to which we refer to as Deep This method offers the possibility to realize a fast interaction between user and machine by improving the classification at higher layers. More research is needed to test this method in real-time applications as well as to compare it with other state-of-the-art techniques.

Comparing Laplacian denoising in Chapter 4 with Deep L1-PCA, it is clear that while the first requires training and a maximum likelihood estimator to optimally separate data, the latter is simple and fast and it is defined in the linear space. LWe remark that the Laplacian denoising algorithm can be improved for BCI purpose, by means of covariance matrix tracking, or learning techniques to

Conclusions

avoid the estimation of the transformation parameters [START_REF] Ali | Biometricnet: deep unconstrained face verification through learning of metrics regularized onto gaussian distributions[END_REF]. One possibility to benefit of the two approaches is to integrate them to have a robust estimation of graph Laplacian and, simultaneously, to preserve the advantages of Deep L1-PCA. This can be done at the classification level, using fusion procedures [START_REF] Corsi | Integrating eeg and meg signals to improve motor imagery classification in brain-computer interface[END_REF][START_REF] Ruta | An overview of classifier fusion methods[END_REF] or the two methods can be integrated to obtain another version of the graph Laplacian matrix, for instance, combining eigenvectors. In real BCI applications, an appropriate way to identify the most suitable features and the most appropriate strategy to obtain robust FC estimates for a specifc user consists in performing preliminary analyses during the calibration phase .

Chapter 6

Conclusions and future perspectives Conclusions

In this thesis, we investigated the potential of functional connectivity networks to characterize BCI related tasks. We proposed an original framework to adapt graph theory, signal processing and graph signal processing tools in order to have insight in brain processes during cognitive tasks.

By applying functional connectivity and network estimators to characterize brain states, we uncovered functional interaction changes between graph nodes during motor imagery tasks.

In addition, we demonstrated that including FC-based features in the classification generally performs better than standard techniques only in detecting the user's mental state [START_REF] Miller | Cortical activity during motor execution, motor imagery, and imagery-based online feedback[END_REF]. Investigating motor imagery tasks, we discovered the presence of a twofold FC mechanism based on phase and amplitude synchronization between EEG signals.

We developed a novel algorithm to robustly estimate graph Laplacian and to improve the separability between two brain states. We introduced a denoising procedure to obtain a more robust graph FC estimation, based on tools from different domains, such as information theory and graph signal processing. In this context, the separability between two brain states was quantified by an original formulation of the Jensen divergence, which highlights the contribution of each Laplacian element to the final discrimination.

Finally, we presented a novel framework to obtain short-time estimates of graph connectivity using an innovative tool, named Deep L1-PCA. This method is able to capture the dynamic of time-varying signals and to rapidly identify the 86 user's mental state with an iterative robust procedure, particularly appropriate for BCI applications.

Taken together, our findings pave the way for the development of alternative BCI systems, based on FC-related features. However, it is important to underline that particular attention should be devoted to the practical implementation of the whole framework. There are some parameters to be set, such as the number of eigenvectors needed to define the different subspaces for the Laplacian denoising algorithm. In the J-divergence analysis, a threshold has to be fixed to separate the Laplacian variables changing only in mean value from those changing also in variance. For what concerns the L1-PCA computation, the number of elements in each batch has to be set and it influences the number of observations necessary to perform the PCA analysis in higher layers. Since these values are intrinsically dependent on the subject and the task, we recommend to train the algorithms and set all the parameters before testing the experiments. This training step can be, for example, performed during the BCI calibration phase.

More in general, we estimated FC networks by computing coherence and imaginary coherence, but several other measures can be adopted to quantify brain interactions [START_REF] Bastos | A tutorial review of functional connectivity analysis methods and their interpretational pitfalls[END_REF][START_REF] Gonzalez-Astudillo | Network-based brain computer interfaces: principles and applications[END_REF][START_REF] Hamedi | Electroencephalographic motor imagery brain connectivity analysis for bci: a review[END_REF]. The choice of the FC estimator can, in principle, depend on the subject, because humans use different strategies to perform motor or cognitive tasks. More research is needed to generalize our entire method to more subjects by using, for instance, different databases [START_REF] Tangermann | Review of the bci competition iv[END_REF]. Although our research demonstrates that FC features capture complementary information with respect to univariate features, the extent to which they can be applied to real BCI systems strongly depends on the subjects. This problem might be solved at the level of BCI design, in fact feature selection procedures can be performed during the calibration phase.

Current work and future perspectives Experimental validation and implementation into the Openvibe software

Our proposed framework provides a practical approach to include graph FC features in the identification and characterization of human mental states. The research in this thesis aims to improve motor imagery-based BCI systems by taking into consideration informative features describing brain interaction mechanisms.

To test the feasibility of our approach we are currently implementing it in Openvibe [START_REF] Renard | Openvibe: An open-source software platform to design, test, and use brain-computer interfaces in real and virtual environments[END_REF], which is a software specifically designed for real-time BCI. It can generally be used to visualize and process (acquire, analyse, classify) brain data. This ongoing work is a collaboration with the Inria engineer Arthur Desbois in the Aramis team.

Markovian models for time-varying graphs

The method introduced in Chapter 5 provided a framework to obtain a shorttime graph connectivity estimation in a robust way, using an iterative algorithm for principal components analysis. This procedure is particularly appropriate for BCI applications because it gives as output both the estimations at each short-

2. 2 3 . 2 3 . 6

 23236 Main steps to estimate graph connectivity networks from time series. The first step is definition of graph nodes (i.e. the yellow circles on the left panel). From each node, brain signals are recorded. After choosing an appropriate FC estimator, functional connectivity is computed for each pair of nodes in order to obtain an adjacency matrix. Nodes and edges, corresponding to the elements in the adjacency matrix, constitute the final brain network. . . . . . . . . . . 2.3 Graphical representation of signals on graph with the underlying graph structure. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4 Example of network changing its edges in time (m) . . . . . . . . . . 3.1 Connectivity estimators as a function of normalized distance between EEG electrodes. Coherence values are in orange while imaginary coherence results are in green. Connectivity estimated are extracted from all the available experimental observations in the β frequency band. [30] ©2018 IEEE . . . . . . . . . . . . . . . . . . . . . Statistical contrast maps between motor imagery and resting states obtained with band power features. Results are shown for onetailed permutation-based t-tests (p<0.01). In Panel A) the obtained t-values are illustrated for individual representative frequency bins within the theta band, B) alpha band, C) beta band, and gamma band. In the beta band, results also remained significant after correction for multiple comparison (p<0.05, FDR corrected). . . . . . . . . . . . 20 3.3 Statistical contrast maps between motor imagery and resting states in the beta band. In panel A, we report the results in beta band for connectivity and node strength values. In panel B, we show the results of node strength for individual representative frequency bins within the beta band. Results for coherence-based features are in top line, while those for imagery coherence-based features are in bottom line. Only the twenty most discriminant values are represented here for illustrative purposes. . . . . . . . . . . . . . . . . . . . . . . . . . . 21 3.4 Relationships of coherence/imaginary coherence with phase difference. In Panel A) coherence is in pink and imaginary coherence in green, showing the functional connectivity between two sine waves at 10 Hz as function of their temporal shift. The shift, corresponding here to a phase difference, varies from 0 to π in steps of π /500. At each shift value, the two connectivity estimators are evaluated. Panel B) shows the sine waves with different phase differences. In panel 1), a positive ∆ of π /4 in panel 2), a negative ∆ of -π/4. . . . 22 3.5 Phase difference properties and discrimination ability. Panel A), results of permutation-based t-tests in the beta band across all subjects are shown for brain networks reconstructed from the phase difference between EEG signals. Panel B) results of permutation-based t-tests obtained with node strength values extracted from the previous brain networks. Panel C), Spearman correlation plot between imaginary coherence and phase difference values considering all the connections including C3 electrode for one representative subject. . 23 List of Figures vii Improvement of classification performance. Bar plots show the percentage of relative increment between the best combination of features (i.e., coherence-based node strength S C , imaginary coherencebased node strength S IC and power spectrum P. The pie diagram in the inset illustrates the percentage of times that a specific combination of features has been selected across subjects.Two different cases are considered, distiguished by a different color code. In blue, we report results of accuracy improvements when selected electrodes are located in controlateral sensorimotor area (Contra), i.e. FC5, FC3, FC1, C5, C3, C1, CP5, CP3, CP1. In red, results relate to classification analysis when electrodes are in both contralateral and ipsilateral sensorimotor areas (Contra+Ipsi), ie. FC5, FC3, FC1, C5, C3, C1, CP5, CP3, CP1, FCZ, CZ, CPZ, FC2, FC4, FC6, C2, C4, C6, CP2, CP4, CP6 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 3.7 Brain features selected by the classification procedure. The color codes for the group-averaged number of times that a specific feature -in the electrode-frequency space -has been chosen during the sequential feature selection algorithm. The results for P features are illustrated in the top line, those for S C in the middle line and those for S IC in the bottom line. Results relates to the classification analysis framework when only electrodes in the contralateral sensorimotor areas are selected. . . . . . . . . . . . . . . . . . . . . . 26 3.8 EEG signals before the preprocessing for one subjects. We show the first 60s of recordings. . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 3.9 Independent component obtain with Infomax algorithm. We show the first 10 ICs and the first 60s of recording. For each component, both the signals and the associated topography are represented to facilitate the preprocessing procedure. . . . . . . . . . . . . . . . . . . 31 3.10 EEG signals after preprocessing for one subjects. We show the first 60s of recording. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 viii List of Figures 3.11 Statistical contrast maps between motor imagery and resting states for connectivity features. In top line, we have results for coherence features and in the bottom line for imaginary coherence ones. In Panel A) results for coherence in theta band, in B) for imaginary coherence in theta band, in C) for coherence in al pha band, in D) for imaginary coherence in al pha band, in E) for coherence in beta band and in F) for imaginary coherence in beta band in G) for coherence in gamma band and in H) for imaginary coherence in gamma band . Only the twenty most discriminant connections are represented here for the sake of simplicity. . . . . . . . . . . . . . . . . . . . . . . . 33 3.12 Statistical contrast maps between motor imagery and resting states for node strength features. In top line, we have results for coherencebased node strength features and in the bottom line for imaginary coherence ones. In Panel A) results for coherence in theta band, in B) for imaginary coherence in theta band, in C) for coherence in al pha band, in D) for imaginary coherence in al pha band, in E) for coherence in beta band and in F) for imaginary coherence in beta band in G) for coherence in gamma band and in H) for imaginary coherence in gamma band . Only the twenty most discriminant connections are represented here for the sake of simplicity. . . . . . . . . . . . . . 34 3.13 FBCSP procedure. EEG signals recorded from N channels are bandpassed in frequency intervals. Then, for each frequency interval, the common spatial pattern is computed. Finally, the posterior probability is estimated in order to weight CSP features to perform a classification fusion to integrate different type of features. . . . . . . 36 3.14 Scheme of the fusion procedure. Each type of feature is computed and used to train a classifier. Then, the posterior probability is computed in each case and a Bayesian approach is used to integrate the different features. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 3.15 Improvement of classification performance. Bar plots show the percentage of relative increment between the best combination of features (i.e., coherence-based node strength S C , imaginary coherencebased node strength S IC and common spatial pattern CSP. The pie diagram in the inset illustrates the percentage of times that a specific combination of features has been selected across subjects. . . . . . . 37 4.1 Example of transformation effect. In a) we have bidimensional Gaussian distributions which differ under mean and covariance matrix . In b) we have data in the transformed domain . . . . . . . . . . 49 List of Figures ix 4.

  2 J-Divergence contributions as function of mean η and standard deviation σ: a) J (σ,η) for variables whose conditional standard deviation differ under H 1 and H 0 , and b) J (η) for variables with invariant conditional standard deviation. . . . . . . . . . . . . . . . . . . . . . . 4.3 Graphic interpretation of the score computed for the first element in the vector l . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.4 Adjacency matrix with synthetic data. Âij is represented under H 1 in A) and under H 0 in panel B). . . . . . . . . . . . . . . . . . . . . . 4.5 Example of eigenvectors on graph. They relate to the Laplacian matrix estimated from synthetic data under H 1 . . . . . . . . . . . . 4.6 Results of Frobenius distance on synthetic data. Several perturbation configurations are represented: in panel a)σ w =0 and σ b =0,in panel b) σ w =1.2 and σ b =0,in panel c) σ w =0 and σ b =2 and in panel d) σ w =1.2 and σ b =2. In the different colors (in the legend) we represent the different subspaces. . . . . . . . . . . . . . . . . . . . . . . . . . . 4.7 Results of J-divergence analysis on synthetic data. Several perturbation configurations are represented: in panel a) σ w =0 and σ b =0, in panel b) σ w =1.2 and σ b =0, in panel c) σ w =1.2 and σ b =2. In the different colors (shown in the legend) we represent the different subspaces for the filtering. . . . . . . . . . . . . . . . . . . . . . . . . . 4.8 Results of J-divergence analysis for real EEG data. We report the cumulative J-divergence CJ n as function of the first 20 variables. In the different colors (shown in the legend) we represent the different subspaces used to filter Laplacian. . . . . . . . . . . . . . . . . . . . . 4.9 Results of score computation for real EEG data. All the eigenvectors are considered (i.e. U ALL ) in panels a,c);and results associated to the proposed denoising algorithm (ie. U L∪H ) in panels b,d). In the first line, score values refer to links (i.e. extra-diagonal elements) and in the second line, they refer to elements in the principal diagonal (i.e. nodes). For sake of clarity, in all the figures we plot the 20 nodes or links with highest score. . . . . . . . . . . . . . . . . . . . . . . . . 4.10 Results of J-divergence analysis over a moving window on real data.

  Deep L1-PCA computation structure. . . . . . . . . . . . . . 5.2 A) Time-variant bidimensional normal data and in B) the underlying probability density function (PDF) orientation. . . . . . . . . . . x List of Figures 5.3 Time-variant data and first order L1-PCA approximation at different layers : layer 1 in A), layer 2 in B) and layer 3 in C). . . . . . . . . . . 73 5.4 Graphical representation of Deep L1-PCA computation with an M = 8 original data Ξ. At each layer k (k = 1, 2, 3) the L1-PCAs r k (m) are computed and they become the input in R k+1 (m) for the next layer. 75 5.5 Mean squared error map computed on simulated data between the degree matrix corrupted by noise with σ w = 0.6 and its reference version (i.e. without perturbations). . . . . . . . . . . . . . . . . . . . 76 5.6 Accuracy of classification framework for synthetic data for each layer and each batch. In the rows we have the different layers and in the columns we have the batches, which are 4 at the first layer, 2 at the second and 1 batch at the final layer. In panel A), we report accuracy results when only the Laplacian diagonal elements are taken as classification features. In panel B) results are obtained when all the Laplacian elements are considered. . . . . . . . . . . . . . . . . . 77 5.7 Accuracy of classification framework for real EEG data associated to each layer and each batch. In the rows we have the different layers and in the columns we have the batches, which are 4 at the first layer, 2 at the second and 1 batch at the final layer. In panel A), we report accuracy results when only the Laplacian diagonal elements are taken as classification features. In panel B) results are obtained when all the Laplacian elements are considered as features. . . . . . 81 5.8 Graphical representation of data recordings. As the time passes, we collect m time-windows in which we can estimate the graph Laplacian. As soon as new time intervals are recorded we can compute the Deep L1-PCs at several layers. Here, we indicate with red arrows the time instant (and the associated time interval) in which

3. 1

 1 Table of main notation in Chapter 3. . . . . . . . . . . . . . . . . . . .3.2 Spearman correlation coefficient between IC and ∆ for links including C3 for each subject. . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3 Average accuracy across cross-validation fold and repetitions is reported the best feature combination for each subject. In the top rows, we have classification results obtained when selected electrodes are located in both in contralateral and ipsilateral sensorimotor areas, ie. FC5, FC3, FC1, C5, C3, C1, CP5, CP3, CP1, FCZ, CZ, CPZ, FC2, FC4, FC6, C2, C4, C6, CP2, CP4, CP6. In the bottom rows, results come from classification framework including only controlateral electrodes, ie.FC5, FC3, FC1, C5, C3, C1, CP5, CP3, CP1. . . . 3.4 Average accuracy across cross-validation is reported for each subject and each combination of feature when CSP is used. . . . . . . . . . . 4.1 Table of main notation in Chapter 4. . . . . . . . . . . . . . . . . . . . 4.2 J-divergence values on synthetic data. We report in bold characters the highest J-divergence value for each perturbation configuration. . 4.3 J-divergence values obtained with real EEG data. We highlight with bold characters the highest J-divergence value. . . . . . . . . . . . . . 5.1 Table of main notation in Chapter 5. . . . . . . . . . . . . . . . . . . . 5.2 Classification accuracy for synthetic data related to the original M = 8 observations. In the first row, we have accuracies when only the Laplacian diagonal elements are taken as classification features. In the second row, accuracy values relates to when all the Laplacian elements are considered. . . . . . . . . . . . . . . . . . . . . . . . . . . xi Chapter 1

Figure 2 . 1 :

 21 Figure 2.1: Graphical representation of the most common functional connectivity estimators used in literature.They are organized in a coloured structure to identify different properties. We report undirected and directed estimators in green and red colours respectively. Then, they can be represented according to their ability to capture bivariate or multivariate interactions, in blue or green boxes. The characterization of possible non linearities is reported in grey boxes.

Figure 2 . 2 :

 22 Figure 2.2: Main steps to estimate graph connectivity networks from time series. The first step is definition of graph nodes (i.e. the yellow circles on the left panel).From each node, brain signals are recorded. After choosing an appropriate FC estimator, functional connectivity is computed for each pair of nodes in order to obtain an adjacency matrix. Nodes and edges, corresponding to the elements in the adjacency matrix, constitute the final brain network.

Chapter 2 . 9 Figure 2 . 3 :

 2923 Figure 2.3: Graphical representation of signals on graph with the underlying graph structure.

Chapter 2 . 11 Figure 2 . 4 :

 21124 Figure 2.4: Example of network changing its edges in time (m)

Fig. 2 .

 2 Fig. 2.4 represents an example of network changing its edges over time windows.

  caused by the increase in amplitude and phase synchronization between the brain signals during MI task. At the individual level, we demonstrate the potential of network connectivity features in an off-line classification framework. Taken together, our findings offer new insights into the oscillatory mechanisms underlying brain network during MI task and open new perspectives to improve BCI performance. In this chapter are reported figures and part of text from a submitted paper: Cattai, T., Colonnese, S., Corsi, M. C., Bassett, D. S., Scarano, G., & De Vico Fallani, F.(2019). Phase/amplitude synchronization of brain signals during motor imagery BCI tasks. arXiv preprint arXiv:1912.02745.[START_REF] Cattai | Phase/amplitude synchronization of brain signals during motor imagery bci tasks[END_REF]

± 4 .

 4 01 * years, 8 women), all right-handed, participated in the study. The subjects were recruited for a BCI training protocol and they did not have any medical or psychological disorder. The study was approved by the ethical committee CPP-IDF-VI of Paris and each subject signed a written informed consent. Each participant received financial compensation for the participation.

*

  mean and standard deviation data were recorded with a 74-channel system, with Ag/AgCl sensors (Easycap, Germany) in a 10-10 standard configuration. The reference for the EEG signals were mastoid signals and the ground electrode was set on the left scalpula. Data were recorded in a shielded room. Impedances were lower than 20 kOhms, the sampling frequency was 1 kHz, then downsampled to 250 Hz. All the subjects were naive BCI users and participated in a training protocol. For each subject we collected 64 trials of motor imagery and 64 trials of resting state. In each trial, the first second corresponded to the inter-stimulus interval (ISI), when a black screen was presented to the subject. During the following 5s, the target appeared on the screen and during this period subjects had to imagine a sustained grasping of their right dominant hand. During the experiments, hand muscular activity was recorded with EMG (electromyogram) to check the presence of involuntary movements during the motor imagery tasks. On-line, the experimenter ensured that subjects were not generating muscular artefacts during the task. Off-line, all the recorded signals have been checked to exclude the presence of evident muscular artefacts. We remind the reader to[START_REF] Corsi | Integrating eeg and meg signals to improve motor imagery classification in brain-computer interface[END_REF] for a detailed description of the experiments.

Figure 3 . 1 :

 31 Figure 3.1: Connectivity estimators as a function of normalized distance between EEG electrodes. Coherence values are in orange while imaginary coherence results are in green. Connectivity estimated are extracted from all the available experimental observations in the β frequency band. [30] ©2018 IEEE

Figure 3 . 2 :

 32 Figure 3.2: Statistical contrast maps between motor imagery and resting states obtained with band power features. Results are shown for one-tailed permutationbased t-tests (p<0.01). In Panel A) the obtained t-values are illustrated for individual representative frequency bins within the theta band, B) alpha band, C) beta band, and gamma band. In the beta band, results also remained significant after correction for multiple comparison (p<0.05, FDR corrected).

Figure 3 . 3 :

 33 Figure 3.3: Statistical contrast maps between motor imagery and resting states in the beta band. In panel A, we report the results in beta band for connectivity and node strength values. In panel B, we show the results of node strength for individual representative frequency bins within the beta band. Results for coherencebased features are in top line, while those for imagery coherence-based features are in bottom line. Only the twenty most discriminant values are represented here for illustrative purposes.

  haviour of C and IC estimators. C is computed from the cross-spectrum of the two signals and it is sensitive to the signals' amplitude synchronization, i.e. when signals oscillate (or vary) at the same frequency. IC is also sensitive to the phase synchronization measuring possible time shifts between the signals.To study these behaviours, we consider an example of two equal sine waves oscillating at 10Hz, but shifted one with respect to the other within the [-π/2, π/2] interval. Fig.3.4 shows that while C value is constant along the entire phase shift range, IC varies in a way that it tends to zero when the two signals are perfectly in phase (i.e. ∆ = 0). In an additional analysis, re, we report in Appendix C an additional analysis to demonstrate that imaginary coherence between those signals can be analytically written as a function of their relative time delay.

Figure 3 . 4 :

 34 Figure 3.4: Relationships of coherence/imaginary coherence with phase difference. In Panel A) coherence is in pink and imaginary coherence in green, showing the functional connectivity between two sine waves at 10 Hz as function of their temporal shift. The shift, corresponding here to a phase difference, varies from 0 to π in steps of π /500. At each shift value, the two connectivity estimators are evaluated. Panel B) shows the sine waves with different phase differences. In panel 1), a positive ∆ of π /4 in panel 2), a negative ∆ of -π/4.

Figure 3 . 5 :

 35 Figure 3.5: Phase difference properties and discrimination ability. Panel A), results of permutation-based t-tests in the beta band across all subjects are shown for brain networks reconstructed from the phase difference between EEG signals. Panel B) results of permutation-based t-tests obtained with node strength values extracted from the previous brain networks. Panel C), Spearman correlation plot between imaginary coherence and phase difference values considering all the connections including C3 electrode for one representative subject.

Figure 3 . 6 :

 36 Figure 3.6: Improvement of classification performance. Bar plots show the percentage of relative increment between the best combination of features (i.e., coherence-based node strength S C , imaginary coherence-based node strength S IC and power spectrum P. The pie diagram in the inset illustrates the percentage of times that a specific combination of features has been selected across subjects.Two different cases are considered, distiguished by a different color code. In blue, we report results of accuracy improvements when selected electrodes are located in controlateral sensorimotor area (Contra), i.e. FC5, FC3, FC1, C5, C3, C1, CP5, CP3, CP1. In red, results relate to classification analysis when electrodes are in both contralateral and ipsilateral sensorimotor areas (Contra+Ipsi), ie. FC5, FC3, FC1, C5, C3, C1, CP5, CP3, CP1, FCZ, CZ, CPZ, FC2, FC4, FC6, C2, C4, C6, CP2, CP4, CP6

Figure 3 . 7 :

 37 Figure 3.7: Brain features selected by the classification procedure. The color codes for the group-averaged number of times that a specific feature -in the electrodefrequency space -has been chosen during the sequential feature selection algorithm. The results for P features are illustrated in the top line, those for S C in the middle line and those for S IC in the bottom line. Results relates to the classification analysis framework when only electrodes in the contralateral sensorimotor areas are selected.

  mechanisms, consisting in a simultaneous amplitude and phase synchronization of large-scale brain activity. Taken together, our results provide new insight into the network behaviours sub-serving brain functional changes during MI, and open new perspectives to improve BCI systems.

Figure 3 . 9 :Figure 3 . 10 :Figure 3 . 12 :

 39310312 Figure 3.9: Independent component obtain with Infomax algorithm. We show the first 10 ICs and the first 60s of recording. For each component, both the signals and the associated topography are represented to facilitate the preprocessing procedure.

Figure 3 . 13 :

 313 Figure 3.13: FBCSP procedure. EEG signals recorded from N channels are bandpassed in frequency intervals. Then, for each frequency interval, the common spatial pattern is computed. Finally, the posterior probability is estimated in order to weight CSP features to perform a classification fusion to integrate different type of features.

Figure 3 . 14 :

 314 Figure 3.14: Scheme of the fusion procedure. Each type of feature is computed and used to train a classifier. Then, the posterior probability is computed in each case and a Bayesian approach is used to integrate the different features.

Figure 3 . 15 :

 315 Figure 3.15: Improvement of classification performance. Bar plots show the percentage of relative increment between the best combination of features (i.e., coherence-based node strength S C , imaginary coherence-based node strength S IC and common spatial pattern CSP. The pie diagram in the inset illustrates the percentage of times that a specific combination of features has been selected across subjects.

44

 44

( 4 . 1 )

 41 , that becomes here:L = D -Â (4.3)

Chapter 4 .

 4 Improving graph connectivity estimation with graph signal processing 49 lated observations in the H 1 hypothesis. A graphical example of the effect of the transform T, is presented in Fig. 4.1 for the case of bidimensional Gaussian data whose mean and covariance matrix differ under the H 0 , H 1 . The original data are represented in Fig. 4.1 a), whereas their transformed versions are plotted in Fig. 4.1 b). The transformed data are unitary variance, zero-centered under H 0 and they are uncorrelated under H 1 .

Figure 4 . 1 :

 41 Figure 4.1: Example of transformation effect. In a) we have bidimensional Gaussian distributions which differ under mean and covariance matrix . In b) we have data in the transformed domain

50 4. 5 .

 505 Jensen divergence of connectivity statesbeing R(x) the Log-Likelihood Ratio * :

(4. 13 )

 13 By adding and subtracting the term|η 2 n |/[σ 2 n (σ 2 n -1)] we rewrite the equation (4.13) as:

Figure 4 . 2 :

 42 Figure 4.2: J-Divergence contributions as function of mean η and standard deviation σ: a) J (σ,η) for variables whose conditional standard deviation differ under H 1 and H 0 , and b) J (η) for variables with invariant conditional standard deviation.

4 . 3 , 53 Figure 4 . 3 :

 435343 Figure 4.3: Graphic interpretation of the score computed for the first element in the vector l

2 :

 2 Compute the denoised Laplacian L by a: Selecting the number N L , of smallest eigenvalues and the number N H of largest eigenvalues to retain b: Computing L

Chapter 4 . 55 Algorithm 4 . 2 J

 45542 Improving graph connectivity estimation with graph signal processing -divergence and score computationInput: Conditional means µ 0 , µ 1 and covariance matrices K 0 , K 1 of l = Vec( l) under H 1 and H 0

, 3 :

 3 Step 3: Score computation a: Compute S n , n = 0, • • • N -1 as

Fig. 4 .

 4 Fig. 4.4 represents the estimated adjacency matrices under the two conditions H 1 and H 0 in presence of perturbations. Under H 1 some strong connections are present and their value is affected by the perturbations. Under H 0 (4.4 B) there are no evident link and  fluctuates around zero because of the perturbations.

Figure 4 . 4 :

 44 Figure 4.4: Adjacency matrix with synthetic data. Âij is represented under H 1 in A) and under H 0 in panel B).

  .5 the first and the 10th eigenvectors on graph under H 1 hypothesis. The first eigenvector, in Fig. 4.5 A) is completely smooth on the graph and within a subset of connected nodes. Fig. 4.5 B) represents the 10th eigenvector on graph. Its variations over the graph highlight another community, but it has higher variability over connected node compared to the first eigenvector.

Figure 4 . 5 :

 45 Figure 4.5: Example of eigenvectors on graph. They relate to the Laplacian matrix estimated from synthetic data under H 1

Figure 4 . 6 :

 46 Figure 4.6: Results of Frobenius distance on synthetic data. Several perturbation configurations are represented: in panel a)σ w =0 and σ b =0,in panel b) σ w =1.2 and σ b =0,in panel c) σ w =0 and σ b =2 and in panel d) σ w =1.2 and σ b =2. In the different colors (in the legend) we represent the different subspaces.

59 Figure 4 . 7 :

 5947 Figure 4.7: Results of J-divergence analysis on synthetic data. Several perturbation configurations are represented: in panel a) σ w =0 and σ b =0, in panel b) σ w =1.2 and σ b =0, in panel c) σ w =1.2 and σ b =2. In the different colors (shown in the legend) we represent the different subspaces for the filtering.

60 4. 7 .

 607 Real BCI measurements

Chapter 4 . 61 Figure 4 . 8 :

 46148 Figure 4.8: Results of J-divergence analysis for real EEG data. We report the cumulative J-divergence CJ n as function of the first 20 variables. In the different colors (shown in the legend) we represent the different subspaces used to filter Laplacian.

Fig. 4 .

 4 [START_REF] Baillet | Electromagnetic brain mapping[END_REF] shows associated results by plotting the J-divergence, averaged across subject, as a function of the time window index m, m = 0 • • • M -1. Our findings show that for almost all the time intervals m, i.e. on 7 intervals out of 9, the denoising algorithm leads to higher J-divergence compared to the absence of denoising. This result is really interesting because it shows that, even with short time-interval, our method succeeds in separating the two cognitive states.

Figure 4 . 10 :

 410 Figure 4.10: Results of J-divergence analysis over a moving window on real data.We plot the J-divergence over M = 9, 1s long, time intervals with 50% overlapping, versus the time interval index. The J-divergence is computed in β band and averaged across subjects.

  ) n = 2η 2 n , n = P, • • • N -Brain signals are time-varying and non-stationary, which means that the dynamic of the signals and their moments change over time. The development of a robust tool to estimate graph functional connectivity is crucial for BCI applications. L1-Principal Component Analysis (L1-PCA) is a state-of-theart computational tool able to identify relevant components in data with noise and outliers. Significant efforts have been made in the direction of adapting the PCA to time-varying data. Here, we analyse a layered version of L1-PCA, which we refer to as Deep L1-PCA. Deep L1-PCA is obtained by recursive application of two steps: estimation of L1-PCA basis and extraction of the first rank projector. Deep L1-PCA is robust to outliers and it produces relevant components with a reduced computational cost. Here, we show how the Deep L1-PCA can be applied in a classification scheme of time-varying data. This framework is preliminarily applied to simulated graphs and then to EEG data recorded during motor imagery-based BCI experiments. The proposed approach shows potential in the application of real BCI systems for short-time connectivity estimation and discrimination of brain states. This chapter contains part of the work from Orru, G., Cattai, T., Colonnese, S., Scarano, G., De Vico Fallani, F., Markopoulos, P., and Pados, D. (2019, September). Deep L1-PCA of Time-Variant Data with Application to Brain Connectivity Measurements. In 2019 27th European Signal Processing Conference (EUSIPCO) (pp. 1 -5). IEEE. [131] Brain signals are time-varying and non-stationary: the dynamic of the signals and

5 .

 5 Before detailing the mathematical formulation, a presentation of the key idea of the Deep L1-PCA is needed. At the first layer, the original data are partitioned, and for each partition the L1-Principal Component (L1-PC) is estimated. The principal components (L1-PCs) of all partitions are then organized to build a new group, which becomes the input of the second layer. The procedure at the second layer naturally follows: the input data are partitioned in subgroups and the L1-PCs extracted for each subgroup become the input for the next layer. This iterative Chapter Framework for short-time graph connectivity at layer k W K batch size at layer k R k (r k ) L1 -PC matrix (or vector) at layer k B k (b k ) optimal binary matrix (or vector) at layer k  estimated adjacency matrix) L estimated Laplacian matrix) D estimated degree matrix)

L1-Figure 5 . 1 : 5 . 2 .

 5152 Figure 5.1: General Deep L1-PCA computation structure.

Figure 5 . 2 :

 52 Figure 5.2: A) Time-variant bidimensional normal data and in B) the underlying probability density function (PDF) orientation.

Figure 5 . 3 :

 53 Figure 5.3: Time-variant data and first order L1-PCA approximation at different layers : layer 1 in A), layer 2 in B) and layer 3 in C).

  m = D(1)REF (m), where D(1) (m) is the degree matrix computed from Â(1) (m) and D(1)REF (m) is obtained from A

Figure 5 . 4 :

 54 Figure 5.4: Graphical representation of Deep L1-PCA computation with an M = 8 original data Ξ. At each layer k (k = 1, 2, 3) the L1-PCAs r k (m) are computed and they become the input in R k+1 (m) for the next layer.

Figure 5 . 5 :

 55 Figure 5.5: Mean squared error map computed on simulated data between the degree matrix corrupted by noise with σ w = 0.6 and its reference version (i.e. without perturbations).

( 1 )

 1 m = D(1) (m), with m = 1, .., M, while R (0) k is obtained when ξ (0) m = D(0) (m), with m = 1, .., M. For this analysis, we remind that M = 8.

( 1 )

 1 m = L(1) (m) and ξ (0) m = L(0) (m) for m = 1, .., M .

Figure 5 . 6 :

 56 Figure 5.6: Accuracy of classification framework for synthetic data for each layer and each batch. In the rows we have the different layers and in the columns we have the batches, which are 4 at the first layer, 2 at the second and 1 batch at the final layer. In panel A), we report accuracy results when only the Laplacian diagonal elements are taken as classification features. In panel B) results are obtained when all the Laplacian elements are considered.

Fig. 5 .

 5 Fig. 5.6 B) presents associated classification results. The first consideration is that, as for Fig. 5.6 A), accuracy generally increases considering higher layers (from k = 1 to k = 3). Interestingly, batches associated to highest classification performances in Fig. 5.6 B) are the same of Fig. 5.6 A), meaning that Laplacian extra-diagonal elements do not add further information about data. Those elements fruitfully contribute to the classification since the accuracies are generally higher when all the Laplacian is considered. For comparison's sake, we perform analogous classification tests for the 8 batches containing the original data (i.e. when no principal component analysis is applied). Results are presented in Table 5.2. The first row shows results obtained when feature vectors are built with degree matrices while the values on the second row are obtained using all the Laplacian matrices. Accuracy values obtained on the batches of original data are generally lower compared to those collected after Deep L1-PC analysis and, even in this case, Laplacian extra-diagonal elements generally improve the classification.

Fig. 5 . 7 A

 57 Fig. 5.7 A) results are related to the case in which the feature vector is composed by only the Laplacian diagonal elements (i.e.node strength values). Fig. 5.7 B)

  3 results for the 8 observations on original data. Accuracy values are generally lower compared to those obtained after Deep L1-PC computations. Even for original data collected in short-time windows the inclusion of extra-diagonal elements helps the classification.We demonstrate the potential of Deep L1-PCA to improve the discrimination of two mental states during motor imagery task. The most important contribution concerns the possibility to consider short-time intervals to obtain a robust connectivity estimation. This aspect is fundamental for BCI applicability in real scenario.

Figure 5 . 7 :

 57 Figure 5.7: Accuracy of classification framework for real EEG data associated to each layer and each batch. In the rows we have the different layers and in the columns we have the batches, which are 4 at the first layer, 2 at the second and 1 batch at the final layer. In panel A), we report accuracy results when only the Laplacian diagonal elements are taken as classification features. In panel B) results are obtained when all the Laplacian elements are considered as features.

Fig. 5 .

 5 Fig.5.8 we have the M time intervals and we indicate with red arrows the time interval in which we have enough data to pass to the higher layers.

Figure 5 . 8 :

 58 Figure 5.8: Graphical representation of data recordings. As the time passes, we collect m time-windows in which we can estimate the graph Laplacian. As soon as new time intervals are recorded we can compute the Deep L1-PCs at several layers. Here, we indicate with red arrows the time instant (and the associated time interval) in which we can reach an higher level. In fact, after 500ms, we have the first Laplacian estimation and the first classification can be performed. After the second interval (i.e. at 1s) the first L1-PCs can be computed. We have to wait until 2s (i.e. the 4th interval) to collect all the samples to complete the first layer with k = 1 and to perform the second L1-PCA, reaching the second layer. At 4s, we have all the samples needed to compute the last L1-PCA.

  L1-PCA. This is an iterative procedure, where in the first layer, original data are partitioned in batches and L1-PCs are extracted in each group. In the next layer, the input data are the L1-PCA bases obtained as output of the previous layer. With the recursive application of this method, original data are represented in a more compact way. The method ends with the possibility to extract one global PC. The Deep L1-PCA is applied to graph synthetic data to verify the robustness to noise and the applicability to classification scenarios. The last part of this work concerns results on real EEG data recorded during motor imagery-based BCI experiments.

  time interval and the whole Deep structure containing the computations at each layer. The proposed method does not model the time evolution of the graph, which can in principle offer new insight in brain (re)organization processes. In this direction, we proposed in [44], a compound Markov random field model, able to jointly model signal values, graph connectivity and the community structure. Our method, based on signal, connectivity and community fields, generalizes classical Markov model to simultaneously address problems of learning, signal recovery and community detection. As future work we can further develop the model to describe time-varying SoG evolution. One possible strategy to adapt the model is to add a process to consider the temporal evolution of graph edges. In principle, the modified version of the Markov random field model has the potential to describe the temporal evolution of the graph and it might be applied in BCI to observe spatio-temporal changes during a motor or cognitive task.

Table 3 . 1 :

 31 Table of main notation in Chapter 3.

		Description
	C	spectral coherence estimate
	IC	imaginary coherence estimate
	∆	phase difference
	N	total number of nodes
	S C	coherence-based node strength
	S IC	imaginary coherence-based node strength
	S IC	phase difference-based node strength

Improving graph connectivity estimation with graph signal processing 41 4.1 Introduction Functional

  connectivity (FC) describes the interactions between brain areas[START_REF] Bastos | A tutorial review of functional connectivity analysis methods and their interpretational pitfalls[END_REF]. It can be modelled as a graph, which represent one possible formalism to represent networked systems[START_REF] Newman | Networks: an introduction[END_REF][START_REF] Torres | The why, how, and when of representations for complex systems[END_REF]. It has been recently demonstrated that graph statistics, such as node strength, efficiency and modularity[START_REF] Gonzalez-Astudillo | Network-based brain computer interfaces: principles and applications[END_REF]. The detection of brain FC-based features can boost several applications, such as brain-computer interfaces (BCIs). BCIs are communication systems enabling a subject to interact with external work without neuro-muscular activity[START_REF] Wolpaw | Brain-computer interfaces: Definitions and principles[END_REF][START_REF] Wolpaw | The wadsworth center brain-computer interface (bci) research and development program[END_REF]. A critical requirement for BCI to work is the correct detection of the user's intent. Although the research significantly advanced on BCIs, there is still a limitation, commonly known as BCI inefficiency[START_REF] Thompson | Critiquing the concept of bci illiteracy[END_REF]. It indicates the fact that there is a non-neglecting percentage of users who cannot be trained to use the interface and cannot correctly use it. This limitation motivated us to propose novel techniques to increase the detectability of the user's cognitive states. Our intent is to develop original tools to have a more

	Then, we perform analyses with real EEG data recorded during motor
	imagery-based BCI experiments. The novel formulation of the J-divergence
	enables to simultaneously quantify the distance between FC networks
	in motor imagery and resting, as well as to highlight the contribution
	of each Laplacian variable to total J-divergence. Experimental results on
	real data demonstrate the potential of the denoising algorithm on short
	time-windows. This approach provides new practical tools to robustly
	estimate FC networks and it opens new possibilities in the implementation

of real-time BCI systems. In this chapter are reported figures and part of the text from a recently submitted paper: Cattai, T., Scarano, G., Corsi, M. C., Bassett, D. S., De Vico Fallani, F., & Colonnese, S. (2020). Improving J-divergence of brain connectivity states by graph Laplacian denoising. arXiv preprint arXiv:2012.11240. [33] Chapter 4.

robust brain FC features to better separate two cognitive states. The practical implementation of FC estimation from signals acquired at graph vertices (e.g. EEG electrodes) is a not an easy task because of the noise, the high number of edges, artefacts (e.g. ocular, cardiac) and non-stationarity of brain signals.

Table 4 .

 4

	Chapter 4. Improving graph connectivity estimation with graph signal
	processing	43
	Notation	Description
	A , Â	adjacency matrix (real, estimated)
	V	set of all nodes
	N	total number of nodes
	E	set of all links
	D, D	degree matrix (real, estimated)
	L, L	Laplacian matrix (real, estimated)
	λ, λ	eigenvalue (real, estimated)
	u, û	eigenvector (real, estimated)
	U L , U M , U H	subset of smallest, central, larger eigenvalues
	L, l	filtered graph laplacian matrix and vector
	T	transformation matrix
	x	vectorized laplacian in the transformed domain
	J	J-divergence
	S	score

authors propose a new technique to ef-1: Table of main notation in Chapter 4f. ficiently estimate the adjacency matrix by creating and modifying embeddings related to each graph vertex. Since the estimation of FC requires a lot of time and computational resources, it is possible to cluster FC in communities of synchronous components. One example is the method introduced i n[63]. It consists in the application of k-means clustering followed by a tensor decomposition in order to reduce the FC data. Classical signal processing operations have been generalized into the graph setting, where signals are associated to graph nodes, giving rise to the research domain of graph signal processing (GSP)

Table 5 . 1 :

 51 Table of main notation in Chapter 5.

3.5. ConclusionsTable3.4: Average accuracy across cross-validation is reported for each subject and each combination of feature when CSP is used.

* All the power spectral estimates are computed with Welch method, with 1s length Hanning windows and overlap of 50%.

* For each subspace configuration, we compute the Frobenius distance F between its noisy and GT versions. See Definition 2 in[START_REF] Baksalary | On subspace distances determined by the frobenius norm[END_REF] for the mathematical formulation.

Contents

Appendix A: Example of Preprocessing

Here, we show an example of the preprocessing procedure on one subject. In Fig. 3.8, we report the first 60s recorded for EEG signals. EEG channels are ordered from frontal to parietal, and we recognize blink effects in the data recorded from the first channels as regular spikes. We report in Fig. 3.10 EEG signals after ICA. Comparing EEG data with those in Fig. 3.8, we remark the significant decrease of high and regular spikes in the first components, caused by the elimination of the blink component.

Appendix B: Brain connectivity changes during MI tasks in all frequency bands

We report the results of permutation t-test in all the frequency bands of interest, i.e. theta, al pha, beta and gamma, in order to compare brain mechanisms in all the frequencies. We consider the case of connectivity changes in Fig. 3.11 while node strength behaviours are in Fig. 3.12 Figure 3.11: Statistical contrast maps between motor imagery and resting states for connectivity features. In top line, we have results for coherence features and in the bottom line for imaginary coherence ones. In Panel A) results for coherence in theta band, in B) for imaginary coherence in theta band, in C) for coherence in al pha band, in D) for imaginary coherence in al pha band, in E) for coherence in beta band and in F) for imaginary coherence in beta band in G) for coherence in gamma band and in H) for imaginary coherence in gamma band . Only the twenty most discriminant connections are represented here for the sake of simplicity.

Real BCI measurements

Hz) or β (14-29 Hz) band, depending on the subject and the task [START_REF] Meirovitch | Alpha and beta band event-related desynchronization reflects kinematic regularities[END_REF].

Here, we perform the analysis in β frequency band, but in a training BCI scenario, the frequency band, or eventually the frequency bin, can be tuned according to the subject response. For this reason, we filter the connectivity matrix in β frequency band band:

After demonstrating that the denoising based on U L∪H subspace allows a better separability of connectivity states under H 1 and H 0 , we limit the analysis to U L∪H and U ALL .

We compute the scores as explained in section 4.5A) and we report the associated results in Fig. 4.9. In the first row, we collect the results referring to extra The first interesting result is that, even in the case of connectivity matrix restricted to β band, the proposed denoising improves the J-divergence of the two connectivity states (from 79.47 to 160.77).

Considering results from score analysis, different observations are in order.

In Fig. We focus on the score related to nodes weights (Laplacian diagonal elements ) in Fig. 4.9(c-d). Firstly, the score range in absence of denoising is smaller compared to when the Laplacian denoising is applied, i.e. maximum values are 0.4 and 0.7.

Then, the scores obtained without denoising are larger on nodes located in frontal, temporal or parietal area, such as FP Z and P 4 . After denoising, the score values are higher on sensory-motor areas, and we can recognize some more relevant nodes, such as C 2 and FC 5 .

Thereby, the proposed score procedure together with the Laplacian denoising provides an original tool for the analysis and the interpretation of the brain cognitive states.

Short-time estimation of Laplacian coefficients in β band

BCI systems aim to realize a real time communication between the user and the interface [START_REF] Wolpaw | Brain-computer interfaces: principles and practice[END_REF][START_REF] Wolpaw | The wadsworth center brain-computer interface (bci) research and development program[END_REF]; thereby, decreasing the observation time T oss for Laplacian computation is beneficial for potential applicability to online motor-imagery BCI.

In this direction, we test the Laplacian denoising when the window length of the observed signals is reduced to T oss = 1s. For this analysis, we consider a moving window of length T oss = 1s and we shift it by m∆t, m = 0 • • • M -1, with M = 9 ∆t = 0.5s, in order to analyse the total available length of 5s over nine 50% overlapping temporal windows. [START_REF] Shenoy | Towards adaptive classification for bci[END_REF].

In order to obtain a global information across the subjects, we take into consideration the 20 subjects of the experimental study. For each subject, we compute the spectral coherence on the m-th temporal interval, m = 0 • • • M -1 as in Eq.

(5.11) both for resting (H 1 ) and motor imagery (H 0 ) state. Then we compute the

Scientific production

Submitted journal papers Peer-reviewed conference papers

Contributions to conferences

• Oral presentation at Signal processing over graphs, with a focus on neuroscience data workshop, September 2019