Prof Pascal Benquet

Prof Mohamad Khalil

Prof Mark Woolrich

Prof. Fadi Karame

Prof Prof Catherine Julie Marque

Peron

Dr Joan Duprez

Judie Tabbal

Aya Kabbara

Mohamad Khalil

Pascal Benquet

Mahmoud Hassan

Maxime Yochum

Dynamics of task-related electrophysiological networks: a benchmarking study

Keywords: Magneto-encephalography (MEG) Electrophysiological brain networks Dynamic functional connectivity Dimensionality reduction Source separation Electroencephalography, dynamic functional connectivity, brain networks states, neural mass models

I would like to express my

List of Figures

 (A and B)are connected to a pool of neurons C. In this example, pool A communicates with C while functional connectivity between B and C is suppressed. b. In CTC, the phase relationship determines the routing: the neurons in A and C oscillate in phase, whereas the neurons in B do not oscillate in phase with C. The waveforms represent oscillatory population activity in the three regions, whereas the small vertical lines represent spiking activity. c. The magnitude of the pulses of inhibition controls the routing: the flow of information is controlled by an increase of alpha-band oscillations which inhibits firing in pool B, and a decrease in alpha oscillations of neurons in A and C, allowing communication by release from inhibition. This figure is extracted from [START_REF] Bonnefond | Communication between Brain Areas Based on Nested Oscillations[END_REF]

Résumé en français

En tant que système complexe, le cerveau traite de manière flexible les informations grâce à une reconfiguration dynamique des réseaux neuronaux sur une échelle de temps de l'ordre de la milliseconde. Un objectif majeur en neurosciences est de décrire l'organisation spatio-temporelle du cerveau comme une série d'«états de connectivité fonctionnelle » transitoires à travers une analyse dynamique des réseaux. Ce domaine prend de l'ampleur car il permet non seulement d'aborder les processus cognitifs, mais aussi d'apporter des informations importantes sur les altérations fonctionnelles des principaux motifs de connectivité dans le cadre des pathologies neurologiques. Dans ce contexte, deux enjeux principaux ont été identifiés : (1) A quel point les techniques de sujets sains) au cours de la tâche de conflit cognitif nommée Simon-Task. Une variante de l'analyse des composantes indépendantes a été utilisée pour dériver des composantes statistiquement indépendantes dans les deux groupes. Les résultats démontrent l'existence de différences spatiotemporelles dans les états dynamiques des réseaux cérébraux entre les sujets sains et les patients.

Track dynamic brain network states using EEG/MEG

Here the objective is to explore the appropriate methodology that allows extracting relevant connectivity patterns, underlying neural activity when performing tasks. First, three independent MEG datasets from 95 healthy subjects were used during motor and working memory tasks operating on variable time scales. We used the "EEG/MEG source connectivity" method to estimate dynamic functional connectivity (dFC) matrices at the cortical level. Then, several data-driven decomposition techniques were applied to reduce dFC dimensionality by deriving principal brain patterns with their temporal activation. The performance of these techniques was evaluated and compared at group and subject levels.

Second, the previous pipeline was tested using a physiologically based ground truth computational model of a human brain to simulate HD-EEG activity during cognitive task driven at a rapid time scale, as a way to assess a quantitative evaluation of decomposition methods along with multiple key factors used in the pipeline. Primarily, both qualitative and quantitative results show promising outcomes of tested methods with some variability in terms of spatial and temporal accuracy, related to task complexity and time scale. Thus, our findings suggest a careful choice of these methods as they may influence results interpretation.

Tracking dysfunctional electrophysiological networks in Parkinson's disease

The main purpose of this work was to identify the major alterations evoked in the extracted dynamic network states for PD patients. For this reason, HD-EEG data was recorded from 31 subjects (21 patients, 10 healthy subjects) during a Simon task. A variant of temporal independent component analysis was used to derive statistically independent components for both groups. Results demonstrate a difference in the spatiotemporal behavior of the dynamic network states between healthy subjects and PD patients.

Chapter 1. GENERAL INTRODUCTION 1.1. Motivation

Primarily, brain imaging techniques have focused on localizing the brain activity at specific brain regions to study neural information processing. However, there is emerging evidence that most behavioral and cognitive mechanisms involve various distributed brain areas working in concert (Menon, 2015). This has shifted neuroscientists' interest toward a deeper understanding of 'largescale' functional organization and its impact on fundamental brain functions [START_REF] Bressler | Large-scale brain networks in cognition: emerging methods and principles[END_REF][START_REF] Dosenbach | A dual-networks architecture of top-down control[END_REF][START_REF] Fox | Spontaneous fluctuations in brain activity observed with functional magnetic resonance imaging[END_REF][START_REF] Greicius | Functional connectivity in the resting brain: A network analysis of the default mode hypothesis[END_REF]. Consequently, the human brain is currently modeled as a complex, large-scale network of functional connections (edges) between distributed brain regions (nodes).

Two prominent features of this network increasingly attract researchers in the neuroscience field. First, the brain network can flexibly and dynamically reconfigure at a sub-second time scale to maintain a real time-efficient control of cognitive activity (Allen et al., 2018(Allen et al., , 2014;;Bassett et al., 2011b;[START_REF] Calhoun | The Chronnectome: Time-Varying Connectivity Networks as the Next Frontier in fMRI Data Discovery[END_REF][START_REF] Chang | Time-frequency dynamics of resting-state brain connectivity measured with fMRI[END_REF]Hassan et al., 2015;[START_REF] Hutchison | Dynamic functional connectivity: promise, issues, and interpretations[END_REF][START_REF] Preti | The dynamic functional connectome: State-of-the-art and perspectives[END_REF]Vidaurre et al., 2018;[START_REF] Zalesky | Time-resolved resting-state brain networks[END_REF]. Second, brain network dynamics reveal recurring connectivity patterns, namely 'dynamic brain network states', which transiently span over time. Hence, the brain can occupy any of several states over time, each with a specific network topology (Allen et al., 2014;[START_REF] Calhoun | Time-Varying Brain Connectivity in fMRI Data: Whole-brain data-driven approaches for capturing and characterizing dynamic states[END_REF]Iraji et al., 2020;[START_REF] Leonardi | Disentangling dynamic networks: Separated and joint expressions of functional connectivity patterns in time[END_REF]Miller et al., 2016;O'Neill et al., 2018O'Neill et al., , 2017)).

Accordingly, high temporal resolution techniques along with accurate time-resolved tools are required to explore the spatiotemporal dynamics of key connectivity patterns that underlie neural activity.

Although considerable effort has been recently dedicated to track fast temporally evolving state sequences, an authentic evaluation of different existing analytical approaches remains missing.

On the other hand, it is now recognized that several neurological disorders, such as Parkinson, Epilepsy, Alzheimer are associated with disruptions in functional brain networks (Fornito et al., 2015a;[START_REF] Fornito | Connectomics: a new paradigm for understanding brain disease[END_REF][START_REF] Helm | Neuronal connectivity in major depressive disorder: a systematic review[END_REF][START_REF] Kaiser | Large-Scale Network Dysfunction in Major Depressive Disorder: A Meta-analysis of Resting-State Functional Connectivity[END_REF][START_REF] Mulders | Resting-state functional connectivity in major depressive disorder: A review[END_REF]. Still, little is known about how these pathologies affect brain network states (BNS) architecture and influence their dynamics.

Overview

This manuscript is organized as follows:

 Chapter 2 covers the background of the work performed in this thesis. In particular, we define the brain connectivity types, the non-invasive neuroimaging techniques with their relative spatiotemporal resolution. Afterward, we outline the 'EEG/MEG source connectivity' approach used to estimate the functional connectivity, and we focus on its dynamic assessment. Then, the methods used to reduce dimensionality and to derive dominant brain network states are described. Finally, we introduce the main objective of this thesis.

 In Chapter 3, we present the thesis results in the form of published, submitted, and in preparation articles. For each article, we provide the corresponding objectives, methodology, results, and discussion.

 Finally, general conclusions are summarized in Chapter 4. An overall discussion with several methodological considerations is given. We also suggest some ideas for future studies.

Chapter 2. BACKGROUND AND OBJECTIVES

In a very networked system like the brain, it is necessary to shift from thinking in terms of isolated brain regions, and adopt instead the language of networks (Bullmore and Sporns, 2009a;[START_REF] Sporns | Networks of the Brain[END_REF]. In that view, the network itself is the unit, not the brain area. Processes that support behavior depend on the interaction of multiple brain areas, which are dynamically recruited into multi-region assemblies. Pathological perturbations of the brain are also rarely limited to a single region. Local dysfunctions often affect other regions, resulting in large-scale network alterations (Fornito et al., 2015b;[START_REF] Stam | Modern network science of neurological disorders[END_REF]. Here we describe the notions of brain connectivity, the available neuroimaging techniques, the crucial dynamic behavior of the brain networks and their potential applications in cognitive and clinical neuroscience.

Brain connectivity

When applied to the brain, the term 'connectivity' refers to several aspects of brain organization that ensure optimal information processing [START_REF] Horwitz | The elusive concept of brain connectivity[END_REF]. In this context, brain connectivity can be defined as a set of physical/structural links ('structural connectivity') [START_REF] Sporns | The Human Connectome: A Structural Description of the Human Brain[END_REF], statistical interactions ('functional connectivity'), or causal relationships ('effective connectivity') (Friston, 1994) between different neural elements. Hereafter, we will describe the main characteristics of these three distinct forms of connectivity, (see Figure 1). Functional connectivity reflects the statistical coupling between brain areas while effective connectivity accounts for directionality by measuring the causal influence of one brain region over another.

Structural Connectivity (SC)

Structural Connectivity (SC), also called Anatomical Connectivity (AC), refers to the physical (synaptic) connections between fiber tracks (white matter, myelinated axons) linking spatially distant brain regions, giving rise to the so-called connectome. It can be mapped using a modern neuroimaging technique, such as the Diffusion Tensor Imaging (DTI) [START_REF] Hagmann | Mapping the Structural Core of Human Cerebral Cortex[END_REF][START_REF] Hagmann | Mapping Human Whole-Brain Structural Networks with Diffusion MRI[END_REF].

Based on this technique, the locations of anatomical highways spanning the human brain cortex have been estimated and collected to form networked large-scale pathways [START_REF] Basser | In vivo fiber tractography using DT-MRI data[END_REF][START_REF] Behrens | Non-invasive mapping of connections between human thalamus and cortex using diffusion imaging[END_REF][START_REF] Hagmann | DTI mapping of human brain connectivity: statistical fibre tracking and virtual dissection[END_REF]. This network architecture has a precise internal structure (Bassett et al., 2011a), which remains relatively stable across individuals, task and context [START_REF] Finn | Functional connectome fingerprinting: identifying individuals using patterns of brain connectivity[END_REF][START_REF] Kaufmann | Delayed stabilization and individualization in connectome development are related to psychiatric disorders[END_REF]. One of the largest open-access datasets in neuroimaging is the human connectome project (HCP) offering a large amount of structural (and functional) connectivity matrices in health and disease (Van Essen et al., 2013). The structural HCP data was used here in the study 2 to upgrade the human computational model.

Even if the main large scale circuits can be found in every individual, inter-subject variability exists in the size/locations of brain anatomical connections [START_REF] Sporns | The Human Connectome: A Structural Description of the Human Brain[END_REF][START_REF] Thompson | Three-dimensional statistical analysis of sulcal variability in the human brain[END_REF], [START_REF] Llera | Inter-individual differences in human brain structure and morphology link to variation in demographics and behavior[END_REF][START_REF] Müller | Intersubject variability in the analysis of diffusion tensor images at the group level: fractional anisotropy mapping and fiber tracking techniques[END_REF]. Moreover, the connectome can be modified with age [START_REF] Dennis | Development of Brain Structural Connectivity between Ages 12 and 30: A 4-Tesla Diffusion Imaging Study in 439 Adolescents and Adults[END_REF][START_REF] Hagmann | White matter maturation reshapes structural connectivity in the late developing human brain[END_REF][START_REF] Lebel | Longitudinal Development of Human Brain Wiring Continues from Childhood into Adulthood[END_REF], pathology [START_REF] Dai | Disrupted structural and functional brain networks in Alzheimer's disease[END_REF][START_REF] Jeong | Functional and anatomical connectivity abnormalities in left inferior frontal gyrus in schizophrenia[END_REF][START_REF] Seeley | Neurodegenerative diseases target large-scale human brain networks[END_REF][START_REF] Zhang | White matter damage in frontotemporal dementia and Alzheimer's disease measured by diffusion MRI[END_REF][START_REF] Zhang | Altered functional-structural coupling of large-scale brain networks in idiopathic generalized epilepsy[END_REF] or with some mental processes such as learning [START_REF] Casey | Structural and functional brain development and its relation to cognitive development[END_REF][START_REF] Paus | Structural maturation of neural pathways in children and adolescents: in vivo study[END_REF][START_REF] Sampaio-Baptista | White Matter Plasticity in the Adult Brain[END_REF][START_REF] Schmithorst | Cognitive functions correlate with white matter architecture in a normal pediatric population: a diffusion tensor MRI study[END_REF][START_REF] Scholz | Training induces changes in whitematter architecture[END_REF].

Functional Connectivity (FC)

Background

Neuronal communication is a fundamental aspect of information processing and brain function. It has classically been determined by structural connectivity. Although anatomy supports the basic board for inter-areal signaling, it does not fully characterize communication. Even in the absence of changes in the structural connectivity, neuronal synchronization has been proposed to play crucial roles for routing the neural activity within and between brain areas. Dynamic changes in synchronization/oscillations can flexibly alter the pattern of interactions between brain regions on a fast time scale [START_REF] Varela | The brainweb: Phase synchronization and large-scale integration[END_REF]. Such flexible modulations are at the heart of perception and cognition [START_REF] Bosman | Attentional Stimulus Selection through Selective Synchronization between Monkey Visual Areas[END_REF][START_REF] Ruff | Simultaneous multi-area recordings suggest that attention improves performance by reshaping stimulus representations[END_REF][START_REF] Salazar | Content-Specific Fronto-Parietal Synchronization During Visual Working Memory[END_REF], and motor behavior (Arce- [START_REF] Arce-Mcshane | Primary motor and sensory cortical areas communicate via spatiotemporally coordinated networks at multiple frequencies[END_REF][START_REF] Wong | Coherent neuronal ensembles are rapidly recruited when making a look-reach decision[END_REF]. In this context, two essential frameworks have been suggested to assess the flexible routing of information based on oscillations [START_REF] Bonnefond | Communication between Brain Areas Based on Nested Oscillations[END_REF][START_REF] Kohn | Principles of Corticocortical Communication: Proposed Schemes and Design Considerations[END_REF] (see Figure 2).

i.

Communication Through Coherence (CTC) [START_REF] Bastos | Visual areas exert feedforward and feedback influences through distinct frequency channels[END_REF]Fries, 2015[START_REF] Fries | A mechanism for cognitive dynamics: neuronal communication through neuronal coherence[END_REF]. In the CTC framework, neuronal communication is established when the oscillations of the source/target areas in the neuronal pools are coherent and appropriately coordinated [START_REF] Fries | A mechanism for cognitive dynamics: neuronal communication through neuronal coherence[END_REF](Fries , 2009)). For instance, modulations in gamma oscillations and in their coherence in the thalamus and visual cortical areas were associated with attention tasks [START_REF] Gregoriou | High-Frequency, Long-Range Coupling Between Prefrontal and Visual Cortex During Attention[END_REF][START_REF] Zhou | Pulvinar-Cortex Interactions in Vision and Attention[END_REF], whereas coherence in other rhythmic oscillations was altered during higher task demands [START_REF] Rubino | Propagating waves mediate information transfer in the motor cortex[END_REF][START_REF] Wong | Coherent neuronal ensembles are rapidly recruited when making a look-reach decision[END_REF]. The reader can also refer to [START_REF] Buzsaki | Rhythms of the Brain[END_REF]Fries, 2015;Wang, 2010) for reviews.

ii.

Gating By Inhibition (GBI) [START_REF] Jensen | Shaping Functional Architecture by Oscillatory Alpha Activity: Gating by Inhibition[END_REF]. GBI framework suggests an active inhibition of the pathway that is not involved in the task, which is reflected mainly through alpha activity [START_REF] Foxe | The Role of Alpha-Band Brain Oscillations as a Sensory Suppression Mechanism during Selective Attention[END_REF][START_REF] Klimesch | EEG alpha oscillations: the inhibition-timing hypothesis[END_REF].

Figure 2. An illustration for the communication through coherence (CTC) (b) and gating through inhibition (GBI) (c) hypotheses. a. Two pools of neurons (A and B) are connected to a pool of neurons C. In this example, pool A communicates with C while functional connectivity between B and C is suppressed. b. In CTC, the phase relationship determines the routing:

the neurons in A and C oscillate in phase, whereas the neurons in B do not oscillate in phase with C. The waveforms represent oscillatory population activity in the three regions, whereas the small vertical lines represent spiking activity. c. The magnitude of the pulses of inhibition controls the routing: the flow of information is controlled by an increase of alphaband oscillations which inhibits firing in pool B, and a decrease in alpha oscillations of neurons in A and C, allowing communication by release from inhibition. This figure is extracted from [START_REF] Bonnefond | Communication between Brain Areas Based on Nested Oscillations[END_REF].

Metrics

Practically the FC is estimated as the statistical dependencies between neural assemblies. It can be estimated by measuring correlation, covariance, spectral coherence, or phase synchronization between pairs of time series [START_REF] Ansari-Asl | Quantitative evaluation of linear and nonlinear methods characterizing interdependencies between brain signals[END_REF][START_REF] Pereda | Nonlinear multivariate analysis of neurophysiological signals[END_REF][START_REF] Sakkalis | Review of advanced techniques for the estimation of brain connectivity measured with EEG/MEG[END_REF]. These time series can be measured using different neuroimaging techniques as electro-encephalography (EEG), magneto-encephalography (MEG), and functional magnetic resonance (fMRI). Unlike structural connectivity, functional connectivity is time-dependent and fluctuates on a very short time scale (milliseconds to seconds). This type of connectivity has proven its powerful ability to track brain activity during rest (Allen et al., 2018;Baker et al., 2014;[START_REF] De Luca | fMRI resting state networks define distinct modes of long-distance interactions in the human brain[END_REF]de Pasquale et al., 2010;[START_REF] Fox | Spontaneous fluctuations in brain activity observed with functional magnetic resonance imaging[END_REF][START_REF] Yuan | Spatiotemporal dynamics of the brain at rest--exploring EEG microstates as electrophysiological signatures of BOLD resting state networks[END_REF] and behavioral tasks [START_REF] Bassett | Learning-induced autonomy of sensorimotor systems[END_REF](Bassett et al., , 2011b;;Hassan et al., 2015;Kabbara et al., 2019;Mheich et al., 2015;O'Neill et al., 2017). Many researchers have also investigated functional networks alterations in neurological diseases [START_REF] Dai | Source connectivity analysis from MEG and its application to epilepsy source localization[END_REF][START_REF] De Haan | Disrupted modular brain dynamics reflect cognitive dysfunction in Alzheimer's disease[END_REF][START_REF] Diessen | Brain Network Organization in Focal Epilepsy: A Systematic Review and Meta-Analysis[END_REF][START_REF] Engel | Temporal Binding, Binocular Rivalry, and Consciousness[END_REF][START_REF] Kabbara | Reduced integration and improved segregation of functional brain networks in Alzheimer's disease[END_REF][START_REF] Stam | Small-world networks and functional connectivity in Alzheimer's disease[END_REF]. In this thesis, we focused on this type of connectivity.

Effective Connectivity (EC)

Effective Connectivity (EC) represents the directional effect of one brain region over another, thus it reflects the causal interaction between activated brain regions (Friston, 1994). This connectivity is a dynamic measure as it changes rapidly over time. It can be measured using Granger-causality [START_REF] Granger | Investigating Causal Relations by Econometric Models and Cross-spectral Methods[END_REF] and Dynamic Causal Modeling [START_REF] Friston | Dynamic causal modelling[END_REF].

Neuroimaging techniques

As we have previously stated, one of the main objectives of this thesis is to study the dynamic behavior of the brain. Therefore, we particularly focus here on the principal non-invasive neuroimaging techniques used to map functional brain activity at a large scale. These techniques can be divided into two main classes: (1) indirect measure of neural activity such as fMRI and (2) direct measure of neural activity such as EEG and MEG.

Functional magnetic resonance imaging (fMRI)

fMRI is considered an indirect measure of neural activity since it measures the Blood Oxygen Level-Dependent (BOLD) signals, under the principle that neural activity is supported by increased local blood flow. fMRI is a frequently used neuroimaging technique for several reasons including (i) the widespread availability of MRI scanner, (ii) the relatively low cost per scan, (iii) the safety for scanned individuals, and more importantly (iv) the excellent spatial resolution (order of millimeters) that allows a deep discovery of brain regions activation [START_REF] Christopher Decharms | Applications of real-time fMRI[END_REF]. However, it suffers from relatively low temporal resolution (order of seconds), which prevents tracking dynamic activity at a sub-second time scale [START_REF] Logothetis | Neurophysiological investigation of the basis of the fMRI signal[END_REF].

Electro-encephalography (EEG) / Magneto-encephalography (MEG)

EEG/MEG directly measures the electrical/magnetic activity associated with the activation of neurons clusters [START_REF] Berger | Über das Elektrenkephalogramm des Menschen[END_REF][START_REF] Cohen | Magnetoencephalography: Detection of the Brain's Electrical Activity with a Superconducting Magnetometer[END_REF]. The main interests of using EEG/MEG techniques are related to (i) the easiness of use, (ii) the direct recording of brain time series, (iii) the non-invasiveness, (iv) the cost-effectiveness, and more importantly (v) the excellent temporal resolution (order of milliseconds) that supports real-time detection of fast neural modulations which lies at the heart of dynamic communication in the brain [START_REF] He | Modeling & Imaging of Bioelectrical Activity: Principles and Applications[END_REF][START_REF] Nunez | Electric Fields of the Brain: The neurophysics of EEG[END_REF]. Nevertheless, these techniques have a low spatial resolution (order of centimeters) [START_REF] Nunez | A theoretical and experimental study of high resolution EEG based on surface Laplacians and cortical imaging[END_REF][START_REF] Walsh | Transcranial magnetic stimulation and cognitive neuroscience[END_REF], and activity in deep brain structures is difficult to be directly measured from scalp-level recordings. In our thesis, we adopted both MEG and HD-EEG with 256 electrodes as shown in Figure 4.

Trade-Offs between Temporal and Spatial Resolution

fMRI has exquisite spatial resolution but poor temporal resolution. Conversely, EEG/MEG have a high temporal resolution but lower spatial precision than fMRI (Figure 3). Several studies showed correlation between BOLD signals and EEG/MEG data [START_REF] Hunyadi | A dynamic system of brain networks revealed by fast transient EEG fluctuations and their fMRI correlates[END_REF][START_REF] Prestel | fMRI BOLD Correlates of EEG Independent Components: Spatial Correspondence With the Default Mode Network[END_REF][START_REF] Rizkallah | Removing zero-lag functional connections can alter EEG-source space networks at rest[END_REF].

Although EEG/MEG are sensitive to brain fluctuations at millisecond time-scale, it is hard to pinpoint the cortical sources at the same space resolution as the fMRI. However, giving our interest in analyzing large-scale networks (a set of brain Regions of Interests -ROIs-interacting together during rest and/or task), EEG/MEG can be good enough.

EEG/MEG source connectivity

Previously, most EEG/MEG studies were executed at the sensor-level [START_REF] Vico Fallani | Graph analysis of functional brain networks: practical issues in translational neuroscience[END_REF][START_REF] Garcés | Multimodal description of whole brain connectivity: A comparison of resting state MEG, fMRI, and DWI[END_REF][START_REF] Huang | Sensor Level Functional Connectivity Topography Comparison Between Different References Based EEG and MEG[END_REF][START_REF] Min | Change of Brain Functional Connectivity in Patients With Spinal Cord Injury: Graph Theory Based Approach[END_REF][START_REF] Sato | Association between abnormal brain functional connectivity in children and psychopathology: A study based on graph theory and machine learning[END_REF][START_REF] Schoonheim | Functional connectivity changes in multiple sclerosis patients: a graph analytical study of MEG resting state data[END_REF][START_REF] Wang | Graph theoretical analysis reveals disrupted topological properties of whole brain functional networks in temporal lobe epilepsy[END_REF]. However, in such analysis, signals are contaminated by many factors, which poses serious concerns that can lead to ambiguous interpretations. Specifically, the main factors are (1) the volume conduction problem caused by the electrical conduction properties of the head and (2) the field spread issue, where several sensors can pick up the activity arising from one brain source [START_REF] Nolte | Identifying true brain interaction from EEG data using the imaginary part of coherency[END_REF]Schoffelen and Gross, 2009;[START_REF] Stam | Phase lag index: assessment of functional connectivity from multi channel EEG and MEG with diminished bias from common sources[END_REF][START_REF] Van Diessen | Opportunities and methodological challenges in EEG and MEG resting state functional brain network research[END_REF]. To address these limitations, recent works have developed a particular interest in functional connectivity analysis at the cortical level. In this context, the transition from the sensor to cortical level improves the spatiotemporal resolution of the identified functional networks. This approach is known as EEG/MEG source connectivity (Hassan and Wendling, 2018). It includes two main technical steps: (1)

solve the inverse problem to reconstruct the temporal activity of cortical sources, then (2) compute the functional connectivity between reconstructed sources. These steps are illustrated in Figure 7.

Inverse Problem

According to the dipole theory, EEG/MEG signals 𝑋(𝑡) recorded from 𝑀 sensors can be expressed as a linear combination of 𝑃 dynamic dipolar sources 𝑆(𝑡):

𝑋(𝑡) = 𝐺. 𝑆(𝑡) + 𝑁(𝑡)

Equation 1

Where 𝑁(𝑡) reflects the additive noise, and 𝐺(𝑀 × 𝑃) represents the lead field matrix that reveals the contribution of each cortical source to the sensors. It consists of estimating the properties of template source space (position, orientation, and magnitude). More precisely, the lead field matrix can be calculated from electrodes positions and from a multiple layer head model that describes the electric/magnetic flow from sources to sensors through head tissues (brain, skull, and skin). In this context, Boundary Element Method (BEM) is a widely used numerical approach to build realistic head models while considering detailed characteristics of the head anatomy [START_REF] Fuchs | An improved boundary element method for realistic volume-conductor modeling[END_REF][START_REF] Oostenveld | Validating the boundary element method for forward and inverse EEG computations in the presence of a hole in the skull[END_REF]. As the inverse problem is ill-posed (P>>M), mathematical and physical constraints should be added to limit the approximate solution. In most cases, sources are considered homogeneously distributed over the cortex and constrained normally to the cortical surface. Hence, the inverse problem is reduced to source magnitude estimation as follows:

𝑆(𝑡) = 𝑊. 𝑋(𝑡) Equation 2
Here, W represents the inverse matrix, often called spatial filters or weights. Several algorithms have been proposed to estimate the inverse solution W [START_REF] Baillet | Electromagnetic brain mapping[END_REF], such as beamforming family methods (Van Veen et al., 1997) and least-squares minimum-norm type estimates, including two commonly used methods: (1) weighted minimum norm estimate (wMNE) (Hämäläinen and Ilmoniemi, 1994) and ( 2) standardized/exact low-resolution brain electromagnetic tomography (s/eLORETA) [START_REF] Pascual-Marqui | Low resolution brain electromagnetic tomography (LORETA) functional imaging in acute, neuroleptic-naive, first-episode, productive schizophrenia[END_REF].

After reconstructing sources within a high-resolution mesh surface (~15000 vertices), a source space with precise regions of interest (ROIs) should be defined to compute the corresponding regional time series. This computation can be done based on several proposed approaches. For instance, many studies adopted the averaging process of cortical time series across ROIs [START_REF] Hasson | A Hierarchy of Temporal Receptive Windows in Human Cortex[END_REF][START_REF] Hay | Multiregional integration in the brain during resting-state fMRI activity[END_REF][START_REF] Yang | Disentangling time series between brain tissues improves fMRI data quality using a time-dependent deep neural network[END_REF], while others chose to select a single dipole within the ROI that represents the whole brain region activity [START_REF] Coito | Directed Functional Brain Connectivity Based on EEG Source Imaging: Methodology and Application to Temporal Lobe Epilepsy[END_REF][START_REF] Ghumare | A Time-Varying Connectivity Analysis from Distributed EEG Sources: A Simulation Study[END_REF][START_REF] Sohrabpour | Noninvasive Electromagnetic Source Imaging and Granger Causality Analysis: An Electrophysiological Connectome (eConnectome) Approach[END_REF]. The representative dipole can be selected based on ground truth data [START_REF] Babiloni | Multimodal integration of EEG and MEG data: a simulation study with variable signal-tonoise ratio and number of sensors[END_REF], or data-driven methods [START_REF] Farahibozorg | Adaptive cortical parcellations for source reconstructed EEG/MEG connectomes[END_REF][START_REF] Rueda-Delgado | Coordinative task difficulty and behavioural errors are associated with increased long-range beta band synchronization[END_REF][START_REF] Sohrabpour | Noninvasive Electromagnetic Source Imaging and Granger Causality Analysis: An Electrophysiological Connectome (eConnectome) Approach[END_REF]. Usually, ROIs can be defined based on anatomical brain parcellations such as Desikan-Killiany atlas (68 ROIs) (Desikan et al., 2006), Destrieux atlas (148 ROIs) (Destrieux et al., 2010), Automated Anatomical

Labeling AAL atlas (78 ROIs) (Gong et al., 2009;[START_REF] Tzourio-Mazoyer | Automated Anatomical Labeling of Activations in SPM Using a Macroscopic Anatomical Parcellation of the MNI MRI Single-Subject Brain[END_REF] or based on functional parcellations such as Brainnetome functional atlas (246 ROIs) [START_REF] Yu | Functional segregation of the human cingulate cortex is confirmed by functional connectivity based neuroanatomical parcellation[END_REF] (see Figure 5). Besides, another data-driven strategy can be adopted to decompose the cortex into specific regions using General Linear Model (GLM) or spatial independent component analysis (sICA) [START_REF] Tsvetanov | Extrinsic and Intrinsic Brain Network Connectivity Maintains Cognition across the Lifespan Despite Accelerated Decay of Regional Brain Activation[END_REF][START_REF] Xu | Spatial ICA reveals functional activity hidden from traditional fMRI GLM-based analyses[END_REF].

Functional Connectivity (FC) estimation

Once cortical time series are reconstructed, FC can be estimated as the statistical couplings between regional signals. In this context, a variety of methods have been proposed to characterize and quantify the large-scale functional organization of the human brain based on linear/non-linear, phase/amplitude, temporal/spectral measures [START_REF] Behrens | Human connectomics[END_REF][START_REF] Craddock | Imaging human connectomes at the macroscale[END_REF][START_REF] Haufe | A critical assessment of connectivity measures for EEG data: A simulation study[END_REF]Lachaux et al., 1999;[START_REF] Stam | Phase lag index: assessment of functional connectivity from multi channel EEG and MEG with diminished bias from common sources[END_REF]Vinck et al., 2011). The reader can refer to [START_REF] Pereda | Nonlinear multivariate analysis of neurophysiological signals[END_REF] for a review and [START_REF] Wendling | From EEG signals to brain connectivity: a model-based evaluation of interdependence measures[END_REF] for a comparative study. As a result, an adjacency matrix (𝑁 𝑟 × 𝑁 𝑟 ) is produced, where 𝑁 𝑟 refers to the number of ROIs. In general, a threshold is applied on this matrix in order to remove spurious connections and preserve relevant ones. The selection of threshold type (proportional/absolute/statistical) and corresponding value is critical and is the subject of many research works (Garrison et al., 2015;[START_REF] Genovese | Thresholding of statistical maps in functional neuroimaging using the false discovery rate[END_REF][START_REF] Van Den Heuvel | Proportional thresholding in resting-state fMRI functional connectivity networks and consequences for patient-control connectome studies: Issues and recommendations[END_REF].

Consequently, the obtained connectivity matrix can be assimilated to a graph, where brain regions are represented by nodes and FC values are represented by edges. Based on the applied threshold, the obtained brain network can be categorized into binary/weighted graphs. In this thesis, we will focus on weighted brain networks (see Figure 6). static Functional Connectivity (sFC) :

Primarily, a large amount of studies has investigated the static connectivity and characterized their strength and fragility (Brookes et al., 2011b;[START_REF] Colclough | How reliable are MEG resting-state connectivity metrics?[END_REF][START_REF] Garcés | Multimodal description of whole brain connectivity: A comparison of resting state MEG, fMRI, and DWI[END_REF]Hassan et al., 2017a;Liuzzi et al., 2019;[START_REF] Pereda | Nonlinear multivariate analysis of neurophysiological signals[END_REF][START_REF] Wang | Graph theoretical analysis reveals disrupted topological properties of whole brain functional networks in temporal lobe epilepsy[END_REF]. Refer to [START_REF] Bastos | A Tutorial Review of Functional Connectivity Analysis Methods and Their Interpretational Pitfalls[END_REF][START_REF] Friston | Functional and effective connectivity: a review[END_REF] for review. Within this framework, the brain is assumed to remain in a state of static equilibrium, where functional connectivity was considered constant over the entire recording period.

Thus, it was calculated and averaged over all time points to produce a single representative measure of FC, yielding to the so-called 'static Functional Connectivity (sFC)'. Despite its ability to explore relevant topological properties of neural activity [START_REF] Hutchison | Dynamic functional connectivity: promise, issues, and interpretations[END_REF], sFC loses information at a temporal scale, since functional brain networks exhibit spontaneous changes across multiple time scales [START_REF] Chang | Time-frequency dynamics of resting-state brain connectivity measured with fMRI[END_REF][START_REF] Hutchison | Dynamic functional connectivity: promise, issues, and interpretations[END_REF]O'Neill et al., 2018). This has opened the gates to the recently emerging technique, known as 'dynamic Functional Connectivity (dFC)'.

dynamic Functional Connectivity (dFC) :

Currently, a growing body of research has been directed to considering FC modulations across time as physiologically informative during mental and attention tasks [START_REF] Esposito | Independent component model of the default-mode brain function: Assessing the impact of active thinking[END_REF]Fong et al., 2019;[START_REF] Fornito | Competitive and cooperative dynamics of large-scale brain functional networks supporting recollection[END_REF][START_REF] Sun | Functional connectivity of cortical networks involved in bimanual motor sequence learning[END_REF], learning [START_REF] Albert | The Resting Human Brain and Motor Learning[END_REF]Bassett et al., 2011b;[START_REF] Lewis | Learning sculpts the spontaneous activity of the resting human brain[END_REF], sleep [START_REF] Horovitz | Low frequency BOLD fluctuations during resting wakefulness and light sleep: a simultaneous EEG-fMRI study[END_REF], and anesthesia [START_REF] Peltier | Functional connectivity changes with concentration of sevoflurane anesthesia[END_REF]. In contrast to sFC, dFC accounts for the dynamic reorganization of functional brain networks. Hence, it can reflect the genuine timevarying characteristics of the complex brain networks at a large scale. The reader can refer to [START_REF] Calhoun | The Chronnectome: Time-Varying Connectivity Networks as the Next Frontier in fMRI Data Discovery[END_REF][START_REF] Hutchison | Dynamic functional connectivity: promise, issues, and interpretations[END_REF]O'Neill et al., 2018;[START_REF] Preti | The dynamic functional connectome: State-of-the-art and perspectives[END_REF] for technical reviews on dFC.

In this framework, the sliding window approach is a widely used technique to compute dFC in both fMRI (Allen et al., 2014;[START_REF] Chang | EEG correlates of time-varying BOLD functional connectivity[END_REF][START_REF] Elton | Task-related modulation of functional connectivity variability and its behavioral correlations[END_REF][START_REF] Hutchison | Dynamic functional connectivity: promise, issues, and interpretations[END_REF][START_REF] Kucyi | Dynamic functional connectivity of the default mode network tracks daydreaming[END_REF] and MEG/EEG literature (Baker et al., 2014;[START_REF] Brookes | Measuring temporal, spectral and spatial changes in electrophysiological brain network connectivity[END_REF]Brookes et al., , 2011a;;[START_REF] Carbo | Dynamic hub load predicts cognitive decline after resective neurosurgery[END_REF]de Pasquale et al., 2016de Pasquale et al., , 2010;;[START_REF] Doron | Dynamic network structure of interhemispheric coordination[END_REF][START_REF] Lee | Diversity of functional connectivity patterns is reduced in propofol-induced unconsciousness[END_REF]O'Neill et al., 2017O'Neill et al., , 2015;;[START_REF] Yang | Differences in early dynamic connectivity between visual expansion and contraction stimulations revealed by an fMRI-directed MEG approach[END_REF]. This approach has offered promising findings to the neuroscience community, relating brain dynamics to cognition [START_REF] Elton | Task-related modulation of functional connectivity variability and its behavioral correlations[END_REF][START_REF] Kucyi | Dynamic functional connectivity of the default mode network tracks daydreaming[END_REF][START_REF] Madhyastha | Age-related differences in the dynamic architecture of intrinsic networks[END_REF] as well as brain disorders [START_REF] Jie | Integration of temporal and spatial properties of dynamic connectivity networks for automatic diagnosis of brain disease[END_REF]Leonardi et al., 2013;[START_REF] Sakoğlu | A method for evaluating dynamic functional network connectivity and task-modulation: application to schizophrenia[END_REF]. Technically, it consists of simple steps. First, a temporal window, parameterized by its predefined length W, is selected, and a functional connectivity matrix (𝑁 𝑟 × 𝑁 𝑟 ) (considering 𝑁 𝑟 brain regions) is computed over temporal window length [1, 𝑊]. Then, the window is shifted in time by a predefined step T, where the same calculation is repeated over the time interval [1 + 𝑇, 𝑊 + 𝑇]. The same procedure is repeated until the sliding window spans the end part of time courses. Finally, the process produces a three-dimensional tensor of windowed FC matrices (𝑁 𝑟 × 𝑁 𝑟 × 𝑁 𝑤 ), where 𝑁 𝑤 refers to the total number of windows (see Figure 7). The resultant tensor represents the so-called dynamic functional connectivity (dFC).

Dynamic Brain Network States

dFC results in hundreds or thousands of functional brain networks fluctuating over time windows.

Therefore, it is important to point out that modeling the full dynamic brain connectivity is difficult to be interpreted and thus has to face the curse of dimensionality to aid analysis on results. Interestingly, several studies showed that dFC organizes itself into reproducible activity patterns. Hence, an increasing amount of research suggests further time-resolved analysis on the high-dimensional dFC results. It involves tracking recurring spatial FC configuration over time and summarizing it into a smaller set of dominant brain patterns with their temporal evolution [START_REF] Calhoun | Time-Varying Brain Connectivity in fMRI Data: Whole-brain data-driven approaches for capturing and characterizing dynamic states[END_REF][START_REF] Iraji | The connectivity domain: Analyzing resting state fMRI data using feature-based datadriven and model-based methods[END_REF]O'Neill et al., 2018). These dynamic patterns are denoted 'dynamic Brain Network States (dBNS)'. To this end, windowed dFC matrices are first concatenated along temporal dimension for all experimental trials. Then, the concatenated matrices represent the input into summary statistics using computational methods of dimensionality reduction. Two broad classes of computational methods may be identified, namely knowledge-based and data-driven approaches. The latter approach, also termed as unsupervised or exploratory methods, can be generally subdivided further into decomposition methods and clustering techniques [START_REF] Lachaux | Brain Connectivity Analysis: A Short Survey[END_REF][START_REF] Li | Review of methods for functional brain connectivity detection using fMRI[END_REF]. Hereafter, we will provide an overview of both model-based and data-driven approaches used to derive dBNS.

Knowledge-based approach

Knowledge-based methods afford strong prior neuroscience knowledge about spatiotemporal patterns activation, as well as a model for the data generation procedure. Some studies primarily select some ROIs as 'seeds', and detect whether other regions are functionally connected to these seeds, using predefined metrics, such as cross-correlation analysis, coherence analysis, and statistical parametric mapping, knowledge-based techniques, namely CCA [START_REF] Fox | The human brain is intrinsically organized into dynamic, anticorrelated functional networks[END_REF][START_REF] Greicius | Default-mode network activity distinguishes Alzheimer's disease from healthy aging: evidence from functional MRI[END_REF][START_REF] Johansen-Berg | Changes in connectivity profiles define functionally distinct regions in human medial frontal cortex[END_REF]. Others define a prior hypothesis to model dFC data as the GLM model (Friston, 1994;[START_REF] Iraji | The connectivity domain: Analyzing resting state fMRI data using feature-based datadriven and model-based methods[END_REF]. Nevertheless, this approach poses some concerns: (1) results depend on different seed choices leading to different connectivity maps [START_REF] Ma | Detecting functional connectivity in the resting brain: a comparison between ICA and CCA[END_REF], ( 2) unexpected yet instructive connectivities related to other brain parts, not considered in the model, may be neglected, and (3) results dependency on the selected fitting model. Therefore, in this thesis, we focus on the full exploration of dFC through data-driven-based studies.

Data-driven approach

To overcome model-based methods limitations, exploratory techniques that are independent of any prior knowledge or hypothesis have been exploited. We refer the reader to [START_REF] Calhoun | Time-Varying Brain Connectivity in fMRI Data: Whole-brain data-driven approaches for capturing and characterizing dynamic states[END_REF] for a review and to (Miller et al., 2016;[START_REF] Xie | Decoding the encoding of functional brain networks: An fMRI classification comparison of non-negative matrix factorization (NMF), independent component analysis (ICA), and sparse coding algorithms[END_REF] for comparative studies on data-driven methods applied in the context of fMRI data. This approach mainly relies on the assumption that the brain is organized in a finite set of functional networks. Hence, it can be formulated as a source separation problem that aims to recover finite hidden states from a set of observations with minimal prior information as follows:

𝑀 = 𝐴 × 𝐵 Equation 3
Where M refers to the concatenated dFC matrices, A represents the mixing matrix illustrating the dominant spatial maps, and B represents the corresponding temporal signatures (see Figure 8). In general, there exists two types of data-driven methods for functional connectivity analysis. The first type involves decomposition techniques that consider any observation as a linear combination of underlying features, such as Principal Component Analysis (PCA), Independent Component Analysis (ICA), and Non-negative Matrix Factorization (NMF). The second type is based on clustering techniques, such as fuzzy/hierarchical clustering analysis. They are all supposed to capture relevant dBNS, however, they mainly differ in the constraints imposed on derived states. Below, we will present a succinct description of these methods with their contribution to dBNS literature.

Principal Component Analysis (PCA): 'Variance maximization'

PCA is a basic technique used for dimensionality reduction issues. It transforms functional images into uncorrelated, orthogonal components, denoted 'eigenvectors'. Technically, this decomposition is based on data variance maximization. It has been used to separate low and high variance components as a way to extract the key features with a minimal redundancy to understand normal neural activity [START_REF] Leonardi | Disentangling dynamic networks: Separated and joint expressions of functional connectivity patterns in time[END_REF](Leonardi et al., , 2013;;[START_REF] Mutlu | A signal-processing-based approach to time-varying graph analysis for dynamic brain network identification[END_REF], as well as altered activity due to brain pathologies such as multiple sclerosis (Leonardi et al., 2013) and Schizophrenia (Miller et al., 2016).

Independent Component Analysis (ICA): 'Temporal Statistical Independency'

While PCA decorrelates second-order dependencies, ICA usually uses higher-order statistics to produce statistically independent features, under the hypothesis that these features are as 'non-Gaussian' as possible. Several algorithms have been proposed to perform ICA decomposition, yielding to different statistical models with specific properties. Some of these algorithms are based on information theory, while others try to optimize contrast functions based on second to fourth-order cumulants. The approach of temporal independent component (tICA) separates brain networks with distinct and independent temporal signatures. It has been applied on dFC to identify dBNS based on common changes in connectivity across multiple connections between brain regions (Miller et al., 2016;O'Neill et al., 2017;[START_REF] Stone | Spatiotemporal Independent Component Analysis of Event-Related fMRI Data Using Skewed Probability Density Functions[END_REF][START_REF] Theis | A robust model for spatiotemporal dependencies[END_REF]Yaesoubi et al., 2015).

Non-negative Matric Factorization (NMF): 'Positivity'

NMF is an unsupervised technique commonly used in machine learning (Lee and Seung, 1999). It has been employed by [START_REF] Chai | Evolution of brain network dynamics in neurodevelopment[END_REF] as a matrix factorization technique that constraints 'positivity' or 'non-negativity' on the extracted dBNS. Specifically, Alternating Least Squares (ALS) algorithm has proven its robustness among existing algorithms in the context of functional brain states identification (Ding et al., 2013).

Kmeans Clustering: 'Sparsity'

Kmeans is one of the simplest and most popular algorithms used to solve clustering problems (Lloyd, 1982). In our framework, it works iteratively to identify clusters that contain spatial maps with similar topologies across time. Then, it assigns functional networks at each temporal window to one of the identified clusters based on features similarity, which generates a probabilistic time course across trials. This technique has often been applied in the context of fMRI (Allen et al., 2014;[START_REF] Gonzalez-Castillo | Tracking ongoing cognition in individuals using brief, whole-brain functional connectivity patterns[END_REF][START_REF] Hutchison | Tracking the Brain's Functional Coupling Dynamics over Development[END_REF][START_REF] Shakil | Evaluation of sliding window correlation performance for characterizing dynamic functional connectivity and brain states[END_REF]), MEG (O'Neill et al., 2015), and EEG analysis (Hassan et al., 2015;Mheich et al., 2015;[START_REF] Zoltowski | A graph theoretic approach to dynamic functional connectivity tracking and network state identification[END_REF]. In addition, several studies revealed the clinical usefulness of clustering approaches to detect alterations in dBNS relative to brain pathologies such as schizophrenia [START_REF] Damaraju | Dynamic functional connectivity analysis reveals transient states of dysconnectivity in schizophrenia[END_REF]Du et al., 2016;Miller et al., 2016), Parkinson [START_REF] Fiorenzato | Dynamic functional connectivity changes associated with dementia in Parkinson's disease[END_REF]Kim et al., 2017) and Autism Spectrum Disorders [START_REF] Wee | Diagnosis of Autism Spectrum Disorders Using Temporally Distinct Resting-State Functional Connectivity Networks[END_REF].

Thesis Objective

One purpose of this thesis is to further explore to what extent can we track the dynamic brain network states (dBNS) when performing tasks at a rapid time scale, in healthy subjects as well as patients with neurological disorders. Specifically, the main objectives are: To accomplish these goals, we used non-invasive neuroimaging techniques with high-temporal resolution (EEG/MEG). Mainly, we applied the dynamic version of a recently developed approach, known as 'EEG/MEG source connectivity', along with decomposition techniques to estimate dynamic brain network states.


Briefly, in study I, we evaluated the performance of nine different source separation techniques using three independent empirical MEG datasets (N=95 subjects) during motor and working memory tasks.

The goal was to explore the ability of different decomposition assumptions to decipher the dynamics of functional connectivity patterns in the context of behavioral tasks. These methods were tested at both group and subject levels. In this work, we reported some variability between methods that show dependency on the task complexity and temporal scale. Promising results were observed using the Independent Component Analysis technique.

Then, a further quantitative evaluation of the essential decomposition methods was investigated in study II, using a physiologically inspired full brain model that provides ground-truth data. In this study, a virtual picture naming scenario was modeled for 20 subjects at the cortical level. We aimed to provide a framework that can be used to optimize the several steps pipeline applied to accurately re-estimate the reference modeled data in terms of spatial and temporal accuracy. Besides decomposition techniques evaluation, different combinations of inverse models and connectivity measures were tested. Our findings suggest the existence of significant variability in the spatio-temporal precision among the tested algorithms.

Finally, the previous approach was applied in study III to detect dysfunctionality induced by Parkinson's disease (PD). Using High Density EEG (HD-EEG) recorded from 10 healthy subjects and 21 patients during the Simon task, we hypothesized the existence of alterations in the BNS of PD patients. This hypothesis was then validated by the means of tICA technique combined with source-level microstate analysis, which revealed some significant differences among the extracted dBNS between HC and PD groups.

Chapter

Introduction

Evolving evidence show that motor, sensory, emotional and cognitive functions emerge from dynamic interactions between cortical and subcortical brain structures. Specific rhythms of neural networks allow synchronization and long-range communication between distant and distributed brain areas. This phenomena was shown crucial during visual ( Bola and Sabel, 2015 ;Hassan et al., 2015 ;Mheich et al., 2018 ), auditory [START_REF] Fontolan | The contribution of frequency-specific activity to hierarchical information processing in the human auditory cortex[END_REF], sensorimotor [START_REF] Pomper | Intersensory selective attention and temporal orienting operate in parallel and are instantiated in spatially distinct sensory and motor cortices: human Intersensory and Temporal Attention[END_REF]Wilkins and Yao, 2020 ) and cognitive [START_REF] Negrón-Oyarzo | Coordinated prefrontal-hippocampal activity and navigation strategy-related prefrontal firing during spatial memory formation[END_REF][START_REF] Rouhinen | Spectral and anatomical patterns of large-scale synchronization predict human attentional capacity[END_REF] tasks. This brain communication is very transient and there is a dynamic reorganization of functional brain networks during behavioral tasks, even at sub-second time scale ( Vidaurre et al., 2018b ). Therefore, the analysis of whole-brain dynamic functional connectivity (dFC) has become a burgeoning field of research in cognitive neuroscience [START_REF] Bassett | Network neuroscience[END_REF]Bullmore and Sporns, 2009 ;Iraji et al., 2020 ;[START_REF] Kabbara | The Dynamic Modular Fingerprints of the Human Brain at Rest bioRxiv[END_REF]. In this regard, Magneto/Electroencephalography (MEG/EEG) provides a unique direct and noninvasive access to the electrophysiological activity of the whole brain, at the millisecond scale. Benefiting from the excellent time resolution of the MEG/EEG (~millisecond), current methods allow of estimating sub-second time-varying functional brain networks in the cortical space through sensor-level signals ( Hassan et al., 2014 ;Hassan and Wendling, 2018 ). The key challenge here is how to characterize and quantify these rapidly changing networks.

In this context, several frameworks have been used to explicitly model/capture dynamics over time such as Hidden Markov Model (HMM) ( Baker et al., 2014 ;Vidaurre et al., 2018aVidaurre et al., , 2018b[START_REF] Vidaurre | Spectrally resolved fast transient brain states in electrophysiological data[END_REF], Autoregressive model (AR) ( Casorso et al., 2019 ) and General Linear Model (GLM) ( Friston, 1994 ). For example, HMM describes the brain activity as a sequence of district states; each represents a unique pattern obtained from an observation model, and a state time course indicating the points in time at which that state is active. Other approaches analyze the time varying signal using data-driven techniques, where 'brain network states' are derived directly from the data without a priori hypothesis on the fitting model. These methods have showed promising results, despite the fact that the selection of the used algorithm is largely empirical. These methods are based on two main steps: (1) sliding window approach, that forms a series of temporal networks, (2) a dimensionality reduction or clustering approach including Kmeans ( Allen et al., 2014 ;[START_REF] Ciric | Contextual connectivity: A framework for understanding the intrinsic dynamic architecture of large-scale functional brain networks[END_REF]Du et al., 2016 ;Fong et al., 2019 ;[START_REF] Liu | Time-varying functional network information extracted from brief instances of spontaneous brain activity[END_REF]Mheich et al., 2015 ;O'Neill et al., 2015 ), component analysis such as temporal Independent Component Analysis tICA ( O'Neill et al., 2017 ), Principal Component Analysis (PCA) ( Leonardi et al., 2013 ) and Nonnegative Matrix Factorization (NMF) [START_REF] Chai | Evolution of brain network dynamics in neurodevelopment[END_REF]. Although the conceptual difference between these methods (and within each family of methods such as different ICA algorithms) is theoretically obvious (as they are based on different assumptions), the studies that investigate the differences between them remained very few. The existing comparative studies are mainly limited to confirming results of differences between two conditions ( Leonardi et al., 2013 ) or to prove that obtained results are unaffected by the method's choice ( Miller et al., 2016 ). However, a throughout quantitative and qualitative comparative study using both simulation-based and data-driven approaches is still missing and there is no clear consensus about the 'best' (if any) source separation or clustering method to be used to adequately tracking dFC, which is the main objective of our study.

Here, we evaluate the performance of nine dimensionality reduction methods used to track functional connectivity states at both group and individual levels. This was done using simulations and three independent MEG datasets (N = 95) recorded during motor and working memory tasks (see Fig. 1 ). The dynamic brain networks were reconstructed using MEG source connectivity method combined with a sliding window technique. The dimensionality reduction algorithms were compared in terms of their temporal and spatial accuracy. These methods include PCA, NMF, Kmeans and six various versions of ICA (Joint Approximation Diagonalization of Eigen-Matrices (JADE), INFOMAX, Second-Order Blind Identification (SOBI), fixed-point algorithm (Fas-tICA), COM2 and Penalized Semi-Algebraic Unitary Deflation (P-SAUD). The motivation behind using several ICA subtypes is that each one has its own definition of statistical independence and several studies showed conceptual differences between them [START_REF] Kachenoura | ICA: a potential tool for bci systems[END_REF][START_REF] Sahonero-Alvarez | A comparison of SOBI, FastICA[END_REF]. We also analyzed the optimal number of subjects needed for each method to reveal significant results. This study aims at providing a framework for researchers interested in studying reconfiguration of functional brain networks during cognitive processes.

Materials and methods

Data

Dataset1: 'self-paced button press task'

This dataset includes 15 healthy righthanded participants (9 male and 6 female, aged 25 ± 4 years (mean ± SD)). They were asked to press a button with the index finger of their non-dominant hand, once every 30 seconds, and should not count the time between presses. More details about this dataset can be found in ( Kabbara et al., 2019 ;O'Neill et al., 2017 ).

Dataset2: 'HCP left hand movement Task'

61 healthy participants (28male and 33 female, aged 22-35) completed the MEG Motor task provided by the Human Connectome Project (HCP) (MEG-1 release) ( Van Essen et al., 2012 ). The corresponding experimental protocol was adapted from Buckner and colleagues ( Buckner et al., 2011 ;[START_REF] Yeo | The organization of the human cerebral cortex estimated by intrinsic functional connectivity[END_REF]. It was performed in two sessions of 14min each, with a small break between them. Each session consisted of 42 total blocks randomly distributed; 32 of them were partitioned into 16 hand movements blocks (8 right and 8 left), and 16-foot movements blocks (8 right and 8 left), and the remaining 10 blocks were interleaved resting/fixation blocks. Each motor effector block was preceded by a 3sec visual cue that prompts participants to either tap their left or right index and thumb fingers or squeeze their left or right toes. The block lasted for 12sec and consisted of 10 sequential movements, each initiated with 150ms pacing stimuli followed by 1050ms black screen for task execution. Here, for simplicity, we were interested in the trials related to the left hand moves only. MEG data was recorded at Saint Louis University at 508.6275Hz sampling frequency and co-registered with the available subject specific MRI. EMG activity was also recorded from each limb.

Dataset3: 'sternberg working memory task'

19 healthy participants (10 male and 9 female, aged 25 ± 3 years (mean ± SD)) performed Sternberg task, in which two example visual stimuli, mainly abstract geometric shapes, were successively presented on a screen; each for 0.6sec and separated by 1sec. Then, a maintenance period of 7sec was left before the presentation of a third probe stimulus. Consequently, subjects were asked to press a button with their right index finger only if the probe stimulus matched either of the two example stimuli and an immediate feedback will be given to show their response correctness. 30 trials were presented separated by 30sec of rest. In both datasets 1 and 3, MEG data were recorded using a 275-channel CTF MEG system at 600Hz sampling frequency and co-registered with subject-specific MRI. Both datasets were approved by the University of Nottingham Medical School Research Ethics Committee ( O'Neill et al., 2017 ;Vidaurre et al., 2018a ).

Methodology

Preprocessing

Both datasets 1 and 3 were received already preprocessed as described in ( O'Neill et al., 2017 ). Briefly, bad segments produced by muscles, eye or head movement were already visually inspected and removed. For dataset 2, we used the preprocessing pipeline offered by the HCP consortium, which includes removing bad channels, segments and bad independent components from task data. Segments were retrieved from the dataset 1 in the interval [-15; + 15sec] relative to the button press onset, and from the dataset 3 in the interval of [-16; + 28sec] relative to stimulus presentation. In HCP analysis, we chose data epochs time-locked to EMG onset as we were concerned in exploring brain networks involved during movement execution. Thus, trials were segmented in [-1.2; + 1.2sec] relative to EMG onset. Then, as functional connectivity was proved to be frequency-dependent ( Baker et al., 2014 ;Hipp et al., 2012 ), each dataset was preprocessed in its appropriate frequency band actively involved in the corresponding cognitive task. While beta band was used for self-paced and HCP left hand motor task, working memory data was filtered in a broader band as it is has been shown to involve multiple frequency bands, according to previous studies ( Brookes et al., 2012 ;O'Neill et al., 2017 ). After these preprocessing steps, an average of 34, 150 and 29 per subject were kept from dataset 1, 2 and 3, respectively.

Source reconstruction and functional connectivity

In order to localize brain sources and reconstruct their activities, we used the Linearly Constrained Minimum Variance Beamforming (LCMV) ( ROBINSON, 1999 ) approach on parcellated cortex using AAL atlas (N = 78 regions of interests -ROIs- ( Gong et al., 2009 )) ( Hillebrand et al., 2016 ). This was done by registering each subject's anatomical MRI to an MNI template ( Smith et al., 2004 ) followed by an inverse registration to the anatomical subject space. Data covariance was computed within the specific frequency band used and a time window spanning the whole experiment ( Brookes et al., 2008 ) with a regularization parameter (5%) using Tikhonov method. The forward model was based upon a dipole approximation [START_REF] Sarvas | Basic mathematical and electromagnetic concepts of the biomagnetic inverse problem[END_REF] and a multiple local sphere head model fitted to the subject-specific MRI scalp surface. Dipole orientation was determined using a non-linear search for optimum 'signal to noise ratio' (SNR) [START_REF] Sekihara | Adaptive spatial filters for electromagnetic brain imaging[END_REF]. Following this, we estimated the functional connectivity by computing the amplitude envelope correlations (using Hilbert transformation) between all ROIs ( Brookes et al., 2012 ;Hipp et al., 2012 ). In order to avoid spurious estimates of functional connectivity, we performed leakage correction on the reconstructed sources signals. We used the multivariate approach based on symmetric orthogonalisation proposed by ( Brookes et al., 2012 ;[START_REF] Colclough | A symmetric multivariate leakage correction for MEG connectomes[END_REF] for datasets 1 and 3, while pair-wise orthogonalization ( Brookes et al., 2016 ;Tewarie et al., 2019b ) was applied to dataset 2 due to the short time period of the task.

Fig. 1. Illustration of the investigation structure for each of the three task-related paradigms. A. The fundamental processing pipeline applied on each subject data from sensor-level (using non-invasive MEG technique) to cortical-level (using beamforming as the inverse problem solution) to dynamic functional connectivity computation (S-dFC) (using the sliding window approach). B. Concatenation of S-dFC of all subjects along time axis to form a group data referred to as G-dFC, C. Comparative analysis between nine different source separation (SS) methods (six variants of ICA, PCA, NMF and Kmeans) applied on both group-level (X group ) and subject-level (X subj ) data in order to derive k dominant task-related spatiotemporal components (the mixing matrices represent brain spatial maps while the extracted sources represent corresponding temporal weights fluctuations).

Dynamic functional connectivity analysis (dFC)

To estimate the dynamic functional brain networks, we adopted the widely used approach of sliding windows for datasets1 and 3. To this end, a time window of length 6sec with 0.5sec was used for datasets 1 and 3 as applied by ( O'Neill et al., 2017 ). Concerning the dataset 2 (HCP dataset), the fast time scale of the task imposes a very small time window width that may be too noisy to extract meaningful information ( Liuzzi et al., 2019 ). Thus, we avoided to apply the sliding window approach, and used instead the high temporal resolution version of amplitude envelope correlation metric; the 'Instanteneous Amplitude Correlation' (IAC) already validated in a recent work for the same dataset ( Tewarie et al., 2019b ). As a result, we obtained, for each subject trial, a 'dynamic functional connectivity (dFC)' matrix of dimension [NxNxT], where T refers to the number of windows for datasets 1 and 3, and number of total time samples for dataset 2 (T = 49, 1221 and 75 for datasets 1, 2 and 3 respectively). Next, due to symmetry, we unfolded this matrix into a 2-D [Nx(N -1)/2 × T] matrix by removing the redundant connections in each time window. Then, the mean of each row of this matrix is subtracted from the data. Finally, all subjects' trials dFC were concatenated along the temporal dimension.We defined this matrix as a 'Group dynamic functional connectivity matrix (G-dFC)', denoted 'X'. Note that for dataset 2, we averaged connectivity matrices of all trials relative to each subject ( Zhu et al., 2020 ) due to memory limitation in Matlab regarding high dimensional data of the temporally concatenated 'sample-by-sample' dFC of all subject's trials.

Task-related functional brain networks 2.2.4.1. Problem statement.

The resultant G-dFC matrix representing the time-varying features can be expressed as as a linear mixture of elementary brain networks that fluctuate dynamically over time. Such issue is the main concern of Source Separation (SS) approach aiming at recovering 'k' hidden sources from a set of observations with minimal priori knowledge about these sources. In this context, the SS problem can be formulated as follows:

X = A × S (1)
Where:

• 'X' is the computed G-dFC matrix of dimension [qxm]: ○ q = Nx(N -1)/2 with N = 78, representing connectivities between all ROIs. ○ m = TxN tot with T is the number of time windows and N tot is the total number of trials for all subjects. • 'A' is the mixing matrix of dimension [qxk] illustrating the contribution weights of each individual connection to the components sources, thus the spatial maps of brain networks (k < min(q,m)). • 'S' is the sources matrix of dimension [kxm] representing temporal sources signatures of G-dFC, collapsed across all connections.

Among existing SS algorithms, we chose nine popular/well-known methods: six different variants of temporal Independent Component Analysis (tICA), Principal Component Analysis (PCA), Non-negative Matrix Factorization (NMF) and Kmeans as a state-of-the-art clustering method. They all transform the desired matrix factorization into spatial maps and time series. However, they differ primarily in the constraints imposed on decomposed components. Below, we will give a succinct description about these methods.

Independent component analysis: 'temporal statistical independence'.

ICA tends to linearly transform multivariate observations into a set of 'statistically mutually independent' latent variables under the hypothesis that these variables are as 'non-Gaussian' as possible. In our study, we examine temporal ICA (tICA) adopted by several previous studies ( O'Neill et al., 2017 ;Yaesoubi et al., 2015 ) in order to obtain states that fluctuate independently in time. In this context, decomposed signals 'S' consist of the 'k' source time courses and the associated mixing matrix 'A' illustrates the contribution of temporally independent maps.

There are several criteria to measure independence such as minimization of mutual information and maximization of non-Gaussianity. Hence, different algorithms are proposed to perform ICA decomposition, each yielding to different ICA model with specific characteristics. Here, we evaluate tICA using six different popular and prominent methods:

(1) JADE, (2) InfoMax, (3) SOBI, (4) FastICA, (5) CoM2 and (6) PSAUD. These methods are chosen in such a way to cover various statistical independence definitions, statistical order and computational process techniques. Briefly, InfoMax and FastICA are based on information theory, while all other selected methods optimize contrast functions based on cumulants of the data. Among them, SOBI uses only Second Order (SO) cumulants in contrast to others that exploit both SO and Fourth Order (FO) cumulants. In addition, FastICA and PSAUD use a deflation process for decomposition while other ICA variants jointly separate sources. Details about ICA subtypes used can be found in supplementary materials.

Principal component analysis: 'variance maximization'.

PCA is a basic linear technique widely used for data dimensionality reduction. It involves a mathematical procedure that transforms a set of observations of possibly correlated variables into smaller number of orthogonal, hence linearly uncorrelated variables called principal components or 'eigenvectors. This procedure is defined in such a way that the variance or 'eigenvalues' of the data is maximized. Then, a fixed number 'k' of eigenvectors and their respective eigenvalues can be chosen to obtain a consistent representation of the data. Here, we apply the Singular Value Decomposition (SVD) algorithm of PCA ( Golub and Reinsch, 1970 ) on our predefined input matrix 'X'. Defining 'A' and 'S' matrices from SVD outputs is more clarified in Supplementary Materials.

2.2.4.4. Non-negative matrix factorization: 'positivity'. Nonnegative matrix factorization (NMF) is an unsupervised machine-learning technique ( Lee and Seung, 1999 ) that imposes 'non-negativity' constraint on the decomposed factors when solving SS problem. When applied to G-dFC data 'X', NMF leads to parts-based representation that captures additive combination of basis subgraphs 'A' at each time window with temporal coefficients 'S' eliminating negative signal variations. Among several existing NMF approaches, we selected Alternating Least Squares (ALS) algorithm that has previously shown good performance in fMRI context ( Ding et al., 2013 ) with 100 times replications.

2.2.4.5. Kmeans clustering: 'sparsity'. Kmeans is one of the simplest unsupervised learning algorithms that solve the SS problem through clustering approach ( Lloyd, 1982 ). The algorithm works iteratively to assign each point to only one of the 'k' groups based on feature similarity. Mathematical computation of Kmeans clusters is defined in Supplementary Materials. In our framework, the sparse coding adopted by Kmeans restricts a single time point to have a unique activated network state. The computed clusters 'A' represents the structure of common connectivity patterns across subjects. For a given trial, each time window is assigned with the corresponding cluster index. Then, the matrix 'S' is calculated as the frequency of reoccurrence of each cluster at each time window across all trials and subjects. Here, we adapted the same procedure of Kmeans used by Allen et al. ( Allen et al., 2014 ): L1 (Manhattan) distance is used, as it was suggested to be more effective than L2 (Euclidiean) distance for high-dimensional data ( Aggarwal et al., 2001 ). The algorithm is replicated 100 times to increase chances of escaping local minima, and centroid positions were randomly initialized. Then, Kmeans returns the solution with the lowest 'SUMD' (within-cluster Sums of points-tocentroids Distances).

Comparative analysis 2.2.5.1. MEG group-level analysis.

2.2.5.1.1. Selection of optimal number of components (NC opt ). In the context of dimensionality reduction methods, the choice of the optimal number of components (NC opt ) to be extracted is still a challenging issue. Here, we used the well-known approach: 'Elbow criterion' ( Allen et al., 2014 ) for Kmeans method (with maximum number of clusters = 10). For all other SS methods, we estimated NC opt based on the goodness of fit approach ( Timmerman and Kiers, 2000 ;Wang et al., 2018 ) previously used by many recent works ( Tewarie et al., 2019b ;Zhu et al., 2020 ). We performed the 'DIFFIT' method that refers to the difference in data fitting with a range of input NC varied from 2 (for motor tasks) and 4 (for working memory task) to 10 components, and selected the NC that gives the largest DIFFIT value as the NC opt . Technical details about these approaches can be found in the Supplementary Materials.

2.2.5.1.2. Selection of significant components. Among the NC opt extracted components, identifying those that reflect genuine brain activity related to the task is critical. In this paper, we followed a testing procedure adopted by ( O'Neill et al., 2017 ) and previously described in ( Hunt et al., 2012 ;Winkler et al., 2014 ) to determine significant components modulated by the tasks. The testing relies on the construction of empirical null distribution based on a 'sign flipping' permutation approach. Therefore, a component was considered significant if, at any time point, the corresponding time signal, averaged over trials, fell outside a threshold defined at 0.05 with corrections. For all datasets, 2tailed distribution was allowed, and Bonferroni corrections were applied for multiple comparisons across the NCopt components and across tem-poral degree of freedom. More details about 'sign-flipping' approach and threshold values setting can be found in Supplementary Materials.

MEG subject-level analysis.

Besides group-level analysis, it is crucial to test the performance of each method when applied directly on individual dFC. To this end, instead of concatenating trials from all subjects as in the final step of 'G-dFC' computation, we perform, for each subject, a dFC concatenation of all trials related only to this subject to form a subject specific dFC, denoted 'S-dFC'. Then, all selected SS methods were applied on 'S-dFC' matrix to extract subject-specific spatial and temporal signatures (k = 10). In order to quantitatively evaluate and compare methods strength at subject-level context, we measure, for each method, both spatial and temporal similarities between each extracted S-dFC component and significant G-dFC components. These parameters are:

(1) Average Distance (AD) for 'spatial similarity':

𝐴𝐷 = ∑ 𝑠 𝑑 ( 𝑛 𝑠 , 𝑛 𝑔 ) 𝑁 𝑠 𝑠 ∈ [ 1 , 𝑁 𝑠 ] ; 𝑔 ∈ [ 1 , 𝑁 𝐺 ] (2)
Where 𝑑( 𝑛 𝑠 , 𝑛 𝑔 ) is the Euclidian distance between the node 𝑛 𝑠 of S-dFC network and the nearest node 𝑛 𝑔 from the significant G-dFC network. 𝑁 𝑠 is the total number of nodes in S-dFC network, and 𝑁 𝐺 denotes the total number of nodes in G-dFC network. All networks were 70% thresholded. Lower values of AD indicate stronger spatial similarity between S-dFC and G-dFC networks.

(2) Correlation Signals (CS) for 'Temporal Similarity':

𝐶𝑆 ( 𝑇 𝑆, 𝑇 𝐺 ) = ∑ 𝑠 ∑ 𝑔 ( 𝑇 𝑆 𝑠𝑔 -𝑇 𝑆 )( 𝑇 𝐺 𝑠𝑔 -𝑇 𝐺 ) √ ( ∑ 𝑠 ∑ 𝑔 ( 𝑇 𝑆 𝑠𝑔 -𝑇 𝑆 ) 2 )( ∑ 𝑠 ∑ 𝑔 ( 𝑇 𝐺 𝑠𝑔 -𝑇 𝐺 ) 2 ) 𝑠 ∈ [ 1 , 𝐿 𝑠 ] ; 𝑔 ∈ [ 1 , 𝐿 𝐺 ] (3)
Where 𝑇 𝑆 is the temporal signal of each S-dFC component of length 𝐿 𝑠 and 𝑇 𝐺 represents temporal signals of G-dFC significant component of length 𝐿 𝐺 . Higher values of CS reveal stronger temporal similarity between S-dFC and G-dFC signals.

We perform this analysis on each subject among the 15 subjects of the MEG dataset 1 (Motor task). Therefore, for each method, we counted the number of subjects that show satisfactory results performance in the context of S-dFC, based on the previously explained measures. Then, to approximate the number of subjects/trials needed for each SS method to give significant results, we follow the same procedure explained above, but instead of single subject S-dFC computation, we increased the number of concatenated subjects in dFC computation from N subj = 2 to 14, progressively. In order to have generalized and reliable results, we considered all possible combinations relative to each N subj ( 𝐶 15

𝑁 𝑠𝑢𝑏𝑗

), where different sets of N subj subjects were selected among the 15 existing data subjects.

Results

In the following, we present our evaluation study on real MEG data, however our methodology was also tested on simulated data. These results can be found in the supplementary material. Briefly, the simulation-based analysis showed that all methods provide satisfactory results in terms of spatial and temporal similarity between reconstructed and simulated components with the best performance for NMF method and the worst for SOBI. All methods, except for FastICA, NMF and Kmeans, provided consistent results. PSAUD and PCA were the fastest. Results revealed that SOBI, NMF and Kmeans converge more slowly than others with the increased value of SNR. Reader can refer to supplementary material to see the detailed quantitative analysis on simulated data.

We firstly ran each algorithm at each value of NC and calculated the corresponding DIFFIT values in order to select the optimal number of components (NC opt ) relative to these SS methods. Results are shown in Fig. 2 for the three empirical datasets. Hereinafter, we set NC to the computed NC opt value relative to each method and task.

Results of different SS methods applied on empirical data are illustrated in Figs. 345. In each Fig., we presented only the components that demonstrated significant task modulation based on the applied null distribution approach (described in the methods section). The networks were thresholded only for visualization purpose (70% for dataset 1 and 3, 85% for dataset 2). Corresponding dynamic reconfiguration of each significant network were plotted together. The temporal fluctuations represent component time signals averaged over trials and subjects.

Self-paced button press task

In this task, participants were asked to press a button with the index of their non-dominant hand every 30 seconds.

Based on literature findings (see table S1 in supplementary materials), we were interested in quantifying SS methods ability to extract a sensorimotor network from significant components. To this end, we defined a brain network with activated AAL regions in both motor cortex (including precentral, paracentral, rolandic and supplementary motor areas) and somatosensory cortex (including postcentral, parietal and supramarginal areas) serving as a mask template for our network of interest (sensorimotor network), illustrated in Fig. 7 . Then, we selected each significant network and computed the strength of each activated node in that significant network (defined as the sum of all edges weights connected to that node). The ratio of the strength of activated AAL nodes that belongs to sensorimotor mask is calculated relative to the strength of all activated nodes in that significant network. In case the ratio is greater than a certain threshold value, the network is considered as a sensorimotor network denoted as 'mot' in the Fig. 3 . Otherwise, the network is denoted 'Aux' referring to auxiliary network. After many trials (threshold = 0.5, 0.6, 0.7), threshold value was set to 0.6 as it has shown more convenient results, when visually inspecting components classification (false positive and false negative). The reader can refer to supplementary materials (Fig. S8) for more details about the computation of the ratio values for all components.

Fig. 3 shows that all SS methods were able to extract at least one significant 'Mot' network. All significant components extracted from the five ICA methods (JADE, InfoMax, FastICA, CoM2 and PSAUD), NMF and Kmeans methods were categorized as 'Mot' network due to the strong participation of sensorimotor nodes in these networks (sensorimotor strength ratio > 0.6) although some of them may involve additional few connections to other regions. On the other hand, SOBI and PCA methods showed 'Aux' networks (sensorimotor strength ratio < 0.6) besides 'Mot' networks, with remarkable activations in frontal regions. Temporal variation was similar for almost all significant components over all methods showing a peak value at 0sec, the button press time, with slight differences in amplitude values, indicating signal intensities relative to each component. Note that negative connectivity, referred to as blue connections in spatial networks and negative temporal values in temporal signals, represents desynchronization between brain regions. Therefore, all studied SS methods were able to extract at least one significant component that highlight strong connections between sensory and motor regions modulated significantly by the task at the exact button press instant ('Mot').

Left-hand movement task

This task is also motor but different than the previous one. Here the participants were asked to rapidly and successively tap their left index and thumb fingers. Similarly to the previous task, the same sensorimotor mask was applied to quantify resultant networks to discriminate 'Mot' from 'Aux' networks. The reader can refer to supplementary materials (Fig. S9) for more details about the computation of the ratio values for all components. . DIFFIT values are plotted against number of component 'J' for all ICA methods, PCA and NMF. The blue plot corresponds to the self-paced button press task, the orange plot for HCP left-hand movement task and the yellow one for Working Memory task. The optimal NC that gives the highest value of DIFIT relative to each task is marked by a small circle on the x-axis. Results of optimal NC relative to Kmeans using the elbow criterion is also shown. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) Fig. 4 shows that not all SS methods were able to extract a sensorimotor 'Mot' network. For example, none of the significant components of SOBI, PSAUD and Kmeans has survived the threshold imposed for the strength of activated sensorimotor nodes (Fig. S9) and are therefore considered as 'Aux' networks as indicated in Fig. 4 . In these methods, 'Aux' networks consist of either a visual network significantly modulated directly after the onset of movement in Kmeans and 0.35sec before onset in SOBI and PSAUD, or a network (one in SOBI and PSAUD and six in Kmeans) involving strong connections between almost all brain areas modulated at 0.2sec before and after onset. It should be noted that although this network shows strong activation of the right precentral and postcentral nodes, it failed to be quantified as a sensorimotor network due to the high coverage of the brain.

All remaining SS methods were able to extract one 'Mot' network among all significant components. The spatial representation of 'mot' network involves sensorimotor with some cingulate nodes from the left cortex in JADE, right cortex in InfoMax and both left and right cortices in other SS methods. These 'Mot' networks show significant drop in connectivity around 0.2sec following the movement onset. Significant increased modulation was also observed at -0.2sec in InfoMax, FastICA, CoM2, PCA and NMF methods. Exact times of components significance are indicated with stars on Fig. 4 . From this Fig., we can see that 'Aux' networks show various spatial patterns between SS methods (such as the integration of areas from visual, motor-frontal, motor-visual, temporal lobes…).

Regarding temporal evolution, the hand movement here are much more frequent than the previous task. Clearly, the fast neural activity due to the short time between successive button presses is expressed as an oscillatory behaviour of brain network activity around the zero time button press, as showed also previously ( Vidaurre et al., 2018a ). This yields to the obtained temporal variation where the motor network state seems to have high connectivity before button press and begin to have a drop-in connectivity to reach its significant peak after ~0.2sec (referred to as a desynchronization in high frequencies ( Vidaurre et al., 2018a )).

Working memory task

This task is much more complex comparing to the other two tasks. Subjects here were asked to visualize and memorize two visual shapes and respond to a third probe stimulus by a button press (with their right index finger) in case of matching. The increased cognitive load evoked by the Sternberg task is expected to induce variations in a greater number of significant brain networks including stimulus visualisation (visual network), semantic processing and pattern recognition (semantic, language networks) and button press response (sensorimotor network).

In a similar context of previous tasks, four masks were defined here related to the most relevant working memory related networks found in literature. These masks are also illustrated in Fig. 7 . The visual mask consists of the activation of primary visual cortex (occipital areas, cuneus, calcarine and lingual) and is denoted as 'Vis' in the Fig. 5 . The semantic mask involves connections between bilateral temporal (including fusiform, heschl, parahippocampal) and parietal lobe (postcentral, supramarginal, angular, precuneus). The language mask (denoted 'Lang') is defined as a left lateralised network with activation of nodes from temporal, frontal and parietal regions from the left cortex. The sensorimotor mask ('Mot') is previously defined in motor tasks. Detailed ratio values of the four masks for all components can be found in details in supplementary materials (Fig. S10). In this task, by applying Fig. 3. Self-paced motor task results. Spatial and temporal distribution of all significant components derived from all compared SS methods applied on G-dFC in the self-paced motor task (N = 15 subjects). All brain networks were thresholded for visualization; lines width indicates connectivity strength between regions. Red lines represent positive connectivity values while the blue ones represent negative values. Integrated AAL nodes are represented by spheres of different sizes that reveal connectivity weights (strength) between that region and the rest of brain. Corresponding temporal evolution is averaged across all trials and subjects. Time values on the x-axis represent the position of the sliding window's center, relative to the button press at t = 0sec (as illustrated by a vertical line). A color code is attributed for each component in space and time. For each SS method, only significant components (p corrected < 0.05) that appear outside the 'sign-flip' based null distribution (as described in methods sections) are shown here. All NC opt extracted components with corresponding null distribution are shown in Supplementary Fig. S4 for an example of ICA-JADE method. Note that sensorimotor network is clearly activated at the button press instant in all SS methods. In this Fig., 'Mot' refers to sensorimotor network and 'Aux' refers to all others 'non-sensorimotor' networks. An interactive version of ICA-JADE results can be found on our github https://github.com/judytabbal/dynbrainSS.git using rotatable MATLAB figures. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) simultaneously four masks on the same network component, there is a possibility that the ratio strength of more than one mask survives the threshold (0.6). In this case, the network belongs to the mask that gives the highest ratio strength value. In case of equality between two masks, we consider that the network belongs to both masks.

Fig. 5 illustrates the spatial distribution of all significant components for all SS methods, and temporal variation for three of these methods (JADE, NMF and Kmeans) for visualisation clarity. The temporal evolution of other SS methods can be found in Supplementary Fig. S7.

Starting from t = 0sec, two visual stimulus (shapes) were presented successively, each for 0.6sec. During this period, all methods, except for PCA, were able to extract one or more significant 'Vis' network. We can notice few additional connections from occipital to parietal or temporal regions in this network. Time variation of this network shows significant peak during the first two seconds period.

Following stimulus presentation, subjects should retain the observed shapes in working memory. During this period, known as the maintenance phase, all methods, except Kmeans, show significant decreased modulation at [4-6sec] of at least one 'Sem' network. The breakdown in Fig. 4. HCP left hand movements task results. Spatial and temporal distribution of all significant components derived from all compared SS methods applied on G-dFC in the left-hand motor task (N2 = 61subjects). All brain networks are thresholded for visualization. Time values on the x-axis represent the position of the sliding window's center, relative to the button press at t = 0sec (as illustrated by a vertical line). A color code is also attributed for each component in space and time. For each SS method, only significant components (pcorrected < 0.05) that appear outside the null distribution are shown here. Exact times of significance relative to each component are indicated with stars., revealing an oscillatory temporal activation of motor component for all SS methods. All NC opt extracted components with corresponding null distribution are shown in Supplementary Fig. S5 for an example of ICA-JADE method. In this Fig., 'Mot' refers to sensorimotor network and 'Aux' refers to all others 'non-sensorimotor' networks. An interactive version of ICA-JADE results can be found on our github https://github.com/judytabbal/dynbrainSS.git using rotatable MATLAB Figs.. Reproducing these results is also possible/available using the MATLAB interface on github. this network's connectivity was previously demonstrated ( O'Neill et al., 2017 ). During the same period, we can notice a drop-in connectivity in 'Sens' network revealed only by NMF method.

At [10-12sec] period, many networks seem to be significantly modulated among all methods: (1) 'Vis' is re-activated at the probe stimulus presentation in all methods except for PSAUD, (2) 'Sens' shows significant increase with FastICA method. This network becomes most strongly connected around the time button press response, (3) 'Lang' is commonly derived by JADE, InfoMax, and PCA, and exhibits an increased connectivity peaking during probe presentation. We can notice that this network is also significantly decreased in CoM2, PSAUD and Kmeans around 5sec. (4) Three methods (JADE, InfoMax and PCA) also showed significant increased modulation of 'Sem' network. (5) A network that belongs equally to both 'Sem' and 'Sens' masks (denoted as 'Sem + Sens' network) is strongly activated during this period in CoM2 and PSAUD. Note that these two masks have common brain regions mainly in parietal lobe responsible for sensory processing which is coherent with the task evolution. This 'Sem + Sens' network is also modulated in PCA method. Besides these components, few 'Aux' networks were also considered significant as shown in Fig. 5 .

In summary, the three SS methods (CoM2,PSAUD and PCA) succeeded to derive all expected components with their appropriate temporal significant modulation. JADE and InfoMax were able to extract visual, semantic and language but not the sensorimotor network. FastICA and NMF missed the language component. However, SOBI was unable to show both sensorimotor and language networks and Kmeans failed to extract semantic and sensorimotor components. For the three tasks, spatial and temporal distribution of the NC opt components derived from JADE method, with corresponding null distribution, can be found in Supplementary Figs. S3, S4, S5. For further clarity in visualisation and interpretation, we illustrated, in Fig. 6 , the spatiotemporal reconfiguration of the functional brain networks as obtained by ICA-JADE.

Furthermore, we discriminated different SS methods performance in terms of the activation of relevant brain networks in each task. For example, in motor tasks, we calculated the 'Mot' network occupancy percentage defined as the number of 'Mot' networks (quantitatively defined by the mask as explained above) divided by the total number of significant components found in the corresponding SS method. Similarly, for working memory task, the occupancy percentage of visual, semantic, language, sensorimotor and auxiliary networks were evaluated for each method. Results are shown in Fig. 7 with the spatial representation of the corresponding relevant brain regions. Therefore, Fig. 7 resumes the overall performance of each SS method showing variability in the methods' ability to directly extract the appropriate task-related components.

Performance of each SS methods at subject-level

Here our objective is to evaluate the performance of the methods at the subject-level. We test i) the capacity of each method to extract significant components related to the task: to do so, we computed the correlation between the components obtained by each method on each subject with the significant network obtained at the group level and ii) the number of subjects needed for each method to detect 'expected' networks: here we tested the overall performance of each method by increasing the number of subjects, going from 1 to 15 as we performed subject-level analysis on the self-paced data. Fig. 8. A summarizes the subject-level analysis scenario. For each method, we chose one of the significant motor components derived from the decomposition of the group-level (N = 15subjects), mostly the one showing little intervention from regions other than sensorimotor ('Mot') and having high temporal coefficients amplitude (supposed to be the best for each method). This component illustrates a 'group' motor network with temporal modulation at the button press time. It will, eventually, serve as a 'mas' component for subject-level analysis, as we are concerned in motor component extraction.

For each SS method, NC = 10 components were derived from each subject data. Then, Average Distance (AD) and Correlation Signals (CS) between each of these components and the 'group' motor component relative to the SS method were computed in order to quantify the ability of the method to extract, from a single subject, a task-related component in space (motor network) and time (temporal modulation at button press time) respectively. Following this, only one of these 10 components is selected for results calculation. This selection is based on two conditions criteria on AD and CS values. In the case where AD component is less than a threshold (set as the average of AD values of all components for all subjects and SS methods), and CS is higher than 0.7 (chosen as a trade-off between moderate and high correlation), then the component is considered to be a motor component. By setting these thresholds, we considered the existence of inter-subject differences, thus, allowing subjects to have different but near spatial distribution of motor network. Therefore, if at least one of the extracted components pass these conditions, corresponding AD and CS values are denoted and the number of subjects that give similar results to group-level is raised by 1. Otherwise, we selected component 𝑖 as the nearest component to the group-level result, with a compromise between spatial and temporal similarities.

A typical example is illustrated in Fig. 8. A showing that PCA decomposition was able to extract a motor component from subject 2 (component 5), whereas no motor component was derived from subject 14, AD and CS values of component 7 were denoted in this case. Spatial and Fig. 7. SS methods performance evaluation for real MEG tasks. For each task, brain regions involved in each relevant network are illustrated on the left side using the AAL atlas, while brain networks occupancy are shown on the right side. The occupancy percentage represents the presence percentage of these defined brain networks relative to all significant extracted components. For motor tasks, motor network ('Mot') was emphasized with auxiliary ('Aux') networks relative to all significant components, while visual ('Vis'), semantic ('Sem'), language ('Lang'), sensorimotor ('Sens') and 'Sem + Sens' networks are highlighted in contrast to other auxiliary ('Aux') networks. Referring to these representations, capabilities of different SS methods in extracting relevant task-related components can be evaluated. temporal distributions of selected subjects' components in both cases are shown in Fig. 8 . Results of the remaining components for these 2 subjects' examples are shown in Supplementary Fig. S11.

As a result, two parameters were collected and represented in Fig. 8. B and 8. C respectively. Group-subject similarity percentage was calculated as the number of subjects that gives a motor component similar to the group-level result relative to the total number of subjects (N = 15). Fig. 8. B illustrates this parameter for all SS methods. We can see that JADE was able to extract a task-related component from 8 out of 15 subjects (53.33%), InfoMax and PSAUD from 7 subjects (46.67%), SOBI and CoM2 from 6 (40%), FastICA from 5 (33%), PCA and NMF from 4 (26.67%) and Kmeans from 3 (20%). The Fig. 8. C shows the distributions of AD and CS values of selected components from each subject over all SS methods. Methods with higher subject-group similarity percentage have lower median values of AD and higher median values of CS. In addition, we can notice from AD and CS median values that similarity in space was much easier to be satisfied than temporal similarity for most SS methods. Interquartile range values show the existence of inter-subject variability results. However, some methods showed higher interquartile range of AD values (CoM2 and PCA), or CS values (JADE, CoM2 and NMF) relative to other methods.

The optimal number of subjects of each SS method

Then, the same procedure was applied with increasing the number of subjects from one subject (single-subject) to 14 subjects. AD and CS are Group motor component is shown for PCA example. AD and CS values computed for all components are presented for two subjects. Based on AD and CS condition limits highlighted in each polar bar, success and failure in extracting motor component are both illustrated by subjects 2 and 14 respectively. Spatial and temporal distribution with corresponding AD and CS values for all components of both subjects 2 and 14 are illustrated in Supplementary Fig. S11. A small table on the right illustrates results of success and failure for all subjects in PCA. B. Results of the number of subjects that successfully extracted a motor component in each SS method relative to the total subject's number, denoted Group-Subject Similarity, are shown. C. distributions of AD and CS values of all selected components for the 15 subjects are illustrated. D. Generalization study with increasing number of subjects and the corresponding results of Group-Subject similarity percentage (when considering all possible combinations). E. Grey shaded values (5,6,10 and 14) represent the critical number of subjects required for each SS method to have a 90% precision in extracting the task related component. computed for all possible combinations. The number of possible combinations is calculated. For example, 7-subjects analysis requires 6435 combinations, hence 6435 values of AD and CS. For each number of subjects, we calculated subject to group similarity as the ratio between number of combinations that succeeded in extracting a motor component relative to the total number of possible combinations. Results are illustrated in Fig. 8. E, F. As expected, the percentage similarity increases with increasing number of subjects for all SS methods. Fluctuations in similarity results are observed in some methods (as FastICA, NMF and Kmeans) due to the non-consistency characteristic of these methods (as previously proved). Some methods required a smaller number of subjects for data analysis to provide satisfactory results (motor component at button press time in our case) than others. For example, the four ICA versions (JADE, InfoMax, SOBI and PSAUD) required 5 subjects in order to attain a minimum similarity level of 90% between subject and group level results. CoM2 and PCA required 6 subjects, while much more subjects were needed for others (10 subjects for FastICA and Kmeans and 14 for NMF) as showed in Fig. 8. E. Overall, ICA methods and specially those based on the high order statistics (such as JADE) outperform other methods in extracting networks at the subject-level.

Discussion

In this study, we have evaluated the robustness of the most popular SS methods applied to extract the main brain networks fluctuating during time in order to help researchers make a rational choice (if any) among the multitude of available methods. Specifically, nine algorithms have been compared using simulated data (see supplementary materials) and three independent MEG datasets (N = 95) recorded during motor and memory tasks. The discrepancy in the datasets size and behavioral tasks performed allows testing SS methods performance on different scenarios. As the evoked responses (analyzed here) last for hundreds of milliseconds, we conducted our comparative analysis based on MEG datasets to benefit from the excellent temporal resolution of this technique. However, the same pipeline study can be applied in task-related fMRI context.

Overall, our results show variability between the evaluated SS methods and even between ICA subtypes. The performance of these methods depends on the nature of the task (simple vs complex, slow vs fast time scale tasks). In a simple and relatively slow time scale task (as self-paced button press task), all methods succeeded in tracking spatially and temporally the dynamic brain activity. However, when it comes to much faster task (HCP motor task) or more complex task (Working Memory), some methods (SOBI and Kmeans for instance) showed lower performance in extracting relevant brain networks (as defined by our masks). Results relative to each task will be discussed later in details.

First, the quantitative comparison performed on simulated dynamic networks showed that all SS methods have successfully separated functional networks based on their connectivity time courses. However, spatial and temporal similarities in SOBI were significantly lower than other SS methods, especially for the fourth simulated state (P4) and the second state (P2) as shown in Supplementary Fig. S3, which involves more complex spatiotemporal activity than other states. As expected, FastICA, NMF and Kmeans methods were proved inconsistent with multiple runs. This is caused by the nature of these algorithms that is based on random input initialisations until solution convergence. The noise effect on the results obtained was also tested and showed an increased performance for all methods with higher SNR value with slower convergence to the optimal accuracy for some methods (SOBI, NMF and Kmeans) relative to others. Regarding computation time, CoM2, PCA and PSAUD were the fastest whereas InfoMax and JADE were the slowest. Still, the executed time of these algorithms is sensible to dataset's features (size, complexity, type…). Other metrics such as the number of floating-point operations (FLOPs) required for the algorithm completion could be also tested. We also suggest for future studies to explore other data simulation approaches that build the desired ground-truth brain states based on more realistic modeling (using Neural Mass for instance), however, this may introduce the effect of other parameters in the comparison (forward problem, inverse solution…). Besides simulation approach, some studies attempt to consider fMRI data as a ground truth to quantify and compare SS methods performance in the context of M/EEG studies [START_REF] Colclough | How reliable are MEG resting-state connectivity metrics?[END_REF][START_REF]Exploring the correlation between M/EEG source-space and fMRI networks at rest[END_REF].

The method's performance were evaluated on three real MEG datasets already published and tested by previous studies ( Casorso et al., 2019 ;O'Neill et al., 2017 ;Tewarie et al., 2019a ;Vidaurre et al., 2018a ;Zhu et al., 2020 ). According to self-paced motor task, results showed that all SS methods have successfully extracted one or more significant network that involve strong connectivity between sensorimotor regions ('Mot'). For HCP data, a similar sensorimotor network was revealed among significant components in all SS methods except for SOBI, PSAUD and Kmeans. Integrated regions in this network mainly include nodes from central and parietal gyrus. The sensorimotor network is strongly coherent with the task [START_REF] Melnik | EEG correlates of sensorimotor processing: independent components involved in sensory and motor processing[END_REF]Yousry, 1997 ) since it requires both movement (through button press or hand movement) and tactile response (as the subject will feel the button or fingers tape). The effect of right-handedness of all participants of self-paced dataset is also revealed by the presence stronger implication of sensorimotor nodes from the right cortex relative to the left one as revealed by the sphere sizes and connections in Fig. 3 . It is noteworthy to mention the existence of a network that highlighted significant connections in the visual lobe in JADE for self-paced button press task and most SS methods for HCP task. This network was previously noticed by Oneill et al. studying the same button press task ( O'Neill et al., 2017 ). This can be interpreted as a cross modal synchronization between visual and sensorimotor cortex as previously studied ( Bauer et al., 2020 ). Regarding temporal evolution, it is clear that all networks modulate significantly with the exact button press time for self-paced task. Differently, the temporal variation related to HCP motor task takes an oscillatory shape, which was also reported by other studies dealing with the same dataset ( Vidaurre et al., 2018a ;Zhu et al., 2020 ). A possible reason for this activity was suggested by ( Vidaurre et al., 2018a ) considering a leakage effect of temporal activity of previous button press into the next trial due to the fast successive trials. In both tasks, there exists auxiliary 'Aux' networks that significantly modulated with the task but not directly related to the motor cortex activity. The occupancy percentage of 'Aux' networks increases for all SS methods in HCP results mainly in CoM2 and Kmeans. The presence of these networks can be related either to the robustness/sensitivity of the SS method relative to spurious networks or to the reliability of the techniques used for selection of optimal number of components (DIFFIT) or significant components (null distribution) that will be further discussed later.

In order to evaluate the spatial and temporal accuracy of SS methods at higher levels of complexity, we tested the methods on Sternberg working memory task. All SS methods detected visual network, which is consistent with the presentation of visual stimuli at two different times. Regions in the primary visual lobe related to stimulus visualisation [START_REF] Grill-Spector | Cue-invariant activation in object-related areas of the human occipital lobe[END_REF] and lateral occipital cortex responsible for object/shape recognition ( Corbetta et al., 1991 ;[START_REF] Grill-Spector | The lateral occipital complex and its role in object recognition[END_REF][START_REF] Kourtzi | Representation of perceived object shape by the human lateral occipital complex[END_REF] were present in these networks. The button press response is reflected by sensorimotor connections consistently with previous working memory studies [START_REF] Metzak | Constrained principal component analysis reveals functionally connected load-dependent networks involved in multiple stages of working memory[END_REF][START_REF] Yamashita | Predicting learning plateau of working memory from whole-brain intrinsic network connectivity patterns[END_REF] only by using CoM2, PSAUD and PCA methods. In order to process and maintain observed stimuli as a way to memorize them, a higher level of cognition is illustrated by a 'semantic network', which mainly encompasses bilateral parietal and temporal areas activation in all SS methods, except for Kmeans. This is coherent with previous studies that demonstrate the evident role of parietal cortex as a workspace for sensory and perceptual processing in working memory framework ( Chai et al., 2018 ) through angular [START_REF] Frackowiak | The anatomy of phonological and semantic processing in normal subjects[END_REF]Vandenberghe et al., 1996 ), precuneus ( Cavanna and Trimble, 2006 ), and hippocampal ( Baddeley et al., 2011 ) areas. Bilateral inferior temporal regions also play important role in semantic processing [START_REF] Nestor | Declarative memory impairments in Alzheimer's disease and semantic dementia[END_REF]Vigneau et al., 2006 ). Fusiform gyri, strongly modulated in our results, has also shown a particular concern in this context [START_REF] Mion | What the left and right anterior fusiform gyri tell us about semantic memory[END_REF]. The detection of the 'language' network by JADE, InfoMax, CoM2, PSAUD, PCA and Kmeans methods, was compatible with previous findings ( Brookes et al., 2011b ;O'Neill et al., 2017 ). Temporal and parietal lobes were remarkably activated by these methods, mainly the parahippocampal and supramarginal gyri respectively. These regions are critical in memory encoding and retrieval and semantic cognition ( Axmacher et al., 2008 ;Caminiti et al., 2015 ;Demb et al., 1995 ;[START_REF] Derrfuss | Cognitive control in the posterior frontolateral cortex: evidence from common activations in task coordination, interference control, and working memory[END_REF][START_REF] Deschamps | On the role of the supramarginal gyrus in phonological processing and verbal working memory: evidence from rTMS studies[END_REF]Vigneau et al., 2006 ). In a similar (abstract shape based) working memory task, the interpretation of this network was related to a verbalisation naming strategy employed by participants as a way to aid in memory encoding ( Caminiti et al., 2015 ;O'Neill et al., 2017 ). Therefore, this network's activation may be possible with the task as it modulates strongly with the probe presentation and response time.

It is important to point here that the resultant networks were denoted objectively in this study using a quantification approach. For a better interpretation of the functional significance of results, we built template brain masks, referring to the literature (see table S1 in supplementary materials), from AAL cortical regions. These masks are used for seeding our networks of interest: 'Mot' in motor tasks, 'Vis'/'Sem'/'Lang'/'Sens' in working memory task. Networks were then classified based on their activated nodes strength relative to each mask. In addition, the used masks have distinct spatial distribution, with some shared regions mainly between 'Sem' and 'Sens' networks. However, when we aim to deeply investigate task-related sub-networks, the 'mask' technique seems to have much more complexity related to the specificity of the integrated brain regions and the precision of an accurate threshold in order to appropriately classify networks results.

Although we performed our study on cognitive tasks, it is a topic of great interest to apply this methodology pipeline on resting-state experiments since many studies have shown the dynamic reconfiguration of the brain during rest as well ( Kabbara et al., 2017 ;[START_REF] Liégeois | Resting brain dynamics at different timescales capture distinct aspects of human behavior[END_REF].

Regarding methodological considerations, first, the optimal number of components to be derived was still a challenging question for all SS methods rather than a direct limitation of our algorithms. In this study, we applied two well-known algorithms (DIFFIT and elbow criterion) for a range of number of components to automatically select the optimal number of components relative to each SS method. The evaluated range of NC values was upper limited by 10 components in order to avoid spurious networks. Moreover, we tried to fix the number of components to 10 for all SS methods in the self-paced motor task at the group-level, as already set in a previous work ( O'Neill et al., 2017 ) that used FastICA algorithm and little difference was observed for the overall results.

Second, we should point that not all these components are necessarily essential, especially in the case of simple tasks as motor tasks. To this end, we followed the approach of the null distribution based on sign flipping algorithm ( Hunt et al., 2012 ;O'Neill et al., 2017 ;Winkler et al., 2014 ;Zhu et al., 2020 ) to select only components whose temporal dynamics significantly modulate with the task. In this way, we ensure that brain dynamics relative to the studied behavioural tasks can be summarized and described by the retained components through an automatic way that allows us to objectively compare the SS methods performance. In addition, the fact that this technique is a purely data-driven procedure that does not require any prior hypothesis or conditions manipulation makes it likely adapted to the specific examined dataset. Two points regarding significant components selection are important to mention here. First, in the applied null distribution, networks were defined to be significant if they fell outside the null distribution in either positive or negative sides because they reflect trial-onset-locked that either increases or decreases in connectivity across subjects (as amplitude envelope correlation was adopted). Second, it should be noted that components significance was evaluated relative to the specific task duration. For example, temporal duration of the entire analysis for self-paced and working memory tasks can include dynamics that should be excluded from the analysis. In this context, we limited our significance interpretation/assessment in the interval of [-2; + 2sec] and [-0.5; + 0.5sec] relative to the button press instant in the case of self-paced and HCP motor tasks respectively, and [-2; + 16sec] relative to the visual stimulus presentation for memory task. In addition, few limitations are to be discussed when dealing with this selection. First, it was not convenient to rely on this technique when the number of trials and subjects was either too small or too big. A small number will not allow to build a reliable null distribution while a huge one will have its computational cost regarding all possible subjects' combinations for sign-flipping procedure, as already executed in the HCP analysis. Moreover, there is no consensus about the thresholds/margins that define well a limit level for component's amplitude. For instance, there exists networks whose temporal variation peaks at the limit of null distribution envelope. These are considered to be critical components that may be integrated in the task but considered not to be following the automatic criteria of this null distribution. Future works should therefore investigate more about this methodological approach in the framework of cognitive tasks, in addition to resting state experiments. Also, it is crucial to seek more methods to use or combine with the applied technique in order to have more robust basis for significant components selection. It is noteworthy to report that null distribution-based technique was applied uniformly for all SS methods, thus our main objective of comparison was built on a unified evaluation framework.

Third, we used the same pipeline supported by the previous studies dealing with the same dataset (cortical parcellation, source reconstruction, functional connectivity metric and source leakage correction, frequency bands and sliding window settings) ( Kabbara et al., 2019 ;O'Neill et al., 2017 ). By applying already tested and validated methodological approaches, we avoid influencing factors on the comparison performed. However, we point out that other methodological solutions could be exploited by other researches using the same pipeline adopted in this work. Regarding cortical parcellation, we chose AAL atlas based on its successful use in previous MEG investigations ( O'Neill et al., 2017 ;[START_REF] Tewarie | Predicting haemodynamic networks using electrophysiology: the role of non-linear and cross-frequency interactions[END_REF]. This atlas also provides good basis for the orthogonalisation procedure adopted since its number of regions is sufficiently low (78 ROIs) and well separated [START_REF] Colclough | A symmetric multivariate leakage correction for MEG connectomes[END_REF]. The beamformer spatial filtering was selected as the inverse problem solution due to its demonstrated efficiency in the measurement of static ( Brookes et al., 2011a ) and dynamic ( Baker et al., 2014 ) functional connectivity. Functional connectivity between ROIs regions was estimated through Amplitude Envelope Correlation (AEC). This technique has been successful in elucidating electrophysiological networks of functional connectivity [START_REF] Colclough | How reliable are MEG resting-state connectivity metrics?[END_REF]. Other methods, such as phase couplings can be considered as an alternative way to probe different type of functional connectivity ( Lachaux et al., 1999 ). Sliding window settings (length and step) were selected carefully as a trade-off between temporal resolution and the accuracy of the derived adjacency matrices (length = 6sec, step = 0.5sec for self-paced and working memory tasks) ( O'Neill et al., 2017 ). However, according to recent works ( Fraschini et al., 2016 ;Liuzzi et al., 2019 ), it can be seen that metrics (including amplitude envelope correlation) perform poorly for very short state durations when combined with the sliding window approach below few seconds, providing noisy results of low correlation with ground truth in simulations. For this reason, we followed the work of ( Tewarie et al., 2019b ) to estimate dynamic functional connectivity by taking sample by sample time series rather than windowed aggregated samples using the Instantaneous Amplitude Correlation (IAC). This high temporal resolution measure of FC has shown great sensitivity to genuine fluctuations in functional connectivity applied in the same context of our study.

It would be interesting to test data-driven windows approach in this context using the recurrence plots of the amplitude envelopes as in ( Tewarie et al., 2019b ) instead of averaging trials in dataset 2 to reduce heavy dFC matrices. It could be also tested against fixed sliding window approach as for datasets 1 and 3.

Concerning frequency bands, it was crucial to preprocess each dataset in its appropriate bandwidth. For example, brain signals in selfpaced and HCP motor tasks were proved to be more active in the beta band, while broader range of frequency bands are integrated in complex cognitive tasks as working memory ( O'Neill et al., 2017 ;Zhu et al., 2020 ).

Conclusion

Deciphering of dynamics of electrophysiological brain networks is one of the most important goals in neuroscience. In this paper, we evaluated and compared nine popular source separation (SS) methods to identify dominant networks of connections with corresponding temporal dynamics at group-level as well as subject-level, using simulation and empirical MEG data (N = 95subjects) recorded during three different tasks: (1) simple button press task, (2) fast finger movement task (HCP) and (3) Sternberg working memory task. Results show close consistency for all SS methods in successfully identifying a transient network of connections linking somatosensory and primary motor regions in the relatively slow and simple button press task. Variability between these methods' performance is revealed in rapid tasks of sub-second timescale (HCP motor task) and in a more complex task (Sternberg). The SOBI and Kmeans algorithms showed the weakest performance among tested methods. CoM2, PSAUD and PCA showed promising results in working memory task, revealing the formation and dissolution of multiple networks that relate to semantic processing, pattern recognition and language as well as vision and movement. At the subject level analysis, ICA methods using high statistical order (JADE, InfoMax, CoM2 and PSAUD) outperform other methods. Our main message is that researchers should be aware to select the appropriate SS methods and other related parameters (epoch length, task complexity and dataset size) when analyzing dynamics of behavioral tasks.

51

Supplementary Materials for Study I Dynamics of task-related electrophysiological networks: a benchmarking study

Materials and Methods

Simulation Data. To quantify the performance of each method, we simulated four dynamic networks following the approach applied in (Kabbara et al., 2019;O'Neill et al., 2017). As a final step, to test the resistance of each method to noisy data, we added a random Gaussian noise to 𝑆𝑠𝑖𝑚 (𝑡) to have seven different levels of noise defined by signal to noise (SNR) ratio varied from 0dB to 3dB with step of 0.2dB yielding to 16 different noise levels.

Optimal Number of Components (Nopt) selection.

As in empirical data, the same quantification methods were used to select the optimal number of components relative to each method (DIFFIT and Elbow criterion for Kmeans).

Quantitative criteria. Several performance criteria were developed in this study to quantitatively compare the SS methods:

 Similarity. Source separation methods were applied on 𝑆𝑠𝑖𝑚 (𝑡) to extract spatial and temporal components. Temporal similarity was measured between 1) simulated and reconstructed timecourses by calculating the mean squared error between them, and spatial similarity between 2) simulated and reconstructed adjacency matrices using 'siminet' algorithm (Mheich et al., 2018).

Then, to describe the overall performance of each method in a single run, we combined both spatial and temporal similarities and average them over the 4 simulated components to obtain a single parameter called 'global similarity'. Note that for each method, the four spatial masks (relative to the four simulated networks) were applied to all NCopt networks in order to assimilate one component to each simulated network based on the highest spatial similarity. Following this, spatial, temporal and global similarities were computed specifically between assimilated networks.

 Reproducibility. To quantify methods consistency, that is the extent to which each method is able to reproduce same results across several runs, each method is repeated 100 times. We obtain 100 values for global similarity. Then, the interquartile range value was estimated to measure reproducibility. The higher the value, the less consistency the method is.

 Computation Time. For each method, mean and standard deviation of the elapsed time over 100 runs were calculated. The higher the mean value, the slower the method is.

 Noise Dependency. Because brain data is inherently noisy, it is important to quantify the ability of methods to extract correct components in presence of different noise levels. For this reason, global similarity average is also calculated at SNR value.

Results

Simulation. Table S1.

Primary visual networks (theta and high alpha).

Connections between right frontal areas and temporal areas (theta band).

Motor network (beta band).

A bilateral temporal connectivity network (alpha band).

Right lateralized temporo-parietal network (alpha).

Language-related network.

Visual to parietal (alpha band).

Connections between left frontal regions and right temporal regions.

Source Separation Methods details

ICA subtypes:

All ICA algorithms start with preprocessing steps including centering, whitening, and dimensionality reduction. Then, each method performs the source separation following its own criteria and hypothesis as described below.

 JADE: Joint Approximate Diagonalization of Eigen-matrices JADE ICA [START_REF] Cardoso | An Efficient Technique For The Blind Separation Of Complex Sources[END_REF] applies the Jacobi technique to perform joint approximate diagonalization on the Fourth Order (FO) cumulant of data sample to separate the 'k' independent source signals. FO cumulants are used as a measure of non-Gaussianity motivated by the fact that Gaussian distributions possess zero excess kurtosis. Thus, JADE seeks to estimate sources that possess high values of excess kurtosis.

 InfoMax

InfoMax ICA(Bell and Sejnowski, n.d.) belongs to the minimization of mutual information family of ICA algorithms that aims to maximize output entropy. It works as a simple learning algorithm for a feedforward neural network. Let 𝑌 = {𝑦 𝑖 } be the output of the neural network. The joint entropy of Y is:

𝐻(𝑌) = ∑ 𝐻(𝑦 𝑖 ) 𝑖 -𝐼(𝑌) (1)
Where 𝐻(𝑦 𝑖 ) are the marginal entropies of the outputs and 𝐼(𝑌) is their mutual information.

Therefore, referring to this equation, we can see that maximization of the joint entropy implies minimization of outputs mutual information.

 SOBI: Second Order Blind Identification

Second Order Blind Identification (SOBI) is based on the second order (SO) cumulant. It assumes that the sample data is collected from a set of temporally uncorrelated sources and tries to extract them by joint diagonalization of an estimate of several correlation matrices.

 FastICA

FastICA [START_REF] Hyvarinen | Fast and robust fixed-point algorithms for independent component analysis[END_REF] is one of the most popular ICA methods that maximizes non-Gaussianity as a measure of statistical independence, in a fixed-point iterative scheme. FastICA enhances non-Gaussianity by 'Negentropy' maximization. This is defined as the difference in entropy between the input distribution and the Gaussian distribution with the same mean and variance:

𝐽(𝑋 𝑤 ) = 𝐻(𝑋 𝑔𝑎𝑢𝑠𝑠 ) -𝐻(𝑋 𝑤 ) (2)
Practically, 'negentropy' measure can be approximated by the mean of a non-quadratic function. In our study, we chose the 'hyperbolic-tangent' non-linearity function. We opt for the original 'deflation' mode of FastICA where only one component is extracted at a time with random initialization. In order to have more consistent and reliable results, we used the ICASSO technique (Himberg and Hyvarinen, 2003) that applies FastICA algorithm several times (100 in our case).

 CoM2: The Contrast Maximization

CoM2 [START_REF] Comon | Independent component analysis, A new concept? Signal Processing[END_REF] explicitly maximizes a contrast function (Eq. 3) built from the square modulus of FO cumulant iteratively by applying a planar Givens rotation to every whitened signal pair until convergence. This procedure involves rooting successive polynomials of degree 4.

𝐹 𝐶𝑜𝑀2 (𝑋 𝑤 ) = ∑ |𝐶𝑈𝑀 𝑖,𝑖,𝑖,𝑖,𝑋𝑤 | 2 𝑘 𝑖=1
(3) where 𝑋 𝑤 represents the whitened input matrix and 𝐶𝑈𝑀 𝑖,𝑖,𝑖,𝑖,𝑋𝑤 represents the FO cumulant of the i th source of 𝑋 𝑤 .

 PSAUD: Penalized Semi-Algebraic Unitary Deflation

PSAUD method is a new delaftion ICA algorithm proposed by Becker et al. [START_REF] Becker | A Penalized Semialgebraic Deflation ICA Algorithm for the Efficient Extraction of Interictal Epileptic Signals[END_REF]. In a similar context to CoM2, this method optimizes a contrast function based on FO cumulant with alternating Jacobi iterations. However, an extra penalization term is added in the contrast function to ensure that sources of interest are extracted first. This term is related to CCA approach and includes autocorrelation of relevant sources to be extracted.

PCA:

Mathematically, to estimate mixing matrix 'A' and temporal signals 'S', eigenvectors and eigenvalues are computed from the correlation matrix of 'X' defined as:

𝑅 = 𝐸[𝑋𝑋 𝑇 ] (4)
In this way,

𝑋 = 𝑈 × ∧ × 𝑉′ (5)
Where, 'Λ' is a diagonal matrix of the same dimension of 'X' with non-negative diagonal elements in decreasing order of variance, 'U' is a unitary matrix containing orthonormal eigenvectors and 'V' the corresponding temporal coefficients. Thus, the mixing matrix 'A' is computed as the first 'k' elements of 𝑈 × ∧, and 'S' as the first 'k' elements of 'V'' coefficients.

Kmeans:

Mathematically, by considering G-dFC spatial connectivity as a set of m-dimensional vectors {C1, C2, ... , Cm}, the Kmeans clustering aims to generate a partition into 'k' clusters A = {A1, A2, ... , Ak} to minimize the sum of within-cluster distances 'J' over multiple iterations J:

𝐽 = ∑ ∑ 𝑑(𝐶𝑗, 𝜇𝑖) 𝐶𝑗⋲𝐴𝑖 𝑘 𝑖=1 (6)
where μi is the mean of spatial connectivity in Ai, and d represents the distance between the two vectors.

Null Distribution -Sign-flipping Approach

In order to construct the null distribution, we selected randomly half of the subjects and flipped the sign (multiplied by -1) of their temporal signals obtained after SS decomposition. To give all possible realizations, each time a different set of subjects was selected for sign flipping, yielding to 𝐶 𝑘 𝑛 = 𝑛! (𝑘! (𝑛 -𝑘)!) ⁄ possible combinations (n refers to the total number of subjects and k refers to half number of subjects) (𝐶 7 15 =6435 combinations for the self-paced data and 𝐶 9 19 =92378 for the Sternberg data). Regarding HCP Motor data, the same combination procedure could not be conducted directly on the total number of subjects for computational reasons. Therefore, we split the 61 subjects

Selection of optimal Number of Components (NCopt)

Elbow criterion

The elbow criterion is widely used to derive the number of clusters. Kmeans algorithm is repeated for different values of number of clusters 'k'. At each time, the ratio between within-cluster distance (dispersion in the cluster) to between-cluster distance (total variability outside that cluster) is calculated. Then, elbow criterion determine the input 'k' that maximizes within-cluster similarity and between-cluster dissimilarity.

DIFFIT Method

DIFFIT method is based on goodness of fit approach used to determine the optimal number of components. DIFFIT (referring to the difference in data fitting) is calculated based on model reconstruction error and the explained variance. Let us denote NC=J as the input number of component to the SS method. The reconstructed error of SS model at NC=J is defined as:

𝑅𝑒𝐸𝑟𝑟(𝐽) = |𝑋-𝑋𝑟𝑒𝑐(𝐽)| 𝐹 |𝑋| 𝐹 (7)
Where X is the original concatenated dFC matrix to be decomposed, 𝑋𝑟𝑒𝑐 (𝐽) = 𝐴𝑥𝑆 is the reconstructed matrix after SS method application at NC=J and || 𝐹 is the Frobenius norm.

Then, the DIFFIT is calculated as:

𝐹𝑖𝑡(𝐽) = 1 -𝑅𝑒𝐸𝑟𝑟(𝐽) (8)
𝐷𝐼𝐹(𝐽) = 𝐹𝑖𝑡(𝐽) -𝐹𝑖𝑡(𝐽 -1) (9)

𝐷𝐼𝐹𝐹𝐼𝑇(𝐽) = 𝐷𝐼𝐹(𝐽) 𝐷𝐼𝐹(𝐽+1) (10)
Then, the model J with the largest DIFFIT value is considered as the optimal NC.

Supplementary References:

Becker, H., Albera, L., Comon, P., [START_REF] Becker | A Penalized Semialgebraic Deflation ICA Algorithm for the Efficient Extraction of Interictal Epileptic Signals[END_REF]

Introduction

The human brain is currently recognized as a complex network in which spatially separated brain regions are functionally and dynamically communicating during tasks (Bola and Sabel, 2015;Hassan et al., 2015;O'Neill et al., 2017) and at rest (Allen et al., 2014;Baker et al., 2014;Kabbara et al., 2017).

The characterization of transient (dynamic) networks is of utmost importance to better understand the brain processes in healthy brains, as well as in neurological disorders (Kim et al., 2017;[START_REF] Sitnikova | Short timescale abnormalities in the states of spontaneous synchrony in the functional neural networks in Alzheimer's disease[END_REF]. Hence, the analysis of large-scale dynamic functional connectivity (dFC) has become a key goal in cognitive and clinical neuroscience. Among existing neuroimaging modalities, electro/magnetoencephalography (EEG/MEG) has major benefits in exploiting the spatiotemporal organization of the brain, since it provides a direct and non-invasive measure of electrical activity at the -sub-millisecond timescale [START_REF] Brookes | Altered temporal stability in dynamic neural networks underlies connectivity changes in neurodevelopment[END_REF]Vidaurre et al., 2018).

The This is partially due to the fact that a quantitative evaluation of the various methods' robustness is still missing, as 'ground truth' is almost impossible to get when studying empirical data. Therefore, validating the above-mentioned 'three-steps' pipeline using simulated EEG data is needed. In this context, some studies used a toy model where brain sources activity is modeled by oscillatory sinusoidal and Gaussian functions [START_REF] Halder | Quantitative Evaluation in Estimating Sources Underlying Brain Oscillations Using Current Source Density Methods and Beamformer Approaches[END_REF] or multivariate autoregressive (MVAR) filters with volume conductor head models generating pseudo-EEG data [START_REF] Anzolin | Quantifying the Effect of Demixing Approaches on Directed Connectivity Estimated Between Reconstructed EEG Sources[END_REF][START_REF] Haufe | A Simulation Framework for Benchmarking EEG-Based Brain Connectivity Estimation Methodologies[END_REF]. However, such models are linear and too simplistic compared to the complexity of real brain activity. On the other hand, few studies executed a preliminary performance quantification of dimensionality reduction using a set of pre-defined connectivity matrices, considered as simulated brain networks without any constraint on how these matrices are computed (Kabbara et al., 2019;O'Neill et al., 2017;Tabbal et al., 2021;Zhu et al., 2020). Here, we address these challenges by simulating the whole brain EEG dynamics by building physiologically inspired networks based on a realistic neural mass model using a human brain computational model, named 'COALIA' [START_REF] Bensaid | COALIA: A Computational Model of Human EEG for Consciousness Research[END_REF]. It contains several subtypes of GABAergic neurons (VIP-SST-and PV-neurons) and Glutamatergic (Pyramidal cells) neurons coupled through Diffusion Tensor Imaging (DTI) based structural connectivity, to simulate realistic HD-EEG data Technically, EEG simulated signals were computed via the forward model as previously shown by recent works [START_REF] Allouch | Mean-field modeling of brain-scale dynamics for the evaluation of EEG source-space networks[END_REF][START_REF] Bensaid | COALIA: A Computational Model of Human EEG for Consciousness Research[END_REF]. In this paper, we aim to investigate to what extent we can effectively track time-varying brain networks that dominate neuronal activity using HD-EEG scalp signals, specifically on a short timescale that commensurate with cognitive tasks. To this end, we benefit from the COALIA model to simulate a set of dynamic brain network states (BNS) related to a fast-scale cognitive task (picture naming). For this purpose, coherent gamma oscillations were produced in cortical areas of the ventral visual pathway following a picture naming dynamic scenario, as observed in real human data (Hassan et al., 2015).

The Local Field Potentials (LFP) of 20 subjects were simulated and the corresponding HD-EEG signals were calculated. In this context, we evaluated the performance of some inverse/connectivity algorithms (wMNE/eLORETA, PLV/wPLI/AEC) in addition to the quantification of many source separation methods (3 variants of ICA, PCA, NMF, and Kmeans). To our knowledge, no previous study has validated the performance of this 'three-steps' pipeline with the existence of ground truth provided by realistic HD-EEG simulation. We hope that this paper can be used as a benchmark for researchers who intend to investigate methodological or experimental issues to optimize EEG/MEG dynamic networks analysis.

Materials and Methods

The full pipeline of our study is divided into six steps as summarized in Figure 1.

The 'upgraded' COALIA Model

Simulations

The simulated cortical-level activity was generated using an updated version of the computational (iii) In order to improve the realism of simulated functional connectivity, COALIA was upgraded using the averaged structural connectivity matrix obtained in the HCP project (Human Connectome Project https://www.humanconnectome.org/) onto 487 adult subjects by DTI (Van Essen et al., 2013). We used the averaged fractional anisotropy measures to connect all neural masses. We used an average velocity of 7.5 m/s and the distance matrix to calculate the delay matrix.

The reader can refer to [START_REF] Bensaid | COALIA: A Computational Model of Human EEG for Consciousness Research[END_REF] for a detailed description of the original COALIA model.

Scenario

As we are interested in tracking task-related dynamic brain networks states (BNS), we consider the picture-naming task dynamic scenario, inspired by a previous publication (Hassan et al., 2015). Thus, the input of the COALIA model is a set of six simulated networks consecutive in time from 0 to 535ms.

At each time interval Ti (i=[1;6]), different regions of interest (ROIs) are activated as detailed in Figure 2. To be inactivated or activated, the parameters of each neural mass are tuned accordingly to generate background or gamma activity, based on the Desikan-Killiany atlas (Desikan et al., 2006).

The activity of brain sources is simulated in the gamma band as it is shown to be the most relevant frequency band in the context of such cognitive task [START_REF] Fell | The role of phase synchronization in memory processes[END_REF]Rodriguez et al., 1999), with a sampling frequency of fs=1024Hz. To perform a group-level brain networks analysis, LFP signals are simulated for nSubs=20 subjects with nTrials=100 trials for each subject. To consider an intra-subjects variability on the EEGs (between the nTrials), the input noise parameters (mean and variance) for all neural masses are randomly adjusted (uniform law) by ± 20% on each simulated LFP.

To consider an inter-subject variability (between the nSubs), each value of the connectivity matrix, used to link the 66 neural masses, has been randomly (uniform law) increase or decrease by 10%. Each LFP epoch lasts for 2 seconds; including 1-second pre-stimulus and 1-second post-stimulus.

We checked that every neural mass involved in the scenario (Figure 2

Forward Model

In order to generate simulated EEG data ('𝑋 𝑡 ') from simulated cortical time series ('𝑆 𝑡 '), we solved the forward problem as follow:

𝑋 𝑡 = 𝐺 × 𝑆 𝑡 (1)
where '𝐺' is the lead field matrix to be computed, which describes the electrical and geometrical characteristics of the head. To this end, T1 magnetic resonance imaging MRI (Colin27 template brain, [START_REF] Holmes | Enhancement of MR images using registration for signal averaging[END_REF]) was realigned and segmented. Then, realistically shaped shells representing brain, skull, and scalp surfaces were prepared to build a realistic head model using the Boundary Element Method (BEM) provided by the OpenMEEG package (Fieldtrip, (Gramfort et al., 2010)). In this work, we used 257 electrodes density (EGI, Electrical Geodesic Inc.) for EEG electrodes configuration.

Therefore, the computed lead field matrix '𝐺' (dimension:66 × 257) describes the physical current

𝑊 = 𝐵𝐺 𝑇 (𝐺𝐵𝐺 𝑇 + 𝛼𝐶) -1 (4)
Where '𝐵' is the diagonal matrix related to source covariance (inversely proportional to the norm of lead field vectors), '𝐺' is the lead field matrix, '𝛼' is the regularization parameter (inversely proportional to signal to noise ratio and set to 0.3 in our case) and '𝐶' is the noise covariance matrix (computed from one-second pre-stimulus baseline). In this work, we used the Matlab function implemented in Brainstorm toolbox codes (https://neuroimage.usc.edu/brainstorm) to compute wMNE (Tadel et al., 2011), with the signal to noise ratio set to 3 and depth weighting value to 0.5 (default values).

Exact low-resolution brain electromagnetic tomography (eLORETA)

The exact low-resolution electromagnetic tomography (eLORETA) [START_REF] Pascual-Marqui | Discrete, 3D distributed, linear imaging methods of electric neuronal activity[END_REF]) is a genuine inverse solution that gives more importance to the deeper sources with reduced localization error. In this way, it does not only consider depth bias but also provides exact localization with zero error even in the presence of measurement and structured biological noise.

𝑊 = (𝐺 𝑇 (𝐺𝐵𝐺 𝑇 + 𝛼𝐶) -1 𝐺) 1 2 (5)
Here, eLORETA was evaluated using FieldTrip toolbox (http://www.fieldtriptoolbox.org) following the approach of (Pascual-Marqui, 2007) with the default regularization parameter set to 0.05.

Source Space Resolution

As previously stated, in this study, we used a source-space including only the centroids of 66 desikan regions. However, another approach adopted consists of solving forward/inverse problems within a cortical mesh of high resolution (15000 vertices in most cases), followed by a projection onto an anatomical framework where all lead field vectors belonging to a common atlas region are added (Hassan et al., 2014). In this case, the inverse problem is ill-posed (15000 cortical sources>>257 EEG sensors). Thus, in our simulation study, we addressed this technical point to assess the impact of source-space resolution used.

Dynamic Functional Connectivity (dFC)

We estimated the functional connectivity between reconstructed regional time series. As we aim to undertake the dynamics of brain states, we used a sliding window approach (𝑤𝑖𝑛𝑑𝑜𝑤. 𝑠𝑖𝑧𝑒 = 𝛿,𝑤𝑖𝑛𝑑𝑜𝑤. 𝑠𝑡𝑒𝑝 = ∆), where connectivity is measured within each temporal window. In the case of phase synchronization and based on (Lachaux et al., 2000), the smallest number of cycles recommended to have a compromise between good temporal resolution and good accuracy is 6. Thus, as we are working in the gamma band (𝑐𝑒𝑛𝑡𝑟𝑎𝑙 𝑓𝑟𝑒𝑞𝑢𝑒𝑛𝑐𝑦 = 35𝐻𝑧), 𝛿 is equal to 0.17sec. ∆ is set to 0.017sec considering 90% overlapping between consecutive windows. Therefore, the total number of windows over the whole epoch duration is nWinds=108 windows.

In this paper, we evaluated two popular modes of connectivity including phase-based metrics: (1) PLV (phase-locking value) and ( 2) PLI (phase-lag index), and one amplitude-based metric: (3) AEC (amplitude envelope correlation), as described below.

Phase-locking value (PLV)

For each trial, the phase-locking value (Lachaux et al., 2000) that characterizes the phase relationship between two signals 𝑥(𝑡) and 𝑦(𝑡) is expressed as follows:

𝑃𝐿𝑉(𝑡) = | 1 𝛿 ∫ 𝑒 𝑗(𝜑 𝑦 (𝑡)-𝜑 𝑥 (𝑡)) 𝑡+ 𝛿 2 𝑡- 𝛿 2 𝑑𝜏| (6)
Where 𝜑 𝑦 (𝑡) and 𝜑 𝑥 (𝑡) represent the instantaneous phases of signals y and x respectively derived from the Hilbert transform at time window t. The Matlab function used to calculates PLV following Equation ( 6) is available on GitHub (https://github.com/judytabbal/dynCOALIA).

Weighted phase-lag index (wPLI)

Unlike PLV, the weighted phase-lag index (wPLI) is insensitive to zero-lag interaction (Vinck et al., 2011). This index is based only on the imaginary component of the cross-spectrum and is thus robust to noise [START_REF] Peraza | Volume conduction effects in brain network inference from electroencephalographic recordings using phase lag index[END_REF]. Here, we used the Fieldtrip toolbox to compute wPLI (multi-taper method, fast Fourier transform, single Hanning tapper, 2Hz frequency resolution). As a certain amount of averaging across trials is required, wPLI was calculated at each temporal window using all trials of interest for each subject.

Amplitude Envelope Correlation (AEC)

The amplitude envelope correlation (AEC) is calculated using the Hilbert transform of regional time series. Pearson correlation is then computed between the amplitude envelopes of two pairs of regions As a result, the output dimension of the dFC tensor was [𝑛𝑅𝑂𝐼𝑠 × 𝑛𝑅𝑂𝐼𝑠 × 𝑛𝑊𝑖𝑛𝑑𝑠] for every single trial in the case of PLV and AEC, and each subject in the case of wPLI. This tensor was unfolded into a 2D matrix [𝑛𝑅𝑂𝐼𝑠 × (𝑛𝑅𝑂𝐼𝑠 -1)/2 × 𝑛𝑊𝑖𝑛𝑑𝑠] by removing redundant connections due to symmetry, followed by mean row subtraction. Finally, all trials/subjects dFC matrices were concatenated along temporal dimension and a group dFC matrix is constructed denoted '𝑃'.

Dimensionality Reduction

This step is crucial to extract task-related brain network states (BNSs). It consists of summarizing all time-varying connectivity features in the constructed matrix '𝑃' into k dominant brain patterns over given time intervals. This problem can be formulated as follows:

𝑃 = 𝐴 × 𝑇 (7)
where '𝐴' is the mixing matrix illustrating the k spatial maps of dominant brain networks and '𝑇'

represents the corresponding temporal source signatures.

Among existing dimensionality reduction algorithms, we chose to investigate: (1) temporal Independent Component Analysis (tICA), ( 2) Principal Component Analysis (PCA), (3) Non-negative Matrix Factorization (NMF), and (4) Kmeans. They mainly differ in the constraints imposed on decomposed components. To reduce the effect of the number of states per method, we imposed k=6 components equal to the number of simulated networks in this work for all algorithms. We discuss this issue in the discussion section.

temporal Independent Component Analysis (tICA)

tICA approach is already used by several studies (O'Neill et al., 2017;Yaesoubi et al., 2015) to derive k brain states that are 'statistically mutually independent' in time. Here, we evaluated tICA using three prominent ICA sub-methods: (1) JADE (Rutledge and Jouan-Rimbaud Bouveresse, 2013), ( 2) FastICA [START_REF] Langlois | An Introduction to Independent Component Analysis: InfoMax and FastICA algorithms[END_REF] and ( 3) PSAUD [START_REF] Becker | A Penalized Semialgebraic Deflation ICA Algorithm for the Efficient Extraction of Interictal Epileptic Signals[END_REF]. Briefly, FastICA is based on information theory, while JADE and PSAUD tend to optimize contrast functions based on high statistical order cumulants of the data. Technically, we adopted the following functions; 'jader' for JADE [START_REF] Cardoso | High-Order Contrasts for Independent Component Analysis[END_REF], 'icasso' for FastICA (Himberg and Hyvarinen, 2003), and 'P_SAUD' for PSAUD [START_REF] Becker | A Penalized Semialgebraic Deflation ICA Algorithm for the Efficient Extraction of Interictal Epileptic Signals[END_REF], implemented in Matlab (The Mathworks, USA, version 2019a).

Principal Component Analysis (PCA) Inverse Model/Functional connectivity combination evaluation

The 'spatial similarity' metric was used to quantify the performance of different inverse model/functional connectivity combination methods at the network level. Briefly, this measure takes into consideration the distribution of weights across edges between and within brain lobes (occipital, parietal, temporal, central, frontal, prefrontal, limbic). The reader can refer to supplementary materials (Figure S1) for a more detailed description of spatial similarity calculation. The proportional threshold was also used to retain the top-weighted edges, varied from top 1% to 2% of total undirected possible edges, which correspond to the possible number of connections between the selected range of nodes density (see section 6.1).

Dimensionality reduction methods evaluation

For each dimensionality reduction method, the 'temporal similarity' (correlation between signals) and 'spatial similarity' were computed to evaluate both the temporal and spatial performance of each method. Then, 'global similarity' was calculated as the average between both similarities.

Statistical Analysis

We also applied our pipeline with performance analysis at subject level data yielding a distribution set of values across 20 subjects. Then, the ANOVA test was used (The Mathworks, USA, version 2019a) to statistically quantify differences between the tested methods, followed by a post-hoc correction for multiple comparisons (using the built-in function 'multcomp') with Bonferroni correction (statistical significance p-value=0.01).

Inverse solutions and connectivity measures combination

To visualize results at the network level, we illustrate in Figure 4 the estimated dynamic functional connectivity matrices at time intervals for each inverse model/functional connectivity method. For visualization purposes, results were averaged over trials and subjects. For results quantification, we calculated the spatial similarity metric (see supplementary Figure S1) at each time interval between each method and corresponding reference networks. In our case, reference networks correspond to the functional connectivity method applied on simulated LFP sources data. For example, brain networks estimated from PLV applied on LFP sources were considered as a reference for wMNE/PLV and eLORETA/PLV. We followed the same concept for other methods. Then, spatial similarity values were averaged over time intervals to obtain a single quantification value for each method. At the group-level, global spatial similarities were 0.65 for wMNE/PLV, 0.62 for wMNE/wPLI, 0.61 for wMNE/AEC, 0.53 for eLORETA/PLV, 0.52 for eLORETA/wPLI and 0.50 for eLORETA/AEC (see Figure 6.A).

At the group level, one can notice that (1) the highest similarity was reached for wMNE/PLV while eLORETA/AEC performed the worst. ( 2) wMNE performed better than eLORETA for all functional connectivity measures and (3) PLV exhibited the highest similarity values followed by wPLI then AEC for both source reconstruction methods. Estimated networks were also averaged over trials for each subject and similarity was computed at subject-level data. Similarity results of the distribution over 20 subjects are displayed on boxplots in Figure 6.C Statistical analysis performed by ANOVA test showed non-significant differences between (1) wMNE/wPLI and wMNE/AEC (p>0.05) and ( 2) eLORETA/PLV and eLORETA/wPLI (p>0.05) while all other combinations showed significant differences (p<0.05). To globally quantify each method, we averaged maximal global similarity values obtained across all time intervals. At the group-level, maximal global similarity values were 0.73 for JADE, 0.74 for FastICA, 0.78 for PSAUD, 0.62 for PCA, 0.68 for NMF and 0.61 for Kmeans (see Figure 6.B).

Dimensionality reduction methods were also tested at subject-level data where all trial dFC matrices were concatenated along temporal dimensions for each subject. Maximal Global Similarity results distribution over 20 subjects are displayed on boxplots in Figure 6.D (mean values: 0.69 for JADE, 0.70 for FastICA, 0.73 for PSAUD, 0.63 for PCA, 0.68 for NMF and 0.63 for Kmeans). Statistical analysis performed by ANOVA test showed non-significant differences between (1) JADE and FastICA (p>0.05),

(2) JADE and NMF, and (3) PCA and Kmeans (p>0.05) while all other combinations showed significant differences (p<0.05).

Discussion

Tracking dynamic brain networks from non-invasive electrophysiological (EEG/MEG) data is a key challenge in neuroscience. A critical issue is to evaluate to what extent extracted 'brain network states' match those that are truly activated during tasks. In this context, the presence of 'ground truth' data is of utmost importance to ensure objective and quantitative evaluation of pipeline processes. The main contribution of this study is the use of simulated EEG data -in a dynamic context -from a physiologically grounded computational model upgraded with a new HBP structural connectivity Matrix as ground truth for evaluating methods performance in 're-estimating' correctly reference states. The dynamic approach was simulated by adopting the scenario of a picture-naming task evolving several brain states fluctuating over time (Hassan et al., 2015;[START_REF] Mheich | HD-EEG for tracking sub-second brain dynamics during cognitive tasks[END_REF]Mheich et al., , 2015)).

Methods evaluation 93

In this paper, we conducted a systematic evaluation analysis of two source reconstruction methods (wMNE/eLORETA), three functional connectivity measures (PLV/wPLI/AEC), and six dimensionality reduction methods (JADE/FastICA/PSAUD/PCA/NMF/Kmeans) used as a three-step pathway to extracting dynamic brain network states from scalp-EEG signals. Overall, results show that wMNE outperforms eLORETA in the context of our simulation study. This is observed when evaluating source reconstruction methods on both cortical-level (higher precision results achieved with wMNE) and network-level (higher spatial similarity values obtained with wMNE independently of the functional connectivity measure used). In this study, the wMNE/PLV combination exhibits the best performance leading to the highest spatial similarity value between reconstructed and reference networks among all inverse model/functional connectivity combinations tested. Interestingly, these results are in line with previous comparative studies showing consistency and robust results for wMNE/PLV combination in the context of EEG source space connectivity using real data from picture naming task (Hassan et al., 2014) and simulated data from epileptogenic-modeled networks [START_REF] Allouch | Mean-field modeling of brain-scale dynamics for the evaluation of EEG source-space networks[END_REF]Hassan et al., 2017). The strength of PLV results may reflect the potential mechanisms of zero-lag synchronization of neural activity already discussed in previous works [START_REF] Gollo | Mechanisms of Zero-Lag Synchronization in Cortical Motifs[END_REF][START_REF] Roelfsema | Visuomotor integration is associated with zero time-lag synchronization among cortical areas[END_REF]. It is worth noting that, in the particular context of the study, functional connectivity methods had less impact on results than inverse models. For instance, wPLI-AEC performed similarly with wMNE, and PLV-wPLI performed similarly with eLORETA.

We also explored an additional technical point related to the source space resolution used for lead field computation. In this context, there is no clear evidence about the ideal source grid resolution to adopt. Our results suggest that computing lead field after cortex parcellation (referred to as LowRes) provides higher accuracy results than lead field computation before parcellation followed by ROI projection (referred to as HighRes). One explanation might be related to the induced blurring effect between close sources when using a high-resolution source grid, which might increase the cross-talk between specific cortical locations [START_REF] Maldjian | Graph Theoretical Analysis of Resting-State MEG data: Identifying Interhemispheric Connectivity and the Default Mode[END_REF]Schoffelen and Gross, 2009).

Overall, performance analysis showed promising results for most dimensionality reduction methods regarding their ability to derive dominant brain states. In this work, global similarity results highlight the performance of ICA subtypes methods (PSAUD, JADE, and FastICA) relative to other methods.

These findings are in line with a recent comparative study between source separation methods applied on three independent real MEG datasets (Tabbal et al., 2021). Particularly, in this work, it is noteworthy to point that PCA, NMF, and Kmeans exhibit fragility in the temporal domain reflected by relatively low temporal similarities at several intervals (Figure S2). This is also reflected by the statistical comparative analysis between the evaluated methods at the level of spatial and temporal domain (Figure S4). For instance, although NMF and Kmeans were able to extract the correct spatial networks with a high precision relative to other techniques, these methods lack the potential to accurately track their fast temporal activity, which notably influenced the global methods performance (Figure 6). In contrast, ICA techniques were significantly more sensitive than other techniques to the temporal fluctuations of the simulated networks (high temporal similarity). We can notice an attractive behavior for PSAUD method in terms of spatial and temporal precision. Therefore, when studying the dynamics of the taskrelated brain activity, it is important to opt for a decomposition technique that assures a good spatial accuracy in the extracted functional brain networks on one hand, and follows delicately their dynamic aspect on the other hand.

Furthermore, it is crucial to note that most methods revealed weaker performance at specific time intervals (T2, T3, and T6) (Figure S3). Interestingly, these intervals are the narrowest (T2 lasts for 30ms, T3 for 40ms, and T6 for 55ms) compared to other simulated intervals (T1 lasts for 120ms, T4 for 130ms, and T5 for 160ms). This point is of particular interest as it shows method limitation when applied on a very fast timescale (~<100ms). The noise level was varied to test the effect of scalp-level noise on extracted brain states. Although most dimensionality reduction methods were robust against noise variation, a generalized stability assessment of each method relative to measurement noise needs more detailed investigation related to the type and level of the realistic added noise.

Finally, the subject-level and group-level results seem to be convergent (Figure 6.C, 6.D). We believe that these results are of significant neuroscientific interest, in particular, for application in clinical neuroscience, that needs consistent and reliable results at a patient-level. Therefore, it can open new avenues for detailed methodology/parameters evaluation and optimization on this framework at the single-subject data.

Methodological Considerations

In this study, we aim to carry out a dynamic analysis of brain network activity using a computational model serving as a ground truth for our pipeline methodology evaluation. To this end, one could define any set of brain network states as input to the COALIA model. Here, simulated states were determined based on a realistic picture naming task scenario, for which a solid background is available concerning activated brain regions and networks. Furthermore, the variety of brain regions activated distinctly from the onset to the reaction time is a major benefit of this scenario choice as we are interested in evaluating the dynamic behavior of methods at a realistic and rapid timescale. However, other scenarios could be also tested simulating other simple tasks (i.e. motor task), or complex tasks (i.e.

working memory task), or even customized scenarios including pre-defined brain states with manipulation in time length, number of activated regions, number of simulated states, etc… Besides, it should be noted that further investigations about the type/value of the noise added to the simulated gamma oscillations are needed to reduce, as much as possible, the gap between our brain model and the real brain physiology.

The main objective of this paper is to assess the ability of well-known existing methodologies to track dominant brain networks states from EEG signals, rather than to perform an exhaustive comparative analysis between all possible combinations of the 'three-step' pipeline methods. Nevertheless, our analysis work could be extended to cover a much more variety of techniques at each step. First, beamformer family methods such as linearly constrained minimum variance (LCMV) (Van Veen et al., 1997) could be evaluated besides minimum norm estimates methods tested in this paper. Second, as we accounted for the volume conduction effect in phase-based coupling measures by applying both PLV and wPLI methods, it would be also interesting to study the effect of leakage correction on amplitude-based coupling measures by testing the corrected version of amplitude envelope correlation (AEC-c). Moreover, we focused on 'functional' connectivity metrics that do not consider directionality property. In this context, the analysis of 'effective' connectivity that explores causality between brain regions may contribute to further insights. Furthermore, we adopted a sliding window approach to study the dynamic behavior of functional connectivity. Although window size was calculated corresponding to the minimal length required (Lachaux et al., 2000), the generalization of this finding to amplitude-based functional connectivity metrics may be critical, and therefore, further efforts are required to ensure a precise evaluation. Moreover, several findings criticize the performance of sliding windows when computed over very short epoch durations (Fraschini et al., 2016;Liuzzi et al., 2019). We, therefore, suggest for future studies to test high-resolution measures of functional connectivity that are sensitive to fast fluctuations (Tewarie et al., 2019) as an extension to our simulation study.

Third, to extract dominant brain states, we tested dimensionality reduction methods including source separation and clustering techniques. Nevertheless, it will be interesting to add other strategies that revealed accurate results in previous studies to such simulation analysis as Hidden Markov Model (HMM) (Vidaurre et al., 2018) and the multivariate autoregressive model (Casorso et al., 2019).

The connectivity matrices were thresholded using proportional thresholding to standardize our comparison across methods, trials, and intervals. Compared to absolute thresholding that depends on average connectivity values, proportional thresholding has been demonstrated to result in more stable network metrics (Garrison et al., 2015;[START_REF] Van Wijk | Comparing brain networks of different size and connectivity density using graph theory[END_REF]. In order to account for the effect of the threshold value, we performed our analysis on several threshold values as explained in the Materials and Methods section. One limitation is related to the fact that a priori knowledge about network density is not available in real experimental data contrary to simulated data. Thus, in the context of real data, the selection of the appropriate threshold values is critical and crucial for further investigation, which is out of the scope of this paper.

In addition, the number of reference states is known a priori in the simulation-based analysis, but critical and often hard to predict in the case of empirical data. The optimal number of brain network states could be defined in numerous ways based on several optimization criteria, such as crossvalidation criterion (Pascual-Marqui et al., 1995), Krzanowski-Lai criterion (Krzanowski and Lai, 1988), kneedle algorithm (Satopaa et al., 2011), and difference of data fitting (DIFFIT) (Timmerman and Kiers, 2000;Wang et al., 2018). However, to keep our analysis compact, we decided to fix this parameter to the exact number of simulated states (6 states) for all decomposition and clustering methods. Although this could pose a limitation concerning the optimum performance of each specific algorithm, the focus of the present work is the direct comparative evaluation of different algorithms based on pre-defined reference networks. For instance, applying a specific optimization criterion to simultaneously all algorithms may fit better some algorithms than others, and hence, influence the interpretation of our results. Therefore, future studies may test different existing approaches used to optimize the number of derived components relative to each algorithm in a unified framework that helps researchers determine the most accurate procedure to follow in such a context as ours.

Finally, in this work, we set the number of EEG channels (257 channels) since it established accurate source localization results relative to lower sensor densities [START_REF] Allouch | Mean-field modeling of brain-scale dynamics for the evaluation of EEG source-space networks[END_REF][START_REF] Song | EEG source localization: sensor density and head surface coverage[END_REF]. It can be however interesting to follow the work of such studies and examine the effect of sensor spatial resolution by decreasing successively the number of electrodes from high-density (257 channels) to low-density (19 channels) EEG signals.

Conclusion

In this study, we propose a complete framework for analyzing the task-related dynamic electrophysiological brain networks. We used a physiologically inspired full-brain model (named COALIA) as 'ground-truth' to systematically validate and optimize the pipeline study including EEGsource space connectivity estimation and dynamic brain network states extraction. As a proof of concept, the study was conducted in a dynamic scenario that emulates the picture naming task. Our findings suggest a good performance of wMNE/PLV combination to elucidate the appropriate functional networks. Results also revealed the promising efficiency of ICA techniques to derive relevant dynamic brain network states. We suggest using these realistic models to go further in the evaluation of the different steps and parameters involved in the EEG/MEG source-space network analysis. This can, to some extent, reduce empirical selection of inverse model, connectivity measure and dimensionality reduction method as some of the methods can have considerable impact on the final results and interpretation. 109

Introduction

Parkinson's disease (PD) can be associated with a broad spectrum of symptoms. Although it is mostly known for its effects on motor function, such as bradykinesia, rigidity, and resting tremor [START_REF] Hayes | Updates and advances in the treatment of Parkinson disease[END_REF], cognitive impairments associated with PD also have a huge impact on the quality of life [START_REF] Lawson | Cognitive decline and quality of life in incident Parkinson's disease: the role of attention[END_REF]. These impairments affect several cognitive functions, such as attention, memory, executive functions, visuospatial abilities, emotional processing, and language [START_REF] Halliday | The progression of pathology in Parkinson's disease[END_REF][START_REF] Papagno | Cognitive and behavioral disorders in Parkinson's disease: an update. I: cognitive impairments[END_REF][START_REF] Reid | Dementia in Parkinson's disease: a 20-year neuropsychological study (Sydney Multicentre Study)[END_REF]. Although cognitive functions are mostly evaluated through neuropsychological assessment, they are also studied through experimental tasks that can give a more precise insight into specific processes. For instance, alterations in inhibitory control have been described in PD by analyzing performances at the stop-signal reaction task [START_REF] Di Caprio | Early-Stage Parkinson's Patients Show Selective Impairment in Reactive But Not Proactive Inhibition[END_REF][START_REF] Gauggel | Inhibition of ongoing responses in patients with Parkinson's disease[END_REF][START_REF] Obeso | Deficits in inhibitory control and conflict resolution on cognitive and motor tasks in Parkinson's disease[END_REF]. Another well-studied function in PD is conflict resolution, a sub-process of cognitive control that allows the suppression of automatic responses in favor of voluntary actions (as considered in light of the dual-route and activation-suppression models, see [START_REF] Hommel | Towards a Unitary Approach to Human Action Control[END_REF]. PD patients have consistently shown alterations in this function, although the nature of the impairment is somewhat inconsistent between studies [START_REF] Cagigas | Flanker compatibility effects in patients with Parkinson's disease: Impact of target onset delay and trial-by-trial stimulus variation[END_REF][START_REF] Duprez | Impulsive oculomotor action selection in Parkinson's disease[END_REF][START_REF] Falkenstein | Effects of stimulus-response compatibility in Parkinson's disease: a psychophysiological analysis[END_REF][START_REF] Wylie | The effect of Parkinson's disease on the dynamics of on-line and proactive cognitive control during action selection[END_REF][START_REF] Wylie | Activation of conflicting responses in Parkinson's disease: evidence for degrading and facilitating effects on response time[END_REF].

PD disturbs cortical-subcortical loops involved in motor, cognitive, and limbic processes, and thus, PDrelated alterations in cognitive functioning are arguably associated with changes at both the cortical and subcortical levels. Several studies investigated the correlations between brain activity and behavioral changes, using a variety of neuroimaging modalities (fMRI, PET, M/EEG). For instance, it is now undisputed that conflict resolution is associated with activity in the dorso-lateral prefrontal cortex (DLPFC), the pre-supplementary motor area (pre-SMA), the anterior cingulate cortex, and the subthalamic nucleus. Although very informative, this localizationist approach suffers the bias of disregarding how brain regions interact with each other, thus precluding from a better understanding of how neurodegenerative diseases such as PD alter cognitive functions. Indeed, a lot of evidence now shows that cognitive functioning emerges from the communication of distant brain regions [START_REF] Bassett | Network neuroscience[END_REF]. Focusing on these brain networks to better understand cognitive (dys)functioning is crucial given their association with neurological disorders (Fornito et al., 2015). Brain networks can be studied via the estimation of functional connectivity (FC) between brain areas estimated on electrophysiological (M/EEG) or metabolic (MRI, PET) signals. FC is not in itself a direct measure of communication between brain areas, but rather reflects statistical dependencies of brain activity between different regions. M/EEG FC is particularly interesting because it is usually inferred through the phase synchronization of neural oscillations, a mechanism that has been proposed to facilitate communication between neuronal assemblies (Fries, 2015).

So far, most studies have reported alterations in the static FC associated with PD using resting-state fMRI [START_REF] Baggio | Cognitive impairment and resting-state network connectivity in Parkinson's disease[END_REF][START_REF] Baggio | Functional brain networks and cognitive deficits in Parkinson's disease[END_REF][START_REF] Lopes | Cognitive phenotypes in parkinson's disease differ in terms of brainnetwork organization and connectivity[END_REF][START_REF] Skidmore | Connectivity brain networks based on wavelet correlation analysis in Parkinson fMRI data[END_REF][START_REF] Wolters | Resting-state fMRI in Parkinson's disease patients with cognitive impairment: A metaanalysis[END_REF], or resting-state M/EEG [START_REF] Bertrand | Brain Connectivity Alterations Are Associated with the Development of Dementia in Parkinson's Disease[END_REF][START_REF] Bosboom | MEG resting state functional connectivity in Parkinson's disease related dementia[END_REF]Hassan et al., 2017). However, the literature focusing on dynamic FC evaluation using M/EEG neuroimaging techniques during cognitive tasks in PD is scarcer, although it would surely provide valuable insights on the neurophysiological alterations caused by the disease and its relationship with behavioral changes. Even if M/EEG doesn't benefit from the same spatial resolution as fMRI, recent advances in cortical source reconstruction allow for enough details in inferring cortical area long-range FC (Hassan and Wendling, 2018). One strong advantage of M/EEG techniques is their important temporal resolution. This aspect is fundamental when studying cognitive processes that are inherently dynamic. Conflict resolution is a process that has dynamic properties that are evident at the behavioral level: action selection and suppression are time-resolved processes allowing conflict resolution, and these aspects were usually masked by focusing on task-averaged behavioral performances in princeps studies.

As cognitive processes, FC is also time-dependent and brain networks dynamically rearrange themselves in rest and task (Baker et al., 2014;Bola and Sabel, 2015;de Pasquale et al., 2016de Pasquale et al., , 2010;;Hassan et al., 2015;[START_REF] Kabbara | The dynamic modular fingerprints of the human brain at rest[END_REF]O'Neill et al., 2018O'Neill et al., , 2017)), with consequences on behavior (Allen et al., 2018). Several new methods now allow investigating how inter-regional communication varies with time [START_REF] Sizemore | Dynamic graph metrics: Tutorial, toolbox, and tale[END_REF]Tabbal et al., 2021). Such methods would tremendously help in understanding how dynamic cognitive processes such as conflict resolution unravel and how this would be affected by neurodegenerative diseases such as PD. In this study, we hypothesize that specific cognitive processes depend on brain networks that dynamically rearrange and that PD is associated with changes in these dynamic properties. We used the example of a classic cognitive task with the Simon task (which informs on conflict resolution) paired with high density EEG (HD-EEG, 256 channels) reconstructed at the cortical source level in a group of PD patients and healthy controls (HC).

We combined the calculation of dynamic FC matrices with a dimension reduction method (independent component analysis, ICA) and source-level microstate approach to investigate and side (left/right) of the correct button press associated with the color (yellow/blue). For example, the trial is considered congruent when a blue circle appears on the right side, and incongruent when a blue circle appears on the left side (Figure 1).

After a training session of familiarization with the task, participants executed 10 blocks of 60 trials, with a pause every three blocks to avoid fatigue and to check EEG electrodes' impedance. In total, 600 trials were performed with 300 congruent and 300 incongruent trials with a pseudo-randomized display. However, in this study, we were only interested in considering the incongruent trials that were responded correctly, as they are the ones associated with efficient control of the strongest conflict.

Data acquisition and Preprocessing

EEG signals were recorded using a high-density EEG (HD-EEG) system (EGI, Electrical Geodesic Inc., 256 channels), with a sampling frequency of 1000 Hz. We kept 199 electrodes as shown in the channel file presented in the GitHub repository (https://github.com/judytabbal/dynCogPD).

Since EEG signals are often contaminated by several sources of noise and artifacts, a preprocessing step is required before analysis. Here, we performed manually the EEG preprocessing using the Brainstorm toolbox (Tadel et al., 2011).

First, we applied a DC offset removal; followed by a notch filter (50Hz) and a band-pass filter (1-100Hz).

Then, EEG data was visually inspected to identify, remove and interpolate the bad channels.

Independent Component Analysis (ICA) was used to remove the remaining eye blinks or muscle artifacts. Here, we segmented the recorded signals into epochs relative to the stimulus appearance (onset), with a 700 ms baseline (pre-stimulus) and 1200 ms post-stimulus. Finally, after another visual inspection, epochs with excessive remaining noise were rejected.

Consequently, for each subject, 174 incongruent trials on average (STD=49.2) of 1900 ms length [-700; +1200ms] were conserved.

After EEG preprocessing, several steps were applied to construct the dynamic brain network states relative to each group (HC and PD) as summarized in Figure 2 and explained in the following sections. Then, the dynamic functional connectivity (dFC) was estimated for each subject and trial using the sliding window approach. Phase Locking Value measure (PLV) was used to compute the statistical coupling between ROIs. Finally, the temporal Independent Component Analysis (tICA) was applied on

In this work, we used the Matlab function implemented in the Brainstorm toolbox to compute wMNE (Tadel et al., 2011). The 𝑆𝑁𝑅 was set to 3 and the depth weighting value to 0.5 (default values).

Dynamic Functional Connectivity: PLV -sliding window

Next, several approaches have been proposed to compute the functional connectivity between reconstructed cortical regions. In this study, we used the phase-locking value (PLV) method (Lachaux et al., 2000). As we aim to assess the dynamics of functional connectivity, we adopted the sliding window approach defined by its length 𝛿 with an overlapping step ∆. Hence, for each trial, PLV measures the phase synchronization between two signals 𝑥(𝑡) and 𝑦(𝑡) within each temporal window t:

𝑃𝐿𝑉(𝑡) = | 1 𝛿 ∫ 𝑒 𝑗(𝜑 𝑦 (𝑡)-𝜑 𝑥 (𝑡)) 𝑡+ 𝛿 2 𝑡- 𝛿 2 𝑑𝜏| (4)
Where 𝜑 𝑦 (𝑡) and 𝜑 𝑥 (𝑡) are the instantaneous phases of signals 𝑦(𝑡) and 𝑥(𝑡) respectively derived from the Hilbert transform. The Matlab function used to calculate the following PLV equation is available on Github (https://github.com/judytabbal/dynCogPD).

The choice of temporal window length was based on (Lachaux et al., 2000), where it is recommended to have at least 6 cycles at the given frequency band as a compromise between temporal and spatial accuracy. In this study, we conducted our analysis separately in both beta band [12-25Hz] (central frequency: 𝐶 𝑓 = 18.5𝐻𝑧) and gamma band [30-45Hz] (central frequency: 𝐶 𝑓 = 37.5𝐻𝑧). Thus, we chose the smallest window length equals to 6 𝐶 𝑓 ⁄ , that is 320ms for the beta band and 160ms for the gamma band. We considered a 90% overlapping between consecutive windows, as a way to track the very fast neural activity. Therefore, the total number of windows over the whole trial duration was 𝑛𝑊𝑖𝑛𝑑𝑠 = 49 windows for beta band and 𝑛𝑊𝑖𝑛𝑑𝑠 = 109 windows for gamma band, and the output dimension of the dynamic functional connectivity (dFC) was [𝑛𝑅𝑂𝐼𝑠 × 𝑛𝑅𝑂𝐼𝑠 × 𝑛𝑊𝑖𝑛𝑑𝑠] for each subject trial.

Dynamic Brain Network States (dBNS)

Temporal Independent Component Analysis (tICA): JADE

For each subject trial, dFC tensor can be unfolded into a 2D matrix of dimension [𝑛𝑅𝑂𝐼𝑠(𝑛𝑅𝑂𝐼𝑠 -1)/ 2 × 𝑛𝑊𝑖𝑛𝑑𝑠] due to symmetry. Then, for each group, the resultant dFC matrices of all trials and subjects were concatenated along the temporal dimension to generate a group-specific dFC matrix denoted 𝑀. For instance, 𝑀 𝑐 refers to the dFC relative to the HC group, while 𝑀 𝑝 refers to the dFC relative to the PD group.

We aim to summarize and extract the most relevant time-varying connectivity patterns in both 𝑀 𝑐 and 𝑀 𝑝 . This problem can be formulated as follows:

𝑀 = 𝐴 × 𝐵 (5)
Where 𝐴 is the mixing matrix that represents the 'k' spatial maps of dominant brain network states and 𝐵 describes their temporal evolution. Among the existing decomposition and clustering techniques, we derived the dynamic Brain Network States (dBNS) using temporal Independent

Component Analysis (tICA) adopted by several previous studies (O 'Neill et al., 2017;Yaesoubi et al., 2015). This technique assumes maximal independence between the time courses of the extracted dBNS. Here, tICA was performed using the JADE algorithm (Joint Approximate Diagonalization of Eigenmatrices) [START_REF] Cardoso | An Efficient Technique For The Blind Separation Of Complex Sources[END_REF]Rutledge and Jouan-Rimbaud Bouveresse, 2013). Briefly, JADE applies the Jacobi technique to optimize contrast functions based on high statistical order (Fourth Order: FO) cumulants of the data. We used the implemented JADE function in Matlab (The Mathworks, USA, version 2019a).

Number of states selection

Determining the optimal number of states to be extracted by tICA is a crucial issue, usually for most decomposition and dimensionality reduction methods [START_REF] Cong | Validating rationale of group-level component analysis based on estimating number of sources in EEG through model order selection[END_REF][START_REF] Mørup | Automatic relevance determination for multi-way models[END_REF].

Here, we used the DIFFIT (difference in data fitting) method based on the goodness of fit approach (Timmerman and Kiers, 2000;Wang et al., 2018), previously used by recent studies (Tabbal et al., 2021;Tewarie et al., 2019;Zhu et al., 2020). DIFFIT is calculated based on the following equations:

𝐹𝑖𝑡(𝐽) = 1 - |𝑀 -𝑀′(𝐽)| 𝐹 |𝑀| 𝐹 (6)
𝐷𝐼𝐹𝐹𝐼𝑇(𝐽) = 𝐹𝑖𝑡(𝐽) -𝐹𝑖𝑡(𝐽 -1) 𝐹𝑖𝑡(𝐽 + 1) -𝐹𝑖𝑡(𝐽)

Where 𝑀 is the original concatenated dFC matrix to be decomposed and 𝑀′ is the reconstructed matrix after tICA decomposition at the number of states J, and || 𝐹 is the Frobenius norm. In this study, we varied J from 3 to 10 states, then chose the number of states J that gives the largest DIFFIT value. For each group, we applied the DIFFIT method on the dFC matrix of each subject separately, then averaged the obtained DIFFIT values across subjects to obtain an average of 5 states for both HC and PD groups (mean=5.3 for HC, mean=4.9 for PD).

Significant task-modulated states

A further step was applied to automatically select, among all extracted ICA states, those that are significantly modulated by the task for each group. First, we followed the procedure adopted by (Hunt et al., 2012;O'Neill et al., 2017;Tabbal et al., 2021;Winkler et al., 2014) to build an empirical null distribution through the generation of a surrogate time course based on a sign-flipping permutation.

An ICA state was deemed significant if its corresponding time course fell outside the distribution for a duration determined by 3 successive cycles (a cycle is determined based on the lower band of the studied frequency interval). 2-tailed distribution was allowed followed by Bonferroni correction for multiple comparisons across the extracted IC states. The reader can refer to (O'Neill et al., 2018(O'Neill et al., , 2017) ) for a more detailed description. Consequently, we will get a set of 𝑁𝑆 𝐶 significant states for HC group and 𝑁𝑆 𝑝 significant states for the PD group.

Between-Group Statistical Differences

As a final stage, we intended to quantify the statistical differences of spatial and temporal features between HC and PD groups. To this end, we projected the approach of EEG microstates, previously used at the sensor-level [START_REF] Khanna | Microstates in resting-state EEG: current status and future directions[END_REF][START_REF] Michel | EEG microstates as a tool for studying the temporal dynamics of whole-brain neuronal networks: A review[END_REF], into the source-level of dynamic Brain Network States (dBNS) as described below and illustrated in Figure 3.

Back-Fitting

First, all selected 𝑁𝑆 𝐶 and 𝑁𝑆 𝑝 states were combined. Then, for each participant (from HC and PD groups), we assigned every temporal window of the individual dFC to the most similar state using the classical microstates back-fitting method [START_REF] Ville | EEG microstate sequences in healthy humans at rest reveal scale-free dynamics[END_REF]. In the following, 𝑁𝑆 𝐶 and 𝑁𝑆 𝑝 states are defined as 'microstates' maps. The spatial similarity (correlation) was calculated between dFC map at every temporal window and each of the microstates maps. Then, using the so-called 'winner-takes-all' algorithm, each temporal window was labeled with the best-fitting microstate map (having the highest spatial similarity). Therefore, we obtained a temporal microstate sequence for each individual dFC for both groups.

Microstates Metrics

With Back-Fitting completed, several features, defined as microstates metrics, can be computed separately for each of the HC and PD subjects [START_REF] Lehmann | EEG microstate duration and syntax in acute, medication-naive, first-episode schizophrenia: a multi-center study[END_REF]. Our outcomes of interest were:

1. Average lifespan or mean duration. The lifespan of a microstate was calculated as the average time duration that a given microstate remains stable for successive segments [START_REF] Lehmann | EEG alpha map series: brain micro-states by space-oriented adaptive segmentation[END_REF]. (Units in seconds).

Fraction coverage time.

The coverage is defined as the ratio of the time frames for which a given microstate is dominant relative to the total recording duration [START_REF] Lehmann | EEG alpha map series: brain micro-states by space-oriented adaptive segmentation[END_REF].

(Units in percentage between 0 and 1).

3. Frequency of occurrence. The frequency of occurrence represents the number of unique appearances of the microstate per second, independently of its duration [START_REF] Lehmann | EEG alpha map series: brain micro-states by space-oriented adaptive segmentation[END_REF]. (Units in Hz).

Global Explained Variance (GEV).

The Global Explained Variance (GEV) of a microstate is the percentage of the total variance explained by this microstate [START_REF] Brodbeck | EEG microstates of wakefulness and NREM sleep[END_REF]. (Units in percentage between 0 and 1).

5. Transition probabilities. The transition probabilities are defined as the sequence of transitions from one microstate to another [START_REF] Lehmann | EEG microstate duration and syntax in acute, medication-naive, first-episode schizophrenia: a multi-center study[END_REF].

Statistical Analysis

To quantify the differences between healthy and PD groups in terms of microstates parameters (mean duration, fraction coverage time, frequency of occurrence, Global Explained Variance, and transition probabilities), statistical tests were performed at the level of each extracted group ICA component. As compared groups have unequal sample sizes (10 HC vs 21 PD), we selected Welsch's t-test to assess the statistical difference between the two groups.

Results

In the following, we denoted Ci the microstate i derived from the HC group, with i= [1; NSc], and Pj the microstate j derived from the PD group, with j= [1; NSp].

Microstates Dynamics -Significant task-modulated dBNS

Results show that in the beta band, one significant dBNS was found in the HC group (NSc = 1) and two in the PD group (NSp = 2). In the gamma band, two significant microstates were derived from the HC group (NSc = 2) and three from the PD group (NSp = 3).

The spatiotemporal dynamics of the several extracted microstates in the beta band and gamma band are illustrated in Figure 4 and Figure 6 respectively. Since we were interested in tracking the evolution of the task-related components, we have plotted on the same time axis, and for each microstate, the group-specific microstates networks and marked the corresponding duration for which they revealed a significant modulation (as described in Materials and Methods section).

For a clear visualization, we plotted the top 0.5% connectivities relative to the total number of unique possible connections, that is, the top 55 edges per network. Using this threshold, in order to globally describe the integrated brain regions and characterize the functional networks for each of the extracted microstates, we calculated the percentage of the number of active nodes with respect to the total number of ROIs, for each of the five microscopic regions (F: Frontal, T: Temporal, P: Parietal, O:

Occipital, and C: Cingulate/insula).

In the beta band, C1 was significantly derived from the HC group at two different times, ranging from 0.13 to 0.46 s (negatively modulated), then, from 0.70 to 1 s (positively modulated). This network mainly involves connections in the frontal lobe (46.43%), the temporal lobe (25%), and the cingulate lobe (28.57%). Concerning the PD group, a microstate P1 was found to have a temporal variation similar to C1 (negative coupling between 0.13 and 0.36 s, followed by positive coupling between 0.58 and 0.87 s). Spatially, connectivities between temporal regions dominate the P1 network (43.48%) mainly in the right cortex, with some connections in the cingulate lobe (30.43%), and few ones in the occipital (21.74%) and parietal lobes (4.75%). Another microstate denoted P2 was positively modulated by the conflict task over the time period 0.39-0.71 s. In this network, connections are distributed over the left hemisphere, including each of temporal (37.04%), frontal (29.63%), cingulate (29.63%), and occipital (3.7%) lobes.

On the other hand, in the gamma band, two microstates were successively identified: C1 from 0.32 to 0.42 s, followed by C2 from 0.56 to 0.72 s. Both were mainly frontal networks with a high contribution of frontal nodes (60.87% in C1 and 62.5% in C2). Temporal regions were also activated, mainly from the right hemisphere for C1 (26.08%), and the left hemisphere for C2 (12.5%), with some cingulate nodes (13.04% for C1 and 25% for C2 For a better investigation of the spatial and temporal results for each microstate, we have displayed the networks from multiple views (top, left, and right) in Figure S1 for beta results, Figure S3 for gamma results. The corresponding trial-averaged time course are also shown along with the built null distribution, with a highlight on the temporal section that surpasses the boundaries of the abovedescribed null distribution. Moreover, to facilitate the understanding of the spatial maps, we highlighted the labels of the active brain regions (among all destrieux ROIs) for each microstate (Figure S2 for beta results, Figure S4 for gamma results). Microstates Analysis -Significant Differences

As described in the Materials and Methods Section, several microstate metrics were computed for both groups and frequency bands, before performing the statistical analysis to quantify significant differences between groups -if any.

Microstates analysis results are presented in Figure 6 for the beta band and Figure 7 for the gamma band. For each of the microstate metrics, including the average lifespan, the frequency of occurrence, the fraction coverage time, and the global explained variance, the mean values (across group subjects) are displayed (colored bars) with the corresponding standard deviation (error bars).

Primarily, for the beta band, we can notice that all microstates metrics relative to P1 were zeros,

showing that P1 was predominated by other microstates, those for C1 were greater among HC participants, contrary to P2 that revealed higher metrics values among PD patients. To detect whether these variations were significant or not, we applied Welsch's t-test. Consequently, significant differences between HC and PD groups were mainly observed in P2 microstate for the following metrics: average lifespan (p-value=0.023), fraction coverage time (p-value=0.019), and GEV (p-value=0.036). A significant difference was also found in C1 for the fraction coverage time (p-value=0.019). The transition probabilities results describe the sequential activity of the obtained microstates. It can be noted that the HC group was found to switch more often between states (C1 and P2), relatively to the PD group that exhibited more stability in the activity of the states, as indicated by the high values of C1-C1 and P2-P2 transitions (Figure 6).

Following the same strategy, microstates analysis results for the gamma band are displayed in Figure 7. It is shown that the P2 microstate had the highest values among all other microstates metrics for both groups. The frequency of occurrence of P2 showed significant differences between HC and PD groups (p-value=0.045). However, none of the other microstates have shown any significant differences among the several metrics (p-value>0.05). From the transition probabilities results, we can see near sequential activity for HC and PD, revealed by an oscillatory activity between P2 and other microstates (Figure 7). representative dynamic states. It was also applied for clinical purposes to characterize the spatiotemporal alterations induced in several brain disorders as in Alzheimer's Disease (AD) (Koelewijn et al., 2017), epilepsy (Koelewijn et al., 2015), and depression (Knyazev et al., 2016;Nugent et al., 2015) using EEG/MEG data during rest and task. Still, to our knowledge, no one has studied the effect of a cognitive (conflict) task on the dynamics of the functional states in the context of Parkinson's Disease and using EEG/MEG techniques, which is the main interest of this work.

Our results show that the cognitive process can be globally described by dominant frontotemporal functional connectivities covering the one-second period that follows the stimulus onset, for most participants and in both frequency bands. Nevertheless, some variability was observed between HC/PD in beta/gamma bands. Specifically, in the beta band, a relevant contribution from the brain regions located in the orbitofrontal cortex (OFC) was identified among HC subjects. In contrast, the intervention of the temporal lobe was stronger in the two PD microstates. Tracking the evolution of these networks reveals the existence of a desynchronization/synchronization process over time.

Primarily, for both groups, the period of 200 to 400 ms after the stimulus offset was characterized by a negative task-modulation, followed by a positive task-modulation for the remaining time. This can be explained by the effect of the stimulus-induced desynchronization followed by a beta rebound during the decision-making (as considered in the event-related (de)synchronization based studies) or as a spontaneous effect of the inhibition/activation process related to the conflict task [START_REF] Panagiotaropoulos | Desynchronization and rebound of beta oscillations during conscious and unconscious local neuronal processing in the macaque lateral prefrontal cortex[END_REF][START_REF] Wu | Attenuated NoGo-related beta desynchronisation and synchronisation in Parkinson's disease revealed by magnetoencephalographic recording[END_REF]. However, understanding the real cognitive mechanism of the obtained temporal fluctuation needs further investigations.

Regarding the gamma band, the microstates of both groups were characterized by strong frontotemporal connections, followed by a more localized activity in the frontal lobe for HC, and more distributed activity over several brain lobes (frontal, temporal, cingulate, occipital) for the PD group.

Temporally, we noticed a delayed microstate activity in PD relative to the HC group (approx. 100ms), which can be related to a slower neuronal activity for PD patients revealed in the high gamma oscillations [START_REF] Johnson | Visual inspection time in Parkinson's disease: deficits in early stages of cognitive processing[END_REF][START_REF] Sawamoto | Cognitive slowing in Parkinson's disease: a behavioral evaluation independent of motor slowing[END_REF].

Beta significant modulation -Association with conflict task

Since we were interested in detecting the fast dynamic reconfiguration of brain network states, we conducted our study within high-frequency bands (beta and gamma). Nevertheless, no significant differences were shown in the gamma band, except for one metric (frequency of occurrence) relative to one network (P2). The main significant differences between HC and PD groups were mainly detected in the beta band. In this frequency band, our microstates analysis findings indicate that functional connections between and within frontal and right temporal areas cover significantly the cognitive activity among HC participants (C1 network, significant duration/coverage/GEV). In contrast to C1, a functional network, with lower activation in the frontal lobe and higher activation in the temporal area, was significantly stable among PD patients (network P2, significant coverage duration).

These results support literature findings, where the prefrontal cortex (PFC) and the orbitofrontal cortex (OFC), were found to be related to cognitive control in the beta frequency range. For example, Stoll et al. have performed a frontal chronic electrocorticography in monkeys, and have demonstrated that frontal beta oscillations are linked to top-down control mechanisms [START_REF] Stoll | The Effects of Cognitive Control and Time on Frontal Beta Oscillations[END_REF]. Spietzer &

Haegens highlighted an increased beta power during the information processing in working memory and decision making in monkeys and humans [START_REF] Spitzer | Beyond the Status Quo: A Role for Beta Oscillations in Endogenous Content (Re)Activation[END_REF]. Furthermore, several studies showed that beta may be observed outside the sensorimotor system, and occurs in PFC during executive control of action and attention [START_REF] Friedman | The role of prefrontal cortex in cognitive control and executive function[END_REF][START_REF] Schmidt | Beta Oscillations in Working Memory, Executive Control of Movement and Thought, and Sensorimotor Function[END_REF][START_REF] Swann | Intracranial EEG Reveals a Time-and Frequency-Specific Role for the Right Inferior Frontal Gyrus and Primary Motor Cortex in Stopping Initiated Responses[END_REF]; for review, see [START_REF] Engel | Beta-band oscillations--signalling the status quo?[END_REF]Wang, 2010).

In addition, the role of temporal regions was also noted in a recent fMRI study using an attention task, where the control attention area located in the temporal lobe has been significantly modulated in macaques and humans [START_REF] Sani | The human endogenous attentional control network includes a ventro-temporal cortical node[END_REF]. Specifically, among these areas, the fusiform region was identified in our HC and PD microstate, and the occipitotemporal sulcus was activated in the PD microstate.

Methodological Considerations and Future Perspectives

The selection of the optimal number of brain network states

Determining the optimal number of derived components is still a challenging question for most decomposition algorithms, including tICA method. In this study, the selection of this parameter was based on two consecutive steps: (1) a primary number of components is selected using the difference in data fitting (DIFFIT) technique adopted by previous studies (Timmerman and Kiers, 2000;Wang et al., 2018). Then, (2) to get reliable results regardless of the estimated number of states, we hypothesized that not all extracted components are necessarily modulated by the task in a significant manner. Following the solution proposed by (Maris and Oostenveld, 2007), a non-parametric null distribution that evolves was generated based on significance testing on the trial-averaged temporal signals.

Here, we adopted a 'sign-flipping' approach using subjects' trials to build surrogate time courses as in previous publications (Hunt et al., 2012;O'Neill et al., 2017;Winkler et al., 2014).

Therefore, only components with trial-averaged time courses that exceed a specific threshold, for a specific continuous duration, were conserved.

Nevertheless, the use of the null distribution may be deemed as too conservative, since we carry out a Bonferroni correction for multiple comparisons, as well as a constraint on the minimal number of significant cycles. For example, we get here only one significant component for HC and two for PD in the beta band. Some would argue that such a restrictive number of components might miss additional information, but in this paper, we focused on studying the alterations induced on the states that are most directly related to the executed task. Still, there is a strong need to seek further reliable and straightforward techniques that automatically and accurately select the significant task-modulated brain states.

The assumption behind microstates and the Back-fitting approach

In the aforementioned microstate analysis procedure, a priori assumption is defined as follows: only one spatial map configuration (microstate) defines the relevant global state of the brain at each time instant. At the back-fitting stage, a spatial correlation was calculated between each microstate map (group-level ICA), and the momentary functional maps (subject-level PLV) are labeled based on the highest correlation [START_REF] Ville | EEG microstate sequences in healthy humans at rest reveal scale-free dynamics[END_REF]. Thus, while the temporal ICA considers in principle overlapping brain network states with proportional weights to summarize the brain activity, the microstate model assumes that all but one of these states is null. Some researchers have argued that the intriguing concept behind this winner-takes-all strategy is based on functional theories that assume that only one global functional state occurs at a given moment in time [START_REF] Baars | The conscious access hypothesis: origins and recent evidence[END_REF][START_REF] Efron | The minimum duration of a perception[END_REF], and remain dominant during a short time period. However, this concept does not deprecate the temporal ICA theory, since it searches for the dominant state, which can correspond to the ICA state with the highest contribution weight. Thus, the zero values obtained for the microstate denoted P1 in the beta band can be explained by the fact that this state exists but is predominated by the remaining states over temporal windows. Future studies can ultimately explore other back-fitting approaches that rely on the notion of proportional, rather than binary, fitting.

Technically, we used the correlation measure to estimate the spatial similarity between the compared networks (without thresholding) at the back-fitting step. We chose to apply this simple -not sophisticated -measure, since we aim to assign each functional network to the nearest (spatially) microstate rather than evaluating the precision of the exact similarity value by itself. Nevertheless, other spatial similarity metrics could be used and exploited in this context, such as those that take into account the spatial locations of the compared networks nodes (Mheich et al., 2018;[START_REF] Pineda-Pardo | Disparate connectivity for structural and functional networks is revealed when physical location of the connected nodes is considered[END_REF], or the distance-based metrics [START_REF] Cao | Natural Sciences Publishing Cor. Measuring Similarity between Graphs Based on the Levenshtein Distance[END_REF][START_REF] Gao | A survey of graph edit distance[END_REF].

Clinical Perspectives

From clinical perspectives, we are currently planning to perform behavioral analysis, including congruence effect measures (reaction time, accuracy), and other cognitive scores to correlate them with the obtained microstates-based metrics relative to each group. This can provide additional electrophysiological quantitative indications useful for neurologists to complement the PD diagnosis based on neuropsychological tests.

As future work, it will be interesting to extend our database and include other groups of PD patients with different levels of cognitive decline [START_REF] Cozac | Quantitative EEG and Cognitive Decline in Parkinson's Disease[END_REF]Hassan et al., 2017) to test the capability of our pipeline methodology to track the dysfunctional alterations that support cognitive impairment, to be ultimately used as an objective predictor of PD dementia.

Conclusion

In conclusion, we reported a study using HD-EEG connectivity states at the source level in PD patients and healthy controls during the Simon task. We showed that PD microstates are characterized by a connectivity reduction in the frontal lobe with an increased temporal area activation. The source-level microstates metrics identified significant differences between HC and PD groups related to the microstates activity in the beta band. We speculate that our strategy could contribute to a better understanding of the cognitive control and ultimate development of EEG-based tests that could help in PD diagnosis. appropriate decomposition method depending on the data and task specifications (dataset size, epoch length, task complexity) when analyzing the dynamics of behavioral tasks. In this context, ICA methods based on high order statistics provided promising results, while SOBI and Kmeans methods exhibited a fragility when applied to rapid and complex data.

 We proposed a novel framework based on a physiologically inspired full brain model to provide a ground-truth for a systematic validation and optimization of the pipeline study including EEGsource space connectivity estimation and dynamic brain networks states extraction. Our findings highlight the benefits and the limitations of several algorithms used to track the dynamics of taskrelated network states in terms of spatial and temporal precision. This framework was used as a proof-of-concept and can be extended to further evaluations to reduce the arbitrary selection of the parameters in the EEG/MEG source space network analysis. Using simulated HD-EEG data, wMNE/PLV/ICA combination has proved a potential capability in deriving the reference dynamic functional states.

 In the context of clinical application, tICA was applied to Parkinson's disease using a conflict task (Simon task). By tracking the spatiotemporal properties of the derived brain network states using source-level microstates approach, we showed the existence of some significant differences between healthy controls and Parkinsonian patients.

Dynamic brain network states: an approach for brain activity representation

Tracking the dynamics of functional connectivity during a cognitive task is crucial in the field of behavioral research and neurological diseases. The dynamic functional connectivity (dFC) is typically estimated between all possible ROIs, within a large number of temporal windows and across the existing trials and subjects. This leads to a set of thousands of connectivity matrices, calling for an automatic summarization approach to help to analyze and interpret dFC results. This approach is usually based on the concept that time-varying functional connectivity is characterized by key recurrent spatiotemporal patterns.

Several approaches were developed to derive the dBNS from the high-dimensional dFC data based on different criteria. For instance, kmeans clustering has been commonly used in several contexts, as fMRI resting-state (Allen et al., 2014;[START_REF] Ciric | Contextual connectivity: A framework for understanding the intrinsic dynamic architecture of large-scale functional brain networks[END_REF], MEG resting-state (O'Neill et al., 2015), and EEG with cognitive tasks (Hassan et al., 2015;Mheich et al., 2015). Here, Kmeans searches for the essential connectivity maps with similar topologies across time (cluster centroids) that are mutually exclusive in time and creates the corresponding probabilistic time course based on clusters occurrences.

On the other hand, component analysis based studies proposed to decompose dFC into temporally overlapping states. For example, the tICA constrains the temporal signatures to be mutually independent. Its ability to elucidate task-induced modulation in MEG analysis was demonstrated in the literature (Brookes et al., 2012a;[START_REF] Luckhoo | Inferring task-related networks using independent component analysis in magnetoencephalography[END_REF]O'Neill et al., 2017). Besides, using PCA after dFC computation, Leonardi et al. have revealed meaningful patterns in FC fluctuations of resting-state networks derived from orthogonal 'eigenconnectivities' decomposition (Leonardi et al., 2013). Matrix

Factorization that imposes the positivity constraint on both connectivity edges and temporal coefficients has also succeeded to uncover transient networks that characterize the neurodevelopment [START_REF] Chai | Evolution of brain network dynamics in neurodevelopment[END_REF].

However, there is no clear evidence about the ideal choice of decomposition techniques. The critical question is whether all these methods with their different constraints detect similar or equivalent spatiotemporal properties of the brain activity? Do they all lead to functionally meaningful states results that give a satisfactory description of the brain activity?

To address these questions, we presented an evaluation analysis of the above-mentioned decomposition techniques to track the dynamics of functional connectivity states. By taking benefit of the excellent temporal resolution of the electrophysiological signals, these algorithms were applied to empirical MEG data (study I) and simulated dense-EEG (study II) data recorded during cognitive tasks.

The key contribution of these comparative analyses is to directly assess the performance of various methods in a novel and unified framework to track the dynamics of task-related brain network states.

Results show some variability between methods performance when applied at group and subject levels. High statistical order methods (ICA) were robust to both data complexity (working memory task) and rapid timescale tasks (order of sub-seconds tasks) more than other tested methods.

Toward EEG network state-based neuromarkers of brain disorders

Neurological pathologies are frequently associated with disruptions in functional brain connectivity.

This allows moving the connectivity analysis from a tool for neuroscientific researches into a novel biomarker for diagnosis. From a clinical perspective, there is a growing need to explore the appropriate methodology to decipher dysfunctionality in the brain connectivity using easy-to-use and non-invasive tools. Neuroimaging tools such as EEG and MEG techniques are of particular interest for the detection of fast neural dynamics with high temporal precision.

behavior in response to a clinical intervention [START_REF] Carbo | Dynamic hub load predicts cognitive decline after resective neurosurgery[END_REF], which offer the opportunity to elaborate novel neuromarker for neurodegenerative diseases.

Methodological considerations

Source Leakage

Mapping electrophysiological data from sensor level to cortical level induces spurious connections between reconstructed adjacent sources. This effect is denoted 'source leakage' (Schoffelen and Gross, 2009). To reduce this effect, some studies proposed to remove edges between very close sources based on a specific definition of inter-node distance (de Pasquale et al., 2016(de Pasquale et al., , 2010)). Others have suggested ignoring zero-lag interactions among signals at the cortical level before performing any connectivity analysis (Brookes et al., 2012b;[START_REF] Colclough | A symmetric multivariate leakage correction for MEG connectomes[END_REF]. However, these approaches may neglect genuine connectivity at zero-lag or among close sources [START_REF] Brookes | Measuring temporal, spectral and spatial changes in electrophysiological brain network connectivity[END_REF][START_REF] Finger | Modeling of Large-Scale Functional Brain Networks Based on Structural Connectivity from DTI: Comparison with EEG Derived Phase Coupling Networks and Evaluation of Alternative Methods along the Modeling Path[END_REF]Schoffelen and Gross, 2009;[START_REF] Singer | Neuronal Synchrony: A Versatile Code for the Definition of Relations?[END_REF]. A classical way to reduce the influence of leakage is to use leakage-invariant functional connectivity metrics as ImCoh metric [START_REF] Nolte | Identifying true brain interaction from EEG data using the imaginary part of coherency[END_REF] and phasebased metrics (Lachaux et al., 1999;[START_REF] Stam | Phase lag index: assessment of functional connectivity from multi channel EEG and MEG with diminished bias from common sources[END_REF]. Another solution is to apply an orthogonalization approach as suggested by [START_REF] Colclough | A symmetric multivariate leakage correction for MEG connectomes[END_REF] when using FC metrics that are sensitive to source leakage effect (refer to (O'Neill et al., 2018) for a review).

It is important to point here that the number of ROIs used may affect the source leakage problem. To date, there is no clear consensus about the optimal choice of ROIs density or distribution (atlas). While a reduced number of large ROIs could reduce the leakage problem, it may damage the spatial resolution. In our studies, we had not limited the choice of ROIs. Instead, we used several brain atlas with distinct resolutions and regions' definitions (78 AAL atlas in Study I, 66 desikan atlas in Study II, and 148 destrieux atlas in Study III), assuming that these distributions were sufficient enough to make a good compromise between spatial resolution and leakage issue.

Inverse Models / Connectivity measures / Dimensionality Reduction methods

A large variety of inverse solutions, functional connectivity metrics, as well as dimensionality reduction techniques, have been proposed and explored in literature. Indeed, choosing the most accurate combination is a challenging issue. To address such a question, we conducted a comparative evaluation analysis in the context of a fast-scale cognitive task in study II. Two inverse models (wMNE and eLORETA) were tested combined with three functional connectivity metrics (PLV, wPLI, and AEC).

wMNE/PLV combination revealed the highest similarity between reference and estimated functional brain networks, in coherence with previous comparative studies in the context of EEG functional connectivity analysis (Hassan et al., 2017b(Hassan et al., , 2014)). Nevertheless, other source reconstruction methods, including beamformer family methods, could be also evaluated along with various existing connectivity metrics, such as coherence, the imaginary part of the coherence (ImCoh), and mutual information. The reader can refer to [START_REF] Pereda | Nonlinear multivariate analysis of neurophysiological signals[END_REF][START_REF] Sakkalis | Review of advanced techniques for the estimation of brain connectivity measured with EEG/MEG[END_REF] for reviews and [START_REF] Wendling | From EEG signals to brain connectivity: a model-based evaluation of interdependence measures[END_REF] for model-based evaluation of the aforementioned metrics.

We should point out that in our MEG analysis (study I), we directly applied the linearly constrained minimum variance (LCMV) (Baker et al., 2014;Brookes et al., 2011b) with the Amplitude Envelope Correlation (AEC) [START_REF] Colclough | How reliable are MEG resting-state connectivity metrics?[END_REF] to estimate functional connectivity, following the work of (Kabbara et al., 2019;O'Neill et al., 2017). Yet, it would be interesting to test the performance of the various (above-cited) inverse models/connectivity families; in the context of MEG dynamic brain networks studies.

Most importantly, the major contribution of this thesis is to track the dynamics of the relevant BNS by the means of the existing decomposition techniques in the context of EEG/MEG analysis. As previously discussed, there exist several approaches to capture the pertinent spatiotemporal information from dFC data and summarize it into a finite set of dynamic -reoccurring -connectivity patterns. However, there is no consensus about the 'ideal' choice of dimensionality reduction method that extracts the 'best' decomposition of dBNS. As an attempt to address this challenge, we conducted a comparative analysis in our studies I and II to evaluate qualitatively and quantitatively the performance of several techniques using real data (study I) and simulated data (study II) respectively. In particular, we focused on ICA, PCA, NMF, and Kmeans clustering approach. Although these different techniques cover a variety of intriguing assumptions used to derive the essential dBNS components, some additional strategies could be also tested simultaneously in a unified robust framework of EEG/MEG analysis. For instance, some studies have demonstrated the ability of Hidden Markov Models (HMM) to elucidate the dynamics of the brain states based on a generative model and Markovian transitioning between states (Baker et al., 2014;Vidaurre et al., 2018). Furthermore, the modularity approach has been used to extract relevant group communities with topological variation behind neural activity [START_REF] Bassett | Robust detection of dynamic community structure in networks[END_REF]Kabbara et al., 2019;[START_REF] Sporns | Modular Brain Networks[END_REF]. The tensor decomposition method has also shown promising results in the field of dynamic brain networks [START_REF] Liu | Network community structure detection for directional neural networks inferred from multichannel multisubject EEG data[END_REF][START_REF] Tang | Temporal-Spatial Patterns in Dynamic Functional Brain Network for Self-Paced Hand Movement[END_REF][START_REF] Williams | Unsupervised Discovery of Demixed, Low-Dimensional Neural Dynamics across Multiple Timescales through Tensor Component Analysis[END_REF]. Indeed, there is a strong need to have a complete EEG/MEG framework, combining most of the existing approaches used to track the evolution of the dominant brain states, since this can bring new methodological insights into the neuroscience field.

Sliding Window approach

The sliding window has been widely applied by the neuroimaging community as a simple and easy-touse technique that illustrate dynamics of functional brain networks related to cognitive abilities [START_REF] Elton | Task-related modulation of functional connectivity variability and its behavioral correlations[END_REF][START_REF] Kucyi | Dynamic functional connectivity of the default mode network tracks daydreaming[END_REF][START_REF] Madhyastha | Age-related differences in the dynamic architecture of intrinsic networks[END_REF] as well as brain disorders (Leonardi et al., 2013;[START_REF] Sakoğlu | A method for evaluating dynamic functional network connectivity and task-modulation: application to schizophrenia[END_REF]. Whilst effective, it carries some critical limitations related to the appropriate selection of window specifications (length, step). While too short windows may decrease specificity through spurious fluctuations in dFC driven by noise (Fraschini et al., 2016;[START_REF] Hutchison | Dynamic functional connectivity: promise, issues, and interpretations[END_REF]Liuzzi et al., 2019), too large windows may decrease sensitivity and fail to detect fast temporal changes of interest.

Therefore, some suggestions have been proposed to reach a satisfactory trade-off between temporal resolution and the accuracy of connectivity matrices. For instance, we followed Lachaux et al.

recommendation in our studies, which proposes to choose the shortest window that assures enough data points to compute phase-based connectivity depending on the frequency of band used (6 cycles) (Lachaux et al., 2000). Still, an accurate quantitative approach to select the temporal width is missing for amplitude-based connectivity metrics. Furthermore, future studies may explore how to dynamically regulate the size of sliding window length relative to the experiment instead of using a fixed window size.

On the other hand, there exists some alternative approaches that escape sliding window constraints.

For instance, Tewarie et al. suggest estimating the dynamics of functional connectivity through hightemporal resolution metrics on a sample-by-sample level instead of windowed aggregated samples in functional connectivity estimation (Tewarie et al., 2019). This approach is of particular interest to detect rapid but genuine connectivity fluctuations as shown by several previous studies [START_REF] Gao | Granger causal time-dependent source connectivity in the somatosensory network[END_REF]Hassan et al., 2015;[START_REF] Martini | The dynamics of EEG gamma responses to unpleasant visual stimuli: from local activity to functional connectivity[END_REF]Tewarie et al., 2019;[START_REF] Wiesman | Oscillatory dynamics and functional connectivity during gating of primary somatosensory responses[END_REF]. Hence, it was applied in Study I with the Human Connectome Project (HCP) dataset since it consists of very short trials, and we suggest deepening the investigation and analysis of these instantaneous functional metrics performance. Besides, time-frequency analysis has been proposed as a novel way to explore temporal dynamics of connectivity at multiple frequencies [START_REF] Chang | Time-frequency dynamics of resting-state brain connectivity measured with fMRI[END_REF]Yaesoubi et al., 2015).

Number of network states selection

Determining the optimal number of derived components is not trivial and still a challenging question for most decomposition algorithms, mainly when working with empirical data. In our studies I and III, the selection of this parameter was based on two techniques used in literature, namely the difference in data fitting (DIFFIT) for decomposition techniques (Timmerman and Kiers, 2000;Wang et al., 2018) and the elbow criterion for clustering methods (Allen et al., 2014). Nevertheless, since there is no consensus about the most accurate optimization criteria, other criteria could be tested such as crossvalidation criterion (Pascual-Marqui et al., 1995), Krzanowski-Lai criterion (Krzanowski and Lai, 1988), and kneedle algorithm (Satopaa et al., 2011). Auxiliary techniques based on model-order selection could be also explored including the Akaike information criterion [START_REF] Mcelreath | Statistical rethinking: a Bayesian course with examples in R and Stan[END_REF] and Bayes information criterion [START_REF] Schwarz | Estimating the Dimension of a Model[END_REF].

It is worthy to note that to get reliable results regardless of the estimated number of states, we hypothesized that not all extracted components are necessarily modulated by the task in a significant manner. Hence, we followed the solution proposed by (Maris and Oostenveld, 2007) that accounts for significant changes in power in the EEG/MEG context. The concept is to execute significance testing on the trial-averaged temporal signals and generate a non-parametric null distribution that evolves. In our work, we adopted a 'sign-flipping' approach using subjects' trials to build surrogate time courses as in previous publications (Hunt et al., 2012;O'Neill et al., 2017;Winkler et al., 2014;Zhu et al., 2020).

Therefore, only components with trial-averaged time courses that exceed a specific threshold were conserved.

However, using non-parametric null distribution is not always convenient since it is built upon the group dataset. First, it may lead to unreliable results when the dataset consists of very few subjects, and it is not compatible with single-subject level studies. On the other hand, computational cost issues may be introduced when dealing with big datasets that consist of too many subjects. Still, there is a strong need to seek further robust techniques that create a straightforward reliable and automatic approach to select significant brain states' number.

Future Directions

In this thesis, we explored the capability of different data-driven dimensionality reduction methods combined with the EEG/MEG source connectivity approach to identify the dynamics of BNS at rapid time scales during cognitive tasks. The pipeline was rigorously evaluated in the context of MEG empirical data, as well as dense-EEG simulated data based on a physiologically inspired computational model, followed by a clinical application on Parkinson's Disease.

In future work, we aim to extend our methodological approach (decomposition techniques) used to derive dBNS into a higher computational level. In particular, we propose to combine two appealing approaches used in literature: (1) time-frequency analysis, and (2) tensor decomposition. In our studies, as in many previous works, we examined our pipeline within a single frequency band. This requires a pre-specification of the frequency band based on the state-of-art of the specific used task.

However, valuable information is failed to be captured given the complex spectral dynamics of electrophysiological activity [START_REF] Chang | Time-frequency dynamics of resting-state brain connectivity measured with fMRI[END_REF]. For instance, cognitive, sensory, motor and emotional tasks may involve different oscillatory rhythms CTC communication. Hence, our methods should be adapted to the specificity of the neuropsychological task. For this reason, we believe that an additional spectral dimension should be added to the parameter space before the application of any decomposition techniques. One of the emerging approaches that allow full characterization of dFC into low dimensional patterns; including spatial (connectivity networks), temporal (dynamic evolution), and spectral (multiple frequency bands) all at once is the tensor decomposition approach (Zhu et al., 2020).

Therefore, this tensorial extension of our methodology can add new insights into the field of dynamic network neuroscience.

In addition, we suggest deepening the evaluation of the high temporal resolution functional connectivity metrics proposed by Tewarie and his colleagues (Tewarie et al., 2019) through a simulation framework based on a biological computational model of the whole-brain (as the one proposed in our study II). This is of particular interest since these measures are sensitive to very fast fluctuations in the temporal brain activity, which can yield a more accurate tracking of the brain dynamics.

Although our codes are openly accessible from the GitHub repository, we got inspired by Iraji et al. (Iraji et al., 2020) and plan to transform our methodology, including source reconstruction, functional connectivity, and dimensionality reduction methods, into a user-friendly interface to allow researchers from different background use/manipulate easily our pipeline with their EEG/MEG datasets.

Finally, we find it interesting to introduce an additional group, known as Mild Cognitive Impairment (MCI) to our clinical application study (study III), besides healthy controls and Parkinsonian patients.

This can help us assess the ability of our tested methodology to uncover the dysfunctional alterations that support cognitive deficits in PD by the means of the extracted dynamic brain network states.

Data and Codes Availability

We believe that open-science can improve the productivity and the quality of the research system by allowing accurate replication and reproducibility of the results. Therefore, for each study, we shared our data/code to the public via GitHub repositories. Abstract : As a complex system, the brain flexibly processes information through dynamic reconfiguration of distributed brain regions at sub-second time scale. A major endeavor in neuroscience is to describe the spatiotemporal organization of the brain as a series of transient "functional connectivity states" using time-resolved analysis. This field is gaining momentum since it not only allows tackling cognitive processes but also holds valuable information about functional alterations of key connectivity patterns in neurological pathologies. In this context, two main challenges have been identified: (1) To what extent can non-invasive neuroimaging techniques with high temporal resolution, namely electro/magnetoencephalography (EEG/MEG), track fast temporally evolving brain states during behavioral tasks? (2) How can neurological diseases affect, spatially and temporally, the identified dynamic brain network states? Therefore, as an attempt to address both challenges, the aim of my thesis is two-folded: 1. Track dynamic brain network states using EEG/MEG. Here the objective is to explore the appropriate methodology that allows extracting relevant connectivity patterns, underlying neural activity when performing tasks. First, three independent MEG datasets from 95 healthy subjects were used during motor and working memory tasks operating on variable time scales. We used the "EEG/MEG source connectivity" method to estimate dynamic functional connectivity (dFC) matrices at the cortical level.

Then, several data-driven decomposition techniques were applied to reduce dFC dimensionality by deriving principal brain patterns with their temporal activation. The performance of these techniques was evaluated and compared at group and subject levels. Second, the previous pipeline was tested using a physiologically based ground truth computational model of a human brain to simulate HD-EEG activity during cognitive task driven at a rapid time scale, as a way to assess a quantitative evaluation of decomposition methods along with multiple key factors used in the pipeline. Primarily, both qualitative and quantitative results show promising outcomes of tested methods with some variability in terms of spatial and temporal accuracy, related to task complexity and time scale. Thus, our findings suggest a careful choice of these methods as they may influence results interpretation. 2. Tracking dysfunctional electrophysiological networks in Parkinson's disease. The main purpose of this work was to identify the major alterations evoked in the extracted dynamic network states for PD patients. For this reason, HD-EEG data was recorded from 31 subjects (21 patients, 10 healthy subjects) during a Simon task. A variant of temporal independent component analysis was used to derive statistically independent components for both groups. Results demonstrate a difference in the spatiotemporal behavior of the dynamic network states between healthy subjects and PD patients

Figure 1 .Figure 2 .

 12 Figure 1. Brain Connectivity Types. Structural connectivity is represented through fiber tracks between brain regions. Functional connectivity reflects the statistical coupling between brain areas while effective connectivity accounts for directionality by measuring the causal influence of one brain region over another............................................................................................................................... 15

  Figure 2. An illustration for the communication through coherence (CTC) (b) and gating through inhibition (GBI) (c) hypotheses. a. Two pools of neurons (A and B) are connected to a pool of neurons C. In this example, pool A communicates with C while functional connectivity between B and C is suppressed. b. In CTC, the phase relationship determines the routing: the neurons in A and C oscillate in phase, whereas the neurons in B do not oscillate in phase with C. The waveforms represent oscillatory population activity in the three regions, whereas the small vertical lines represent spiking activity. c. The magnitude of the pulses of inhibition controls the routing: the flow of information is controlled by an increase of alpha-band oscillations which inhibits firing in pool B, and a decrease in alpha oscillations of neurons in A and C, allowing communication by release from inhibition. This figure is extracted from (Bonnefond et al., 2017). ................................................................................ 17

Figure 3 . 20 Figure 4 . 20 Figure 5 .Figure 6 .Figure 7 .Figure 8 .

 3204205678 Figure 3. Main neuroimaging techniques including EEG, MEG, and fMRI acquisition systems. The corresponding range of spatial and temporal resolutions of the modalities are illustrated (Pfister et al., 2014). ............................................................................................................................................... 20

  /magnétoencéphalographie (EEG/MEG), peuvent-elles suivre l'évolution temporelle rapide des états cérébraux essentiels durant l'exécution d'une tâche? (2) Comment les maladies neurologiques peuvent-elles affecter, spatialement et temporellement, les états dynamiques des réseaux cérébraux? Par conséquent, pour tenter de relever ces deux défis, les deux objectifs de ma thèse sont les suivants : 1. Estimer les états dynamiques des réseaux cérébraux à l'aide des techniques EEG/MEG Le premier objectif consiste à explorer la méthodologie appropriée qui permet d'extraire des motifs de connectivité pertinents relatifs à l'activité neuronale lors de l'exécution d'une tâche. Tout d'abord, trois ensembles indépendants de données MEG chez des sujets sains ont été utilisés pendant des tâches motrice et de mnésique exécutées sur des échelles de temps variables. Nous avons utilisé la méthode de « EEG/MEG source connectivity » suivie d'une estimation dynamique des réseaux fonctionnels afin d'estimer la connectivité fonctionnelle dynamique au niveau cortical. Ensuite, plusieurs techniques de décomposition basées sur les données ont été appliquées pour réduire la dimension des réseaux dynamiques, et ceci en dérivant les principaux états cérébraux avec leur activation temporelle. La performance relative de ces techniques a été évaluée et comparée au niveau du groupe et au niveau individuel. Dans un second temps, une démarche similaire à la précédente a été testée sur des EEG virtuels produits par un modèle computationnel de cerveau humain dans lequel une tâche cognitive de dénomination d'images a été simulée en respectant une échelle de temps très rapide, afin d'évaluer quantitativement les méthodes de décomposition ainsi que certains facteurs clés utilisés. Principalement, les résultats qualitatifs et quantitatifs montrent les effets prometteurs des méthodes testées avec néanmoins une certaine variabilité en termes de précision spatiale et temporelle, liée à la complexité du scénario et à l'échelle temporelle. Cette étude basée sur une vérité terrain indique que le choix des méthodes peut influencer l'interprétation des résultats. 2. Détecter les anomalies de connectivité fonctionnelle au sein des réseaux cognitifs dans la maladie Parkinson L'objectif principal de ce travail était d'identifier les principales altérations dans les états dynamiques des réseaux cérébraux cognitifs chez les patients Parkinsoniens. Pour cette étude, des données EEG de haute résolution (HD-EEG, 256 éléctrodes) ont été enregistrées à partir de 31 sujets (21 patients, 10

Figure 1 .

 1 Figure 1. Brain Connectivity Types. Structural connectivity is represented through fiber tracks between brain regions. Functional connectivity reflects the statistical coupling between brain areas while effective connectivity accounts fordirectionality by measuring the causal influence of one brain region over another.

Figure 3 .

 3 Figure 3. Main neuroimaging techniques including EEG, MEG, and fMRI acquisition systems. The corresponding range of spatial and temporal resolutions of the modalities are illustrated (Pfister et al., 2014).

Figure 4 .

 4 Figure 4. High-Density EEG system (EGI, 256 electrodes) used in Study III

Figure 5 .

 5 Figure 5. Some examples of available brain atlases used to define Regions of Interest (ROIs).

Figure 6 .

 6 Figure 6. Example of binary vs weighted representations of functional connectivity matrices between regions of interest (Ri) (A and B). These matrices can be illustrated by brain networks (C and D), where nodes represent the activated brain regions,and edges represent the functional connectivity values between brain regions.

Figure 7 .

 7 Figure 7. A schematic diagram describing the fundamental processing steps to estimate the cortical functional networks from the recorded electrophysiological data. A. EEG/MEG data is measured at the scalp level. B. Cortical-level signals are obtained by solving the inverse problem using structural MRI and parcellated atlas.C. The statistical coupling between parcellated brain regions is calculated to compute functional connectivity. FC can be either averaged over recording time toobtain the static FC, or computed using a dynamic approach (as sliding window technique) to get a set of dynamic FC networks.

Figure 8 .

 8 Figure 8. The process of dimensionality reduction methods used to extract k dynamic Brain Network States (dBNS) from the estimated dynamic functional connectivity (dFC) denoted 'M'. dBNS represent the decomposed components that define the k dominant spatial maps 'A' with their corresponding temporal evolution 'B'.

  Evaluate and compare the performance of decomposition methods used to estimate electrophysiological-based dBNS during cognitive tasks. (Study I)  Exploit a physiologically based model that simulates the spatiotemporal brain activity underlying a virtual task to provide a quantitative assessment on the various methodological analysis of brain dynamics studies. (Study II)  Characterize the influence of Parkinson's Disease (PD) on spatial and temporal properties of the derived dBNS using a non-invasive technique. (Study III)

Fig. 2 .

 2 Fig. 2. Optimal Number of Components (NC opt ) results. DIFFIT values are plotted against number of component 'J' for all ICA methods, PCA and NMF. The blue plot corresponds to the self-paced button press task, the orange plot for HCP left-hand movement task and the yellow one for Working Memory task. The optimal NC that gives the highest value of DIFIT relative to each task is marked by a small circle on the x-axis. Results of optimal NC relative to Kmeans using the elbow criterion is also shown. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 5 .

 5 Fig. 5. Sternberg working memory task results. Spatial and temporal distribution of all significant components derived from all compared SS methods applied on G-dFC in the working memory task (N 3 = 19subjects). All brain networks are thresholded for visualization. Time values on the x-axis represent the position of the sliding window's center, relative to the first visual stimulus presentation at t = 0sec (as illustrated by a vertical line). The first two vertical lines illustrate the instant of successive visual examples presentation at t = 0s and 1.6sec and the third vertical line at t~9sec separates between the maintenance period that lasts for 7sec and the probe presentation followed by a possible button press and feedback. A color code is attributed for each component in space and time. For each SS method, only significant components (p corrected < 0.05) that appear outside the null distribution are shown here. Exact times of significance relative to each component are indicated with stars. Temporal variation of only JADE, NMF and Kmeans is illustrated, whereas the rest are shown in Supplementary Fig. S7. All NC opt extracted components with corresponding null distribution are shown in Supplementary Fig. S6 for an example of ICA-JADE method. Note that in this task, much larger variety of significant networks are extracted among SS methods, including visual, sensorimotor, language, semantic, and other networks at different temporal activation. In this Fig., 'Vis' refers to Visual network, 'Sem' to Semantic, 'Sens' to Sensorimotor, 'Lang' to Language and'Aux' to other networks. An interactive version of ICA-JADE results can be found on our github https://github.com/judytabbal/dynbrainSS.git using rotatable MATLAB figures.

Fig. 6 .

 6 Fig. 6. Typical example of the spatiotemporal reconfiguration of brain networks during working memory task using ICA-JADE. All significant networks extracted from JADE are collected and presented sequentially relative to each event and period time. The nomination and the exact temporal period of significant activation of each network is clearly indicated. Corresponding cognitive functions are also specified. In this figure, 'Vis' refers to Visual network, 'Sem' to Semantic, 'Sens' to Sensorimotor, 'Lang' to Language and 'Aux' to other networks.

Fig. 8 .

 8 Fig. 8. Subject-level analysis and results relative to the self-paced motor experiment. A. description of subject component selection procedure based on Average Distance (AD) and Correlation Signals (CS) values between subject's component and the group motor component, when SS methods are applied on S-dFC data.Group motor component is shown for PCA example. AD and CS values computed for all components are presented for two subjects. Based on AD and CS condition limits highlighted in each polar bar, success and failure in extracting motor component are both illustrated by subjects 2 and 14 respectively. Spatial and temporal distribution with corresponding AD and CS values for all components of both subjects 2 and 14 are illustrated in Supplementary Fig.S11. A small table on the right illustrates results of success and failure for all subjects in PCA. B. Results of the number of subjects that successfully extracted a motor component in each SS method relative to the total subject's number, denoted Group-Subject Similarity, are shown. C. distributions of AD and CS values of all selected components for the 15 subjects are illustrated. D. Generalization study with increasing number of subjects and the corresponding results of Group-Subject similarity percentage (when considering all possible combinations). E. Grey shaded values (5,6,10 and 14) represent the critical number of subjects required for each SS method to have a 90% precision in extracting the task related component.

  The spatial patterns of connectivity were represented by four (N*N) adjacency matrixPj (j=1,2,3,4), where N represents the number of ROIs. In our case, we used AAL atlas with N=78 ROIs. These matrices represent two networks with a single involved module: visual network (P1) and motor network (P4), one network with two interconnected brain regions: frontotemporal network (P3) and one network with two separated brain regions: occipito-cingulate network (P2). Besides spatial distribution, time evolution of each network was simulated by a combination of a modulation function 𝑓1(𝑡) (Hanning window of unit amplitude) and an uncorrelated Gaussian noise 𝑓2(𝑡): 𝑚 𝑗 (𝑡) = 𝑎. 𝑓1 𝑗 (𝑡) + 𝑏. 𝑓2(𝑡) a and 𝑏 are scalar values set to 0.45 and 0.15. In our study, the total duration is 60 seconds and 't' is sampled at 2 Hz to obtain 120 windows. We set the onset and duration of each network as presented in FigureS2, in a manner to test the ability of each method to successfully extract components with different durations, repeated in time (as P2), overlapped in time (as P2 and P4), and unique in time (as P1 and P3). Finally, the four dynamic network matrices were combined in order to generate single adjacency matrix at each time point 𝑡 over a time-course spanning 60 seconds.𝑆 𝑠𝑖𝑚 (𝑡) = ∑ 𝑚 𝑗 (𝑡)𝑃 𝑗 4 𝑗=1

  Figure S1 shows NCopt results relative to each SS methods. NCopt=4 for SOBI, PSAUD, PCA and NMF in agreement with the exact number of simulated states. Whereas, NCopt=5 for JADE, FastICA and CoM2, and NCopt=6 for InfoMax and Kmeans. Results were then computed relative to these NCopt values.

Figure

  Figure S2.A summarizes the generated simulation showing spatial and temporal of four simulated components. We test the performance of different SS methods in extracting the spatiotemporal characteristics of the four networks (Figure S2.B). Consistency test over multiple runs was also considered by running each algorithm 100 times. We also assessed the computation time and the methods dependency on the noise level. The GS denotes the averaged spatial and temporal similarities over simulated components. The GS values (over 100 runs) are summarized in Figure S2.C. SOBI

Figure S1 .

 S1 Figure S1. Optimal Number of Components results for simulation data. DIFFIT values are plotted against number of components 'J' for ICA, PCA and NMF methods.The optimal NC that gives the highest value of DIFFIT is marked by a small circle on the x-axis. Result of optimal NC relative to Kmeans using the elbow criterion is also shown.

Figure S2 .

 S2 Figure S2. Simulation data analysis and results. A. Simulation diagram describing the four spatial networks (P1-P4) spanning over 60sec and results. B. Example of spatial, temporal and global similarity computation of all simulated components using PCA method for one run. Distributions of global similarity over 100 runs are shown in C. Average elapsed computation time with standard deviation relative to each method are presented in D. The mean value of global similarity is plotted in function of Signal to Noise Ratio (SNR) between 0 to 3dB in E.

Figure S3 .

 S3 Figure S3. Spatial and Temporal Similarity values over 100 runs for all simulated components computed for optimal number of components for each method and at noise level 1.

Figure S8 .

 S8 Figure S8. Sensorimotor network strength proportion for each significant network for all methods in self-paced button press task. Blue, Orange and Yellow colors refer to the corresponding components of Figure 3 in the article. A line is plotted at threshold 0.6 and components with sensorimotor strength ratio above this threshold are denoted 'Mot', otherwise 'Aux'.

Figure S9 .

 S9 Figure S9. Sensorimotor network strength proportion for each significant network for all methods in HCP Left hand motor task. Colors below x-axis refer to the corresponding components of Figure 4 in the article. A line is plotted at threshold 0.6 and components with sensorimotor strength ratio above this threshold are denoted 'Mot', otherwise 'Aux'.

Figure S10 .

 S10 Figure S10. Relevant networks masks('Vis', 'Lang', 'Sem' and 'Sens') strength proportion for each significant network for all methods in working memory task. Bars colors refer to corresponding network masks as indicated on the top of the figure. Colors below x-axis refer to the corresponding components of Figure5in the article. A line is plotted at threshold 0.6 and components with a network mask strength ratio above this threshold are denoted relative to the mask.

Figure S11 .

 S11 Figure S11. Spatial networks and temporal variation results of all extracted components for PCA method applied on subjectlevel data in the self-paced motor experiment, for two subjects (2 and 14) that succeeded and failed, respectively, to extract motor component. Corresponding Average Distance (AD) and Correlation Signals (CS) values are indicated.

into 20 ,

 20 20 and 21 subjects and performed separately three combinations, two consisting of 𝐶 10 20 = 184756 repetitions and one of 𝐶 10 21 = 352716 repetitions. This null distribution can find out the existence of trial-onset-locked increases or decreases in connectivity. Therefore, if the average temporal signal of a component fell outside the null distribution at a specific time point, we consider a significant modulation of this component at this time point (O'Neill et al., 2017). Temporal signals were normalized for both NMF and Kmeans methods before applying the null distribution as their temporal weights are always positive. Threshold values were determined relative to NCopt and temporal degree of freedom. It is assumed that a single temporal degree of freedom ('dof') was added each time the window shifts by more than half of its width. This meant a total of 8 dof in the self-paced data and 12 in the Sternberg data. The 'dof' applied in HCP was the same as in the self-paced motor task. Thresholds were therefore set at (0.05/(2*NCopt*dof)).

Figure 1 .

 1 Figure 1. Structure of the investigation. (1) COALIA Model was used to simulate cortical sources from a set of simulated networks. (2) 257 scalp level signals were generated using the forward model. (3) Cortical-level sources (66 desikan regions) were reconstructed using wMNE and eLORETA methods. (4) Dynamic Functional Connectivity between sources was calculated using PLV, wPLI, and AEC methods. (5) Dimensionality reduction methods including ICA, PCA, NMF, and Kmeans were applied to reconstruct dominant brain networks states. (6) Spatial and Temporal similarities were finally calculated to measure the performance of brain networks estimation.

  model named COALIA. Briefly, this model is based on interconnected Neural Mass Model (NMM) respecting human structural connectivity based on DTI, which has proven its ability to produce realistic EEG as compared to real EEG recordings obtained in humans during awake and deep sleep[START_REF] Bensaid | COALIA: A Computational Model of Human EEG for Consciousness Research[END_REF].Briefly, each microcircuit based on neural mass can produce distinct brain oscillations such as alpharhythms through Pyramidal somatostatin/positive (SST) loop, gamma-rhythms through Pyrparvalbumin (PV) loop, delta-rhythms through increased thalamocortical connectivity and disinhibition through VIP-SST microcircuits. Therefore, using distinct local oscillations cortical-level simulated signals taking into account macro and micro-circuitry of the human cortex, it becomes possible to generate coherent oscillations in different brain structures to simulate dynamic functional connectivity during a cognitive task. (i) At the local level, each neural mass consists of subsets of a glutamatergic pyramidal neuron population and three types of GABAergic interneuron populations (VIP-, PV-, SSTinterneurons) with physiologically based kinetics (fast/slow) and interconnection. Increased excitability of the perisomatic targeting PV-interneurons onto pyramidal neurons produced gamma oscillations (30-45Hz). Increased excitability of the dendritic targeting SSTinterneurons onto pyramidal neurons produced alpha oscillations (8-12Hz). Remote activation of the VIP-interneurons onto SST induced disinhibition and strengthened coherent gamma oscillations between coactivated interconnected cortical areas. (ii) Then, at the global level, the large-scale model is constructed based on nROIs=66 regions of interest from the standard anatomical parcellation of Desikan-Killiany atlas (Desikan et al., 2006) (right and left insula were excluded, leaving 66 brain regions). In this case, each neural mass represents the local field potential (LFP) of one atlas region, in which the activity is assumed to be homogeneous. The template brain morphology (Colin) is used to spatially distribute neural masses over the cortex.

Figure 2 .

 2 Figure 2. Simulation Scenario. All activated ROIs are shown for each time interval (from T1 to T6). Different color codes refer to different lobe affiliations. Orange color stands for occipital lobe, blue for temporal lobe, pink for central lobe, red for prefrontal lobe, and brown for limbic lobe.

Figure 3 .

 3 Figure3. Inverse model evaluation (wMNE and eLORETA) using different source space resolutions. LowRes refers to lead field computation on 66 desikan regions and HighRes refers to lead field computation on high-resolution cortex 15002 followed by projection on 66 desikan regions. Reconstructed brain sources are shown at the group level (averaged over trials and subjects). The top 10 activated regions are visualized for each method and at each time interval. Brain regions' colors are based on their weights (blue for the less activated, red for the most activated). The precision value was calculated relative to LFP reference sources at each time interval with a proportional threshold ranging from 10% to 15% of the total node number. The precision value is shown on the right side of each method as the average value over all time intervals and threshold values.

Figure 4 .

 4 Figure 4. Inverse model/functional connectivity combination evaluation. Reconstructed brain networks are shown at the group level (averaged over trials and subjects). All brain networks were thresholded for visualization by keeping the top 1% edges

Figure 6 .

 6 Figure 6. Group-level and Subject-level overall methodology evaluation. Spatial similarity averaged over all time intervals was used as an overall evaluation measure for different inverse model/functional connectivity combinations at the group-level (A.) and subject-level (C.). Maximal Global similarity averaged over all time intervals was used as an overall evaluation measure for different dimensionality reduction methods at the group-level (B.) and subject-level (D.). Each method is represented by a specific color. Circular bars and boxplots were used to visualize results at both group and subject levels respectively. Results of each subject are represented by small circles on the left side of each boxplot. We also overlaid results obtained in A. and B. in the boxplots in C. and D. respectively as centered bold circles to show subject-level results relative to the group level.

Figure S5 .

 S5 Figure S5. Noise variation effect on dimensionality reduction methods performance. The maximal global similarity is computed at the group-level and averaged over all time intervals for each noise level represented by λ value ranged from 0.7 (for the most noisy data) to 1 (for the least noisy data). Colors refer to different dimensionality reduction method.

Figure 1 .

 1 Figure 1. Overview of the Simon Task used. The central fixation point was displayed randomly from 1750 to 2170 ms, with a 30 ms step. Then, the stimulus display lasted 200 ms. Participants had 1000 ms to answer by pressing the button. Two conditions could occur: congruent, when color and location of the circle led to the same answer; and incongruent, when color and location didn't lead to the same answer.

Figure 2 .

 2 Figure2. Outline of the dynamic FC pipeline. First, HD-EEG data was recorded during the Simon task (only correct incongruent trials were considered). After preprocessing, cortical-level sources were reconstructed using the weighted Minimum Norm Estimate (wMNE) and the Destrieux atlas (148 Regions of Interests ROIs). Then, the dynamic functional connectivity (dFC) was estimated for each subject and trial using the sliding window approach. Phase Locking Value measure (PLV) was used to compute the statistical coupling between ROIs. Finally, the temporal Independent Component Analysis (tICA) was applied on

Figure 39

 39 Figure39. A. Temporal ICA Applied on HC and PD groups. The resultant significant NSc and NSp ICA networks for HC and PD groups respectively are illustrated by specific colors. B. Back-Fitting approach assigns PLV network at each temporal window with the group ICA network (microstate) having the highest spatial similarity value. This was applied on all HC and PD subjects, followed by a calculation of the main microstates metrics as shown in C.

Figure 4 .Figure 5 .

 45 Figure 4. Spatiotemporal dynamics of the significant task-modulated dBNS for both HC(A.) and PD(B.) groups in the beta band. The time 0 s corresponds to the stimulus offset. Each significant dBNS is illustrated as a brain network with a specific color. All brain networks were thresholded for visualization. Spheres of different sizes proportional to their strength represent the activated brain nodes. A color code is attributed for all nodes belonging to the same brain lobe (yellow for frontal, blue for temporal, light blue for parietal, green for occipital, and red for cingular and insula). For each state, we indicated the temporal duration on which it is significantly modulated by the task (positively modulated are plotted above the time axis, negatively modulated are plotted below the time axis). On the right side, the percentages of nodes relative to each brain lobe are illustrated on the colored bars for each state. The reader can refer to Supplementary FigureS1for a detailed view representation of the network (top, left, and right) with the corresponding averaged-trial temporal signals plotted over the whole temporal duration along with the null distribution to reveal the temporal significance. In the supplementary FigureS3, the labels of the activated destrieux ROIs are highlighted.

Figure 6 .

 6 Figure 6. Microstates parameters and statistical analysis between HC and PD groups in the beta band. The average lifespan, the frequency of occurrence, the fraction coverage time, and the global explained variance results are represented by colored bars (mean value across subjects) (blue for HC group and orange for PD group) in A, B, C, and D. The standard deviations of the corresponding metrics are displayed as error bars. A red asterisk with the corresponding p-value illustrates the presence of statistical differences(Welsch's Test;) between the two groups in terms of microstates parameters. The transition probabilities between all microstates are also shown in E.

Figure 7 .

 7 Figure 7. Microstates parameters and statistical analysis between HC and PD groups in the gamma band. The average lifespan, the frequency of occurrence, the fraction coverage time, and the global explained variance results are represented by colored bars (mean value across subjects) (blue for HC group and orange for PD group) in A, B, C, and D. The standard deviations of the corresponding metrics are displayed as error bars. A red asterisk with the corresponding p-value illustrates the presence of statistical differences(Welsch's Test;) between the two groups in terms of microstates parameters. The transition probabilities between all microstates are also shown in E.

Figure S2 .

 S2 Figure S2. Top, Left and Right views of the extracted microstates networks with the corresponding trial-averaged signal in the gamma band. The highlighted section of the null distribution indicates the time when the microstate surpasses the limits of the null distribution. C1 refers to the microstate extracted from HC group, P1 and P2 to those derived from PD group.

Le

  cerveau humain est un réseau complexe à grande échelle où les régions cérébrales sont hautement interconnectées et distribuées. Pour mieux comprendre les systèmes neuronaux, la neuroscience des réseaux ('Network Neuroscience') a fourni des outils nécessaires qui étudient et caractérisent les réseaux cérébraux. Ceci a contribué à des découvertes cruciales liées à notre compréhension de la fonctionnalité du cerveau chez les humains sains et malades. Les études ont montré que l'organisation fonctionnelle du réseau cérébral se reconfigure d'une manière flexible et dynamique à une échelle de temps inférieure à la seconde pour maintenir un contrôle efficace en temps réel de l'activité cérébrale au repos et durant des tâches cognitives/comportementales. De plus, il est désormais reconnu que les maladies neurodégénératives sont caractérisées par des altérations fonctionnelles dynamiques des réseaux cérébraux, considérées ainsi comme 'maladies de réseau'. Du point de vue clinique, il existe donc une forte demande pour de nouveaux outils de diagnostic ou 'neuromarker' basés sur des techniques non invasives, et des méthodes permettant d'identifier les réseaux pathologiques. Par conséquent, une caractérisation systématique de la dynamique des réseaux cérébraux fonctionnels sains et pathologiques est cruciale dans les domaines des neurosciences et cliniques. Cela nécessite l'utilisation de techniques de neuroimagerie non invasives avec une excellente résolution temporelle comme l'électro/magnétoencéphalographie « EEG/MEG », ainsi que l'application de méthodes avancées d'analyse de données capables de suivre les fluctuations spatio-temporelles très rapide des réseaux cérébraux dominants.Des études antérieures sur les analyses de la connectivité fonctionnelle à partir de l'EEG/MEG ont été principalement réalisées au niveau des électrodes (scalp-level). Cependant, l'interprétation des réseaux EEG obtenus au niveau des électrodes n'est pas simple, car les signaux sont altérés par le problème de « volume de conduction ». Ainsi, pour réduire l'effet de ce problème, la méthode de « connectivité de source EEG/MEG » a été proposée pour étudier de la connectivité fonctionnelle du cerveau large échelle avec une grande précision temporelle. Ceci permet une reconstruction fiable des sources EEG/MEG responsables de la génération des signaux EEG au niveau des électrodes.Par ailleurs, un objectif majeur en neurosciences est de décrire l'organisation spatio-temporelle du cerveau comme une série d'« états de connectivité fonctionnelle » transitoires à travers une analyse dynamique des réseaux. Par conséquent, le cerveau peut occuper plusieurs états au fil du temps, avec une topologie de réseau spécifique relative à chaque état. Ce domaine prend de l'ampleur car il permet non seulement d'aborder les processus cognitifs, mais aussi d'apporter des informations importantes sur les altérations fonctionnelles des principaux motifs de connectivité dans le cadre des pathologies neurologiques.C'est dans ce contexte que s'inscrivent les travaux de cette thèse qui prolongent les développements méthodologiques et cliniques de notre équipe de recherche sur la dynamique de la connectivité fonctionnelle au niveau des sources cérébrales. Dans cette thèse, nous avons exploré les performances de l'analyse dynamique de la méthode "connectivité de source EEG/MEG" suivie de techniques de décomposition pour élucider les «états dynamiques du réseau cérébral» élémentaires qui changent au cours du temps et dominent l'activité cérébrale lors des tâches cognitives. Bien que des efforts considérables aient été récemment consacrés pour suivre la séquence d'états évoluant rapidement dans le temps, i) une évaluation authentique des différentes approches analytiques existantes d'une part, et ii) une exploration de l'effet des pathologies sur l'architecture des états dynamiques du réseau cérébral d'autre part, restent manquantes. Par conséquent, l'objectif de mes travaux de thèse est double : 1) progresser sur les aspects méthodologiques des techniques de décomposition appliquées sur la connectivité de sources en EEG/MEG et 2) utiliser cette méthodologie dans une application clinique en lien avec la maladie de Parkinson. Dans la première partie (aspects méthodologiques), des études comparatives entre neuf différents méthodes de décomposition sont effectuées afin d'évaluer les capacités et les limitations de chaque méthode à extraire les états dynamiques dominants dans le contexte de plusieurs taches cognitives rapides en utilisant des données électrophysiologiques empiriques (première étude) et simulées (deuxième étude). Dans la deuxième partie (application clinique), les altérations dans les états dynamiques des réseaux cérébraux chez des Parkinsoniens sont explorées (troisième étude). Aspects méthodologiques Le premier objectif consiste à explorer la méthodologie appropriée qui permet d'extraire des motifs de connectivité pertinents relatifs à l'activité neuronale lors de l'exécution d'une tâche. Dans ce contexte, des évaluations qualitatives et quantitatives ont été exécutées aux niveau groupe et individuel comme suit. Tout d'abord, trois ensembles indépendants de données MEG chez des sujets sains ont été utilisés pendant des tâches motrice et de mnésique exécutées sur des échelles de temps variables (première étude). Nous avons utilisé la méthode de « connectivité de source EEG/MEG » suivie d'une estimation dynamique des réseaux fonctionnels afin d'estimer la connectivité fonctionnelle dynamique au niveau cortical. Ensuite, plusieurs techniques de décomposition basées sur les données ont été appliquées pour réduire la dimension des réseaux dynamiques, et ceci en dérivant les principaux états cérébraux avec leur activation temporelle. La performance relative de ces techniques (fonctionnalités et limites) a été évaluée et comparée au niveau du groupe et au niveau individuel. Principalement, les résultats montrent les effets prometteurs des méthodes testées avec néanmoins une certaine variabilité en termes de précision spatiale et temporelle, liée à la complexité du scénario et à l'échelle temporelle. Parmi les méthodes de décomposition/regroupement testées, on peut voir que l'indépendance statistique (d'ordre élevé) décrit le mieux l'activité cérébrale temporelle des états dynamiques des réseaux les plus représentatifs au niveau du groupe et du sujet. Dans ce contexte, les méthodes ICA basées sur des statistiques d'ordre élevé ont fourni des résultats prometteurs, tandis que les méthodes SOBI et Kmeans ont montré une fragilité dans leur précision spatiotemporelle lorsqu'elles sont appliquées à des données rapides et complexes. Notre message principal est que les chercheurs doivent être conscients de sélectionner la méthode de décomposition appropriée en fonction des données et des spécifications de la tâche (taille de l'ensemble de données, durée des essais, complexité de la tâche) lors de l'analyse de la dynamique des tâches comportementales. Dans un second temps, une démarche similaire à la précédente a été testée sur des EEG virtuels produits par un modèle computationnel de cerveau humain dans lequel une tâche cognitive de dénomination d'images a été simulée en respectant une échelle de temps très rapide, afin d'évaluer quantitativement les méthodes de décomposition ainsi que certains facteurs clés utilisés (deuxième étude). Dans cette étude, nous avons proposé un nouveau cadre basé sur un modèle cérébral complet inspiré de la physiologie afin de fournir une vérité de terrain pour une validation et une optimisation systématiques de l'étude du pipeline, y compris l'estimation de la connectivité spatiale de la source EEG et l'extraction des états dynamiques des réseaux cérébraux. Nos résultats mettent en évidence les avantages et les limites de plusieurs algorithmes utilisés pour suivre la dynamique des états de réseau liés aux tâches en termes de précision spatiale et temporelle. Cette étude basée sur une vérité terrain indique que le choix des méthodes peut influencer l'interprétation des résultats et peut être aussi étendue à d'autres évaluations pour réduire la sélection arbitraire des paramètres dans l'analyse du réseau spatial source EEG/MEG. Application Clinique L'objectif principal de ce travail était d'identifier les principales altérations dans les états dynamiques des réseaux cérébraux cognitifs chez les patients Parkinsoniens (troisième étude). Dans le cadre de cette application clinique, des données EEG de haute résolution (HD-EEG, 256 éléctrodes) ont été enregistrées à partir de 31 sujets (21 patients, 10 sujets sains) au cours de la tâche de conflit cognitif nommée Simon-Task. Une variante de l'analyse des composantes indépendantes a été utilisée pour dériver des composantes statistiquement indépendantes dans les deux groupes. Les résultats démontrent l'existence de différences spatiotemporelles significatives dans les états dynamiques des réseaux cérébraux entre les sujets sains et les patients. En particulier, la connectivité beta préfrontale est quasi absente chez les patients parkinsoniens. Cette étude souligne que la connectivité fonctionnelle dynamique basée sur les tâches est une approche prometteuse pour comprendre les dysfonctionnements cognitifs observés dans la maladie de Parkinson et au-delà. Titre : Dynamique des réseaux cérébraux électrophysiologiques (dys)fonctionnels Mots clés : Neurosciences, Réseaux cérébraux, Maladies neurodégénératives, Connectivité fonctionnelle, Séparation de sources, Electro-Magnétoencéphalographie (EEG/MEG) Résumé : En tant que système complexe, le cerveau traite de manière flexible les informations grâce à une reconfiguration dynamique des réseaux neuronaux sur une échelle de temps de l'ordre de la milliseconde. Un objectif majeur en neurosciences est de décrire l'organisation spatio-temporelle du cerveau comme une série d'«états de connectivité fonctionnelle » transitoires à travers une analyse dynamique des réseaux. Ce domaine prend de l'ampleur car il permet non seulement d'aborder les processus cognitifs, mais aussi d'apporter des informations importantes sur les altérations fonctionnelles des principaux motifs de connectivité dans le cadre des pathologies neurologiques. Dans ce contexte, deux enjeux principaux ont été identifiés : (1) A quel point les techniques de neuroimagerie non-invasives à haute résolution temporelle, tel que l'électro/magnétoencéphalographie (EEG/ MEG), peuvent-elles suivre l'évolution temporelle rapide des états cérébraux essentiels durant l'exécution d'une tâche? (2) Comment les maladies neurologiques peuvent-elles affecter, spatialement et temporellement, les états dynamiques des réseaux cérébraux? Par conséquent, pour tenter de relever ces deux défis, les deux objectifs les états dynamiques des réseaux cérébraux à l'aide des techniques EEG/MEG. Le premier objectif consiste à explorer la méthodologie appropriée qui permet d'extraire des motifs de connectivité pertinents relatifs à l'activité neuronale lors de l'exécution d'une tâche. Tout d'abord, trois ensembles indépendants de données MEG chez des sujets sains ont été utilisés pendant des tâches motrice et de mnésique exécutées sur des échelles de temps variables. Nous avons utilisé la méthode de « EEG/MEG source connectivity » suivie d'une estimation dynamique des réseaux fonctionnels afin d'estimer la connectivité fonctionnelle dynamique au niveau cortical. Ensuite, plusieurs techniques de décomposition basées sur les données ont été appliquées pour réduire la dimension des réseaux dynamiques, et ceci en dérivant les principaux états cérébraux avec leur activation temporelle. La performance relative de ces techniques a été évaluée et comparée au niveau du groupe et au niveau individuel. Dans un second temps, une démarche similaire à la précédente a été testée sur des EEG virtuels produits par un modèle computationnel de cerveau humain dans lequel une tâche cognitive de dénomination d'images a été simulée en respectant une échelle de temps très rapide, afin d'évaluer quantitativement les méthodes de décomposition ainsi que certains facteurs clés utilisés. Principalement, les résultats qualitatifs et quantitatifs montrent les effets prometteurs des méthodes testées avec néanmoins une certaine variabilité en termes de précision spatiale et temporelle, liée à la complexité du scénario et à l'échelle temporelle. Cette étude basée sur une vérité terrain indique que le choix des méthodes peut influencer l'interprétation des résultats. 2. Détecter les anomalies de connectivité fonctionnelle au sein des réseaux cognitifs dans la maladie Parkinson. L'objectif principal de ce travail était d'identifier les principales altérations dans les états dynamiques des réseaux cérébraux cognitifs chez les patients Parkinsoniens. Pour cette étude, des données EEG de haute résolution (HD-EEG, 256 éléctrodes) ont été enregistrées à partir de 31 sujets (21 patients, 10 sujets sains) au cours de la tâche de conflit cognitif nommée Simon-Task. Une variante de l'analyse des composantes indépendantes a été utilisée pour dériver des composantes statistiquement indépendantes dans les deux groupes. Les résultats démontrent l'existence de différences spatiotemporelles dans les états dynamiques des réseaux cérébraux entre les sujets sains et les patients. Title : Dynamics of electrophysiological (dys)functional brain networks Keywords : Neurosciences, Brain networks, Neurodegenerative diseases, Functional connectivity, Sources Separation, Electro-Magnétoencéphalography (EEG/MEG)

3. RESULTS Study I: Dynamics of task-related electrophysiological networks: a benchmarking study

	Judie Tabbal, Aya Kabbara, Mohamad Khalil, Pascal Benquet, Mahmoud Hassan.
	Published in: Neuroimage journal (2021) 231: 117829
	Highlights
		We explore the fast reconfiguration of electrophysiological brain networks during motor and
		working memory tasks,
		We provide a quantitative evaluation of the performance of nine source separation methods
		applied on dynamic functional connectivity,
		Results show variability between the methods at both group and subject levels in terms of
		space/time accuracy,
		Independent Component Analysis (ICA) methods based on high order statistics provide promising
		results, while SOBI and Kmeans exhibit a fragility related to data complexity and timescale
		resolution.

  Aggarwal, C.C.,Hinneburg, A., Keim, D.A., 2001. On the surprising behavior of distance metrics in high dimensional space. In: Van den Bussche, J., Vianu, V. (Eds.), Database Theory -ICDT 2001, Lecture Notes in Computer Science. Springer, Berlin, Heidelberg, pp. 420-434. doi: 10.1007/3-540-44503-X_27 . Allen, E.A., Damaraju, E., Plis, S.M., Erhardt, E.B., Eichele, T.,Calhoun, V.D., 2014. Track- ing whole-brain connectivity dynamics in the resting state. Cereb Cortex 24, 663-676. doi: 10.1093/cercor/bhs352 . Axmacher, N., Schmitz, D.P., Wagner, T., Elger, C.E.,Fell, J., 2008. Interactions between medial temporal lobe, prefrontal cortex, and inferior temporal regions during visual working memory: a combined intracranial eeg and functional magnetic resonance imaging study. J. Neurosci. 28, 7304-7312. doi: 10.1523/JNEUROSCI.1778-08.2008 . Baddeley, A., Jarrold, C., Vargha-Khadem, F., 2011. Working memory and the hippocampus. J. Cogn. Neurosci. 23, 3855-3861. doi: 10.1162/jocn_a_00066 . Baker, A.P., Brookes, M.J., Rezek, I.A., Smith, S.M., Behrens, T., Probert Smith, P.J., Woolrich, M., 2014. Fast transient networks in spontaneous human brain activity. eLife 3, e01867. doi: 10.7554/eLife.01867 . Bassett, D.S. , Sporns, O. , 2017. Network neuroscience. Nature Neurosci. 20, 353 . Bauer, A.-K.R., Debener, S., Nobre, A.C., 2020. Synchronisation of neural oscillations and cross-modal Influences. Trends Cognit. Sci. 24, 481-495. doi: 10.1016/j.tics.2020.03.003 . Bola, M. , Sabel, B.A. , 2015. Dynamic reorganization of brain functional networks during cognition. Neuroimage 114, 398-413 . Brookes, M.J., Hale, J.R., Zumer, J.M., Stevenson, C.M., Francis, S.T., Barnes, G.R., Owen, J.P., Morris, P.G., Nagarajan, S.S., 2011a. Measuring functional connectivity using MEG: methodology and comparison with fcMRI. Neuroimage 56, 1082-1104. doi: 10.1016/j.neuroimage.2011.02.054 . Brookes, M.J., Tewarie, P.K., Hunt, B.A.E., Robson, S.E., Gascoyne, L.E., Liddle, E.B., Liddle, P.F., Morris, P.G., 2016. A multi-layer network approach to MEG connectivity analysis. NeuroImage 132, 425-438. doi: 10.1016/j.neuroimage.2016.02.045 . Brookes, M.J., Vrba, J., Robinson, S.E., Stevenson, C.M., Peters, A.M., Barnes, G.R., Hillebrand, A., Morris, P.G., 2008. Optimising experimental design for MEG beamformer imaging. NeuroImage 39, 1788-1802. doi: 10.1016/j.neuroimage.2007.09.050 . Brookes, M.J., Wood, J.R., Stevenson, C.M., Zumer, J.M., White, T.P., Liddle, P.F., Morris, P.G., 2011b. Changes in brain network activity during working memory tasks: a magnetoencephalography study. Neuroimage 55, 1804-1815. doi: 10.1016/j.neuroimage.2010.10.074 . Brookes, M.J., Woolrich, M.W., Barnes, G.R., 2012. Measuring functional connectivity in MEG: a multivariate approach insensitive to linear source leakage. Neuroimage 63, 910-920. doi: 10.1016/j.neuroimage.2012.03.048 . Buckner, R.L., Krienen, F.M., Castellanos, A., Diaz, J.C., Yeo, B.T.T., 2011. The organization of the human cerebellum estimated by intrinsic functional connectivity. J Neurophysiol 106, 2322-2345. doi: 10.1152/jn.00339.2011 . Bullmore, E., Sporns, O., 2009. Complex brain networks: graph theoretical analysis of structural and functional systems. Nat. Rev. Neurosci. 10, 186-198. doi: 10.1038/nrn2575 . Caminiti, S.P., Siri, C., Guidi, L., Antonini, A., Perani, D., 2015. The neural correlates of spatial and object working memory in elderly and Parkinson's disease subjects. Behav. Neurol. 2015, 1-10. doi: 10.1155/2015/123636 . Casorso, J., Kong, X., Chi, W., Van De Ville, D., Yeo, B.T.T., Liégeois, R., 2019. Dynamic mode decomposition of resting-state and task fMRI. NeuroImage 194, 42-54. doi: 10.1016/j.neuroimage.2019.03.019 . Cavanna, A.E., Trimble, M.R., 2006. The precuneus: a review of its functional anatomy and behavioural correlates. Brain 129, 564-583. doi: 10.1093/brain/awl004 . Chai, L.R., Khambhati, A.N., Ciric, R., Moore, T.M., Gur, R.C., Gur, R.E., Satterthwaite, T.D., Bassett, D.S., 2017. Evolution of brain network dynamics in neurodevelopment. Network Neurosci. 1, 14-30. doi: 10.1162/NETN_a_00001 . Chai, W.J., Abd Hamid, A.I., Abdullah, J.M., 2018. Working memory from the psychological and neurosciences perspectives: a review. Front Psychol. 9. doi: 10.3389/fpsyg.2018.00401 . Ciric, R., Nomi, J.S., Uddin, L.Q., Satpute, A.B., 2017. Contextual connectivity: a framework for understanding the intrinsic dynamic architecture of large-scale functional brain networks. Sci. Reports 7, 1-16. doi: 10.1038/s41598-017-06866-w . Colclough, G.L., Brookes, M.J., Smith, S.M., Woolrich, M.W., 2015. A symmetric multivariate leakage correction for MEG connectomes. Neuroimage 117, 439-448. Timmerman, M.E., Kiers, H.A., 2000. Three-mode principal components analysis: choosing the numbers of components and sensitivity to local optima. Br. J. Math. Stat. Psychol. 53 (Pt 1), 1-16. doi: 10.1348/000711000159132 . Van Essen, D.C., Ugurbil, K., Auerbach, E., Barch, D., Behrens, T.E.J., Bucholz, R., Chang, A., Chen, L., Corbetta, M., Curtiss, S.W., Della Penna, S., Feinberg, D., Glasser, M.F., Harel, N., Heath, A.C., Larson-Prior, L., Marcus, D., Michalareas, G., Moeller, S., Oostenveld, R., Petersen, S.E., Prior, F., Schlaggar, B.L., Smith, S.M., Snyder, A.Z., Xu, J., Yacoub, E.WU-Minn HCP Consortium, 2012. The human connectome project: a data acquisition perspective. Neuroimage 62, 2222-2231. doi: 10.1016/j.neuroimage.2012.02.018 . Vandenberghe, R., Price, C., Wise, R., Josephs, O., Frackowiak, R.S.J., 1996. Functional anatomy of a common semantic system for words and pictures. Nature 383, 254-256. doi: 10.1038/383254a0 . Vidaurre, D., Abeysuriya, R., Becker, R., Quinn, A.J., Alfaro-Almagro, F., Smith, S.M., Woolrich, M.W., 2018a. Discovering dynamic brain networks from big data in rest and task. NeuroImage, Brain Connect. Dynamics 180, 646-656. doi: 10.1016/j.neuroimage.2017.06.077 . Vidaurre, D., Hunt, L.T., Quinn, A.J., Hunt, B.A.E., Brookes, M.J., Nobre, A.C., Woolrich, M.W., 2018b. Spontaneous cortical activity transiently organises into frequency specific phase-coupling networks. Nature Commun. 9, 2987. doi: 10.1038/s41467-018-05316-z . Vidaurre, D., Quinn, A.J., Baker, A.P., Dupret, D., Tejero-Cantero, A., Woolrich, M.W., 2016. Spectrally resolved fast transient brain states in electrophysiological data. Neu-roImage 126, 81-95. doi: 10.1016/j.neuroimage.2015.11.047 . Vigneau, M., Beaucousin, V., Hervé, P.Y., Duffau, H., Crivello, F., Houdé, O., Mazoyer, B., Tzourio-Mazoyer, N., 2006. Meta-analyzing left hemisphere language ar-eas: phonology, semantics, and sentence processing. Neuroimage 30, 1414-1432. doi: 10.1016/j.neuroimage.2005.11.002 . Wang, D., Zhu, Y., Ristaniemi, T., Cong, F., 2018. Extracting multi-mode ERP features using fifth-order nonnegative tensor decomposition. J. Neurosci. Methods 308, 240-247. doi: 10.1016/j.jneumeth.2018.07.020 . Wilkins, K.B., Yao, J., 2020. Coordination of multiple joints increases bilateral connectivity with ipsilateral sensorimotor cortices. Neuroimage 207, 116344. doi: 10.1016/j.neuroimage.2019.116344 . Winkler, A.M., Ridgway, G.R., Webster, M.A., Smith, S.M., Nichols, T.E., 2014. Permutation inference for the general linear model. Neuroimage 92, 381-397. doi: 10.1016/j.neuroimage.2014.01.060 . Yaesoubi, M., Miller, R.L., Calhoun, V.D., 2015. Mutually temporally independent connectivity patterns: A new framework to study the dynamics of brain connectivity at rest with application to explain group difference based on gender. NeuroImage 107, 85-94. doi: 10.1016/j.neuroimage.2014.11.054 . Yamashita, M., Kawato, M., Imamizu, H., 2015. Predicting learning plateau of working memory from whole-brain intrinsic network connectivity patterns. Sci. Rep. 5, 1-8. doi: 10.1038/srep07622 . Yousry, T., 1997. Localization of the motor hand area to a knob on the precentral gyrus. A new landmark. Brain 120, 141-157. doi: 10.1093/brain/120.1.141 . Zhu, Y., Liu, J., Ye, C., Mathiak, K., Astikainen, P., Ristaniemi, T., Cong, F., 2020. Discovering dynamic task-modulated functional networks with specific spectral modes using MEG. NeuroImage 218, 116924. doi: 10.1016/j.neuroimage.2020.116924 .

doi: 10.1016/j.neuroimage.2015.03.071 . Colclough, G.L., Woolrich, M.W., Tewarie, P.K., Brookes, M.J., Quinn, A.J., Smith, S.M., 2016. How reliable are MEG resting-state connectivity metrics? Neuroimage 138, 284-293. doi: 10.1016/j.neuroimage.2016.05.070 . Corbetta, M. , Miezin, F.M. , Dobmeyer, S. , Shulman, G.L. , Petersen, S.E. , 1991. Selective and divided attention during visual discriminations of shape, color, and speed: functional anatomy by positron emission tomography. J. Neurosci. 11, 2383-2402 . Demb, J.B., Desmond, J.E., Wagner, A.D., Vaidya, C.J., Glover, G.H., Gabrieli, J.D., 1995. Semantic encoding and retrieval in the left inferior prefrontal cortex: a functional

(Yamashita et al., 2015)

		fMRI	3-back Working	Static approach	Left fronto-parietal network.
			Memory task
						Supplementary motor and premotor
						areas and the frontal eye field and
						primary motor network.
						Connectivity	between	'lateral
						temporal' network 'middle frontal
						and parietal' network.
	(Heinrichs-	MEG	Sternberg	Time	Frequency	Bilateral occipital cortices
	Graham and		Working	analysis on sensor
	Wilson, 2015)		Memory task	and cortical levels	Left dorsolateral prefrontal cortex
						(DLPFC) and Left superior temporal
						areas.

Acknowledgment

Acknowledgement

This work was financed by the Rennes University and the Institute of Clinical Neuroscience of Rennes (project named EEGCog). The study was also funded by the National Council for Scientific Research (CNRS) in Lebanon. The authors would also like to thank the the Lebanese University, the Lebanese Association for Scientific Research (LASER) and Campus France, Programme Hubert Curien CEDRE (PROJECT No. 42257YA), for supporting this study.

Acknowledgment

The authors would like to thank Martjin van den Heuvel for providing the HCP DTI matrices used to upgrade the COALIA model and Laurent Albera for providing the PSAUD code.

My gratitude extends to each of AZM Saade Association, National Council for Scientific Research in Lebanon (CNRS-L), Lebanese University, University of Rennes 1,

Data availability

Data supporting the findings of this study are available in the link ( https://github.com/judytabbal/dynbrainSS.git ). All HCP data are available on https://www.humanconnectome.org/software/hcp-megpipelines . The datasets 1 and 3 are available upon request.

Code availability

Codes supporting the findings of this study are available in the link ( https://github.com/judytabbal/dynbrainSS.git ). All analysis codes necessary to produce the results here were performed in MATLAB software, using FieldTrip Toolbox http://www.fieldtriptoolbox.org for data segmentation, filtering and source reconstruction steps, EEGLAB toolbox for some SS methods as JADE, InfoMax and SOBI and other MAT-LAB implemented functions as detailed and provided in the previous link. A graphical user interface is made freely available allowing other researchers to test the methods on our (or their) simulated and real data.

Data availability

All the simulated data supporting the findings of this study are available in the Github (https://github.com/judytabbal/dynCOALIA).

Code availability

All codes supporting the results of this paper can be found at (https://github.com/judytabbal/dynCOALIA). All analysis codes were implemented and performed in Code Availability All the above analysis codes are publicly available on GitHub (https://github.com/judytabbal/dynCogPD).

List of Abbreviations AEC

Credit author statement

JT, AK, MH and PB contributed to the design and implementation of the research, to the analysis of the results and to the writing of the manuscript. MK and PB were involved in providing funding sources for the study and supervised the work. All authors provided critical feedback and helped shape the study.

Supplementary materials

Supplementary material associated with this article can be found, in the online version, at doi:10.1016/j.neuroimage.2021.117829 .

Study II: Assessing HD-EEG functional connectivity states using a human brain computational model Judie Tabbal, Aya Kabbara, Maxime Yochum, Mohamad Khalil, Mahmoud Hassan, Pascal Benquet.

Submitted

propagation from the sources located at centroids of the 66 brain atlas regions to the 257 EEG sensors.

The corresponding sources orientation was constrained to be normal to the surface.

Finally, a spatially and temporally uncorrelated white noise was added to the scalp EEG signals to mimic measurement noise [START_REF] Anzolin | Quantifying the Effect of Demixing Approaches on Directed Connectivity Estimated Between Reconstructed EEG Sources[END_REF].

Where 𝑋(𝑡) are the scalp EEG signals and 𝑛(𝑡) is the white uncorrelated noise. ∥∥ 𝐹 correspond to the Frobenius norm and 𝜆 determine noise level added (signal to noise ratio: SNR) varied between 0.7 and 1 (no added noise) with 0.05 step.

Inverse Solutions

Next, temporal dynamics of the cortical sources were reconstructed by solving the inverse problem. It consists of evaluating the parameters of the template source space, including position, orientation, and magnitude of current dipoles. As in the forward model, we constrained the position of cortical sources at the centroid locations of 66 Desikan-Killiany regions, and the orientation to be normal to the cortical surface. Therefore, dipole moment at time t ('𝑆 𝑡 ') can be calculated from sensor EEG time series '𝑋 𝑡 ' as follows:

Where '𝑊' denotes the inverse matrix, often called spatial filters (or weights). Algorithms to estimate '𝑊' can be divided into beamforming and least-squares minimum-norm type estimates. In this study, we are interested in evaluating inverse methods based on the latter type as they are widely used in EEG source connectivity analysis. In particular, we selected two methods: (1) weighted minimum norm estimate (wMNE) and (2) exact low-resolution electromagnetic tomography (eLORETA). These two methods mainly differ in the prior assumptions of the source covariance.

Weighted Minimum Norm Estimate (wMNE)

The weighted minimum norm estimate (wMNE) (Hämäläinen and Ilmoniemi, 1994) searches for a solution that fits measurements with a least square error. It compensates for the bias of the classical minimum norm estimate (MNE) of favoring weak and surface sources. Technically, depth weighting is implemented by introducing a diagonal weighting matrix 𝐵 in order to boost the impact of deep sources:

PCA is a widely used technique that tends to reduce data dimensionality through a variance maximization approach. Hence, k orthogonal variables called 'eigenvectors' are extracted from a set of possibly correlated variables. Here, we applied the Singular Value Decomposition (SVD) algorithm of PCA (Golub and Reinsch, 1970) implemented in Matlab.

Non-Negative Matrix Factorization (NMF)

Nonnegative matrix factorization (NMF) is an unsupervised machine-learning technique (Lee and Seung, 1999) that imposes a 'positivity' constraint on the decomposed factors. In this work, we selected the Alternating Least Squares (ALS) algorithm that has previously shown reliable performance in the fMRI context (Ding et al., 2013) with 100 times replications. For this purpose, the 'nnmf' Matlab function was used [START_REF] Berry | Algorithms and applications for approximate nonnegative matrix factorization[END_REF].

Kmeans Clustering

Kmeans is one of the simplest clustering approaches (Lloyd, 1982). Based on feature similarity, Kmeans assigns each time point to one of the k centroids clusters and the frequency of occurrence of each cluster at each time window is then calculated across all trials/subjects (Allen et al., 2014). In this study, we used L1 distance with 100 times replications and random centroid positions initialization. We adopted the 'kmeans' function incorporated in Matlab [START_REF] Mucha | Cluster dissection and analysis: theory, FORTRAN programs[END_REF].

Performance Analysis

We estimate the similarity between reconstructed and reference sources/networks for each method.

To this end, we defined four metrics: the 'precision', the 'spatial similarity', the 'temporal similarity' and the 'global similarity', to quantify the performance at each step of the pipeline.

Inverse Model evaluation

The 'precision' metric was used to quantify the performance of inverse models at the cortical level.

'Precision' was calculated as the number of 'correct' regions divided by the total number of 'active' regions. In order to ensure equal density across trials and subjects, we used a proportional threshold to keep the top 𝑥% regions with the highest source weights values. As we can see from the simulation scenario (Figure 2), the maximal number of simultaneously activated regions at a time interval is 7 corresponding to approximately 𝑥 = 10% of the total region number (66 desikan). In this work, we chose to vary the threshold from 10% to 15% for 'precision' calculation, and then take the average value across threshold values.

Matlab software using several toolboxes as Fieldtrip (for BEM, eLORETA, and wPLI computation), EEGLAB (for JADE computation), and other hand-written/customized Matlab scripts and functions.

Results

Source localization: wMNE vs. eLORETA

To first visualize results at the cortical-level, we illustrate reconstructed brain sources in Figure 3 at simulated time intervals for each inverse model (wMNE and eLORETA) along with both source space resolution (LowRes vs. HighRes), where 'LowRes' denotes the lead field calculation on 66 desikan regions directly while 'HighRes' refers to lead field calculation on high-resolution cortex followed by projection (average) on 66 desikan regions (see Materials and Methods). Results shown were averaged over trials and subjects.

For inverse methods quantification, we calculated the precision metric averaged over subjects' data (see Materials and Methods) at each simulated time interval between each method and reference sources. In our case, reference sources correspond to the simulated LFP cortical-level data generated by the COALIA model (as illustrated in the first row of Figure 3). Then, precision values were averaged over time intervals to obtain a single quantification value for each method. wMNE showed 80% precision with LowRes and 55% with HighRes while eLORETA showed 66% with LowResres66 and 28% with HighRes (see Figure 3).

We can notice that (1) wMNE outperformed eLORETA for both resolutions and (2) the lead field computed directly on the regions of interest seems to led to higher precision results for both inverse models.

with the highest connectivity values. Edges line width indicates connectivity strength and nodes sphere size reveals brain region strength. Spatial Similarity was calculated between reconstructed and reference brain networks at each time interval with a proportional threshold ranging from 1% to 2% of all possible connectivities. The spatial similarity value is indicated on the right side of each network as the average value over threshold values. PLV networks applied on LFP sources were considered as reference brain networks for wMNE/PLV and eLORETA/PLV. wPLI networks applied on LFP sources were considered as reference brain networks for wMNE/wPLI and eLORETA/wPLI. AEC networks applied on LFP sources were considered as reference brain networks for wMNE/AEC and eLORETA/AEC. A color code was attributed to each method (blue for wMNE/PLV, orange for wMNE/wPLI, yellow for wMNE/AEC, purple for eLORETA/PLV, green for eLORETA/wPLI, and cyan for eLORETA/AEC).

Dimensionality reduction methods evaluation

As we are dealing with six simulated networks, we imposed 'k=6' states for each dimensionality reduction method. As a result, six dynamic brain states (denoted 'Si' in Figure 5), including spatial maps and corresponding temporal signals, were extracted from dynamic functional connectivity networks.

Here, we chose estimated networks from the best combination of previously evaluated inverse model/functional connectivity methods (wMNE/PLV, with 'LowRes') (denoted 'Ci' in Figure 5) to be considered as reference networks for spatial maps evaluation. The six simulated time intervals denoted 'Ti' were considered as reference temporal occupancy for temporal signals evaluation. In this context, spatial similarity and temporal similarity were computed between each connectivity reference and all extracted states. Then, the global similarity was calculated as the average between both spatial and temporal similarities. Following this, each reference connectivity 'Ci' was matched with the most representative state 'Sj' having the highest global similarity value among all extracted states. For Similarly, we evaluated separately the Spatial and Temporal Similarity distribution over 20 subjects of the best matching states, as shown in the boxplots of Figure S4 in Supplementary Materials, followed by the statistical ANOVA test. Our aim is to deepen the performance examination of our tested methods to consider each of the spatial and temporal modes. For the spatial similarity, there were no significant differences (p-value>0.05) between PSAUD (mean=0.7371), NMF (mean=0.7269) and Kmeans (mean=0.7357), having the highest values among all methods, followed by FastICA (mean=0.6955) and PCA (mean=0.6783) with non-significant differences, then JADE (mean=0.6419).

On the other hand, PSAUD (mean=0.7294) exhibited similar temporal performance to each of JADE (mean=0.7355) and FastICA (mean=0.7088) (p-value>0.05). Weaker temporal similarities (with significant differences) were shown for NMF (mean=0.6257), PCA (mean=0.5787), and Kmeans (mean=0.5343).

Moreover, we tested dimensionality reduction methods' dependency on measurement noise (see Figure S5). The maximal global similarity is plotted against different noise levels expressed by λ value for JADE (𝑚𝑒𝑎𝑛 ± 𝑠𝑡𝑑 = 0.68 ± 0.0216), FastICA (𝑚𝑒𝑎𝑛 ± 𝑠𝑡𝑑 = 0.72 ± 0.0060), PSAUD (𝑚𝑒𝑎𝑛 ± 𝑠𝑡𝑑 = 0.73 ± 0.0228), PCA (𝑚𝑒𝑎𝑛 ± 𝑠𝑡𝑑 = 0.64 ± 0.0212), NMF (𝑚𝑒𝑎𝑛 ± 𝑠𝑡𝑑 = 0.68 ± 0.0255), Kmeans (𝑚𝑒𝑎𝑛 ± 𝑠𝑡𝑑 = 0.64 ± 0.0276). Results were slightly influenced by noise with Kmeans being relatively the most dependable method (with the highest std value) and FastICA as the least dependent method.

Supplementary Materials for Study II Assessing HD-EEG functional connectivity states using a human brain computational model

Figure S1. (A) Lobes distribution on the brain cortex used for spatial similarity calculation. The brain cortex is divided into seven main lobes (Brainstorm distribution) including occipital, parietal, temporal, central, frontal, prefrontal and limbic. In total, 28 possible connections types can be established between different brain lobes (i.e. O-O refers to the existing connections between nodes of occipital lobe, while P-T refers to existing connections between nodes in parietal lobe and nodes in temporal lobe). (B) Two exemplar brain networks to

Study III: Dynamic Functional Connectivity during cognitive control in

Parkinson's Disease

Joan Duprez, Judie Tabbal, Mahmoud Hassan, Aya Kabbara, Ahmad Mheich, Paul Sauleau, Julien Modolo, Mark Verin, Pascal Benquet, Jean-François Houvenaghel To be submitted

Dynamic Functional Connectivity during cognitive control in

Parkinson's Disease

Joan Duprez 1 , Judie Tabbal 1,2 , Mahmoud Hassan 3,4 , Aya Kabbara 3,5 , Ahmad Mheich 6 , Paul Sauleau 7,8 , Julien Modolo 1 , Mark Verin 7,9,10 , Pascal Benquet 1 , Jean-François Houvenaghel 7,9

1 Univ Rennes, LTSI -U1099, F-35000 Rennes, France

Abstract

Parkinson's disease (PD) is a progressive neurodegenerative disorder characterized by motor abnormalities, as well as cognitive impairments. Among the various cognitive domains affected by the disease, PD alters cognitive action control which is a critical process in adapting behaviors.

Understanding the dynamics of the cognitive processes underlying cognitive action control alterations in PD patients is a challenging issue. Here, using scalp High Density Electroencephalography (HD-EEG) recorded from 31 participants (10 healthy controls (HC), and 21 PD patients), we explored the cortical dynamic functional connectivity during a conflict task (Simon task). We adopted the EEG source connectivity approach followed by a sliding window technique to reconstruct the dynamic cortical functional networks within both beta and gamma bands. Then, we applied the temporal Independent Component Analysis (tICA) method to derive the relevant functional connectivity states relative to HC and PD groups. Finally, the spatiotemporal alterations between groups were assessed by the means of source-level microstates metrics. Our results reported the spatiotemporal activity of the resultant microstates and showed significant differences in the beta-modulated states between groups. These states mainly involved frontal and temporal connections with variable density among groups. In this study, we aim to explore the dynamic reconfiguration of the functional connectivity states induced by the cognitive control, and examine its relation to the cognitive impairments in Parkinsonians compared to healthy subjects. quantify the time-varying changes in brain networks induced by the conflict task for both HC and PD subjects.

Materials and Methods

Participants

Ten HC (4 males and 6 females), aged between 45-57 years (mean=52.8, std=4.3) and twenty-one patients diagnosed with idiopathic PD (10 males, 11 females), aged between 48-69 years (mean=59.4, std=6.7) participated in this study.

All participants (HC and PD) underwent a comprehensive neuropsychological assessment for global cognition: the Montreal Cognitive Assessment (MoCA) [START_REF] Nasreddine | The Montreal Cognitive Assessment, MoCA: A Brief Screening Tool For Mild Cognitive Impairment[END_REF]. Additional standardized tests were assessed for PD patients, representing several cognitive abilities. These tests included: Symbol Digit Modalities Test (SDMT), Digit span test, Stroop test, the judgment of line orientation test [START_REF] Benton | Visuospatial judgment. A clinical test[END_REF], Boston naming test [START_REF] Graves | Boston naming test short forms: a comparison of previous forms with new item response theory based forms[END_REF]

Experimental Task

In this study, we used a color version of the Simon task [START_REF] Simon | Auditory S-R compatibility: The effect of an irrelevant cue on information processing[END_REF] Participants were placed 80 cm in front of a 22 inches' computer screen. At the beginning of a trial, a central dark fixation cross was presented on a white screen, during a variable period (pseudo-randomly defined from 1750 ms to 2170 ms). Then, a blue or a yellow circle (3.9 cm diameter) was displayed either on the right or the left side during 200 ms (Figure 1). Participants were asked to press the colored button (a blue one in right hand, a yellow one in left hand) corresponding to the color of the circle displayed while ignoring its location. They had to respond within 1000 ms after the stimulus offset. The location side of the colored circle stimulus could match (congruent) or not match (incongruent) the the dFC tensor to extract the dynamic brain network states, including spatial network maps and temporal activity. A null distribution was generated to assess the temporal moments of significant modulation for each of the extracted states (as highlighted in red).

EEG Source Connectivity

Forward Model

Following the equivalent current dipole model, EEG signals 𝑋(𝑡) measured from Q channels can be expressed as linear combinations of P time-varying current dipole sources 𝑆(𝑡) as follows:

Where 𝐺(𝑄 × 𝑃) represents the forward model, often called the lead field matrix, and 𝑁(𝑡) denotes the additive noise.

The lead field matrix 𝐺 is computed from a realistic head model along with the position of electrodes.

In our case, we adopted the Boundary Element Method (BEM) head model fitted to the MRI template (ICBM) [START_REF] Mazziotta | A probabilistic atlas and reference system for the human brain: International Consortium for Brain Mapping (ICBM)[END_REF], downloaded from https://www.mcgill.ca/bic/software/tools-dataanalysis/anatomical-mri/atlases/icbm152lin) using the OpenMEEG toolbox (Gramfort et al., 2010), and selected the Electrical Geodesic Inc (EGI) configuration for the EEG electrodes used.

Inverse Solution: wMNE

The EEG inverse problem consists in estimating the unknown parameters of dipolar source 𝑆(𝑡) at the cortical level (position, orientation, and magnitude), from the measured EEG signals 𝑋(𝑡) at the scalp level. Here, we used the destrieux atlas parcellation (𝑛𝑅𝑂𝐼𝑠 = 148 Regions of Interests) (Destrieux et al., 2010) to locate cortical sources, and constrained their orientation normally to the cortical surface [START_REF] Dale | Improved Localizadon of Cortical Activity by Combining EEG and MEG with MRI Cortical Surface Reconstruction: A Linear Approach[END_REF].

Therefore, the EEG inverse problem was reduced to the estimation of sources magnitude:

To compute the inverse matrix W, we used the weighted minimum norm estimate (wMNE) [START_REF] Lin | Assessing and improving the spatial accuracy in MEG source localization by depth-weighted minimum-norm estimates[END_REF] that compensates for the tendency of the classical minimum norm estimate (MNE) (Hämäläinen and Ilmoniemi, 1994) of favoring weak and surface sources.

Where 𝐵 is the diagonal weighting matrix (inversely proportional to the norm of lead field vectors), 𝐺 is the lead field matrix, 𝛼 is the regularization parameter (computed based on signal to noise ratio: 𝛼 = 1/𝑆𝑁𝑅) and 𝐶 is the noise covariance matrix (calculated from our 700 ms pre-stimulus baseline).

Discussion

In this study, we aim to explore i) how a specific (fast-scale) cognitive processes, as cognitive action control, are associated with dynamically reconfiguring functional connectivity networks, and ii) how this dynamic structures can relate to cognitive impairments in Parkinsonians compared to healthy subjects. Both objectives are challenging and highly required in the field of cognitive and clinical neuroscience. To this end, we used scalp HD-EEG recorded from 31 participants (10 HC subjects, and 21 PD patients) during a popular conflict task, namely the Simon task (visual color version). The cortical functional networks were estimated using the EEG source connectivity method within two frequency bands (beta and gamma). We applied the combination of wMNE/PLV, followed by a sliding window approach to track the dynamics of FC networks. To summarize dFC into a set of relevant connectivity patterns, known as dynamic brain network states (dBNS), a variant of temporal independent component analysis (tICA) was applied, providing a set of group-specific functional states, denoted HC/PD microstates. To examine and quantify the time-varying alterations in these states, we used source-level microstates metrics, followed by a statistical analysis to assess the significance between the two groups (HC vs PD).

Source-level microstates approach -Application to neurological diseases

Recently, there has been much interest in the concept of functional brain states, characterized by a limited number of functional patterns with a temporarily stable activity followed by a fast transition to another state (Baker et al., 2014;[START_REF] Khanna | Microstates in resting-state EEG: current status and future directions[END_REF][START_REF] Michel | EEG microstates as a tool for studying the temporal dynamics of whole-brain neuronal networks: A review[END_REF]O'Neill et al., 2018).

Typically, sensor-level EEG microstate analysis has demonstrated its potential use as a tool to derive the functional brain states [START_REF] Michel | EEG microstates as a tool for studying the temporal dynamics of whole-brain neuronal networks: A review[END_REF], to understand healthy cognition [START_REF] Britz | EEG-microstate dependent emergence of perceptual awareness[END_REF][START_REF] Brodbeck | EEG microstates of wakefulness and NREM sleep[END_REF][START_REF] Seitzman | Cognitive manipulation of brain electric microstates[END_REF], as well as brain diseases [START_REF] Khanna | Microstates in resting-state EEG: current status and future directions[END_REF][START_REF] Musaeus | Microstates as Disease and Progression Markers in Patients With Mild Cognitive Impairment[END_REF][START_REF] Schumacher | Dysfunctional brain dynamics and their origin in Lewy body dementia[END_REF][START_REF] Tait | EEG microstate complexity for aiding early diagnosis of Alzheimer's disease[END_REF]. Here, instead of performing the microstates analysis at the sensor-level (due to anatomical interpretation issues), we projected our HD-EEG data to source space to investigate the dynamics of cortical functional connectivity at a subsecond time scale (Tait and Zhang, 2021).

Although Kmeans clustering has been commonly used in most aforementioned studies, other clustering/decomposition algorithms can also decipher the dynamic brain states. Among these methods, we were particularly interested in the temporal Independent Component Analysis (tICA), as it has proven its potential ability to track the fast temporal variation in brain connectivity (O'Neill et al., 2017;Tabbal et al., 2021;Yaesoubi et al., 2015). This method is robust with our data sample size as we it consists of temporal windows concatenation over trials/subjects to elucidate the group-level

Supplementary Materials for Study III Dynamic Functional Connectivity during cognitive control in

Chapter 4. DISCUSSION AND FUTURE WORK

There is now growing evidence that the human brain is a large-scale complex network of highly interconnected and distributed regions [START_REF] Buckner | The organization of the human cerebellum estimated by intrinsic functional connectivity[END_REF]Bullmore and Sporns, 2009b;[START_REF] Sporns | Organization, development and function of complex brain networks[END_REF][START_REF] Stam | Graph theoretical analysis of complex networks in the brain[END_REF][START_REF] Van Den Heuvel | Functionally linked restingstate networks reflect the underlying structural connectivity architecture of the human brain[END_REF][START_REF] Van Den Heuvel | Small-world and scale-free organization of voxel-based resting-state functional connectivity in the human brain[END_REF]. Its functional organization reconfigures dynamically and flexibly during resting state (Baker et al., 2014;[START_REF] Damaraju | Dynamic functional connectivity analysis reveals transient states of dysconnectivity in schizophrenia[END_REF]de Pasquale et al., 2016) and behavioral tasks (Bola and Sabel, 2015;[START_REF] Bressler | Large-scale brain networks in cognition: emerging methods and principles[END_REF]Hassan et al., 2015;O'Neill et al., 2017) to guarantee a continuous and adaptive information processing between neural assemblies. The investigation of such networks has contributed to pivotal findings related to our understanding of the brain functionality in healthy humans [START_REF] Fox | Spontaneous fluctuations in brain activity observed with functional magnetic resonance imaging[END_REF][START_REF] Greicius | Functional connectivity in the resting brain: A network analysis of the default mode hypothesis[END_REF][START_REF] Van Den Heuvel | Exploring the brain network: A review on resting-state fMRI functional connectivity[END_REF], as well as the neural dysfunctionality due to brain disorders (Du et al., 2016;Fornito et al., 2015a;[START_REF] Pascual-Marqui | Discrete, 3D distributed, linear imaging methods of electric neuronal activity[END_REF].

Therefore, a systematic characterization of the dynamics of healthy and pathological functional brain networks is highly needed in both neuroscience and clinical fields. This requires the use of non-invasive neuroimaging techniques with excellent temporal resolution as EEG/MEG [START_REF] Hutchison | Dynamic functional connectivity: promise, issues, and interpretations[END_REF][START_REF] Kucyi | Dynamic Brain Network Correlates of Spontaneous Fluctuations in Attention[END_REF][START_REF] Vanrullen | The Time Course of Visual Processing: From Early Perception to Decision-Making[END_REF], along with the application of advanced data analysis methods able to track the sub-second temporal fluctuations of the principal brain networks.

The 'EEG/MEG source connectivity' has opened the pathway for the study of the whole-brain functional connectivity with high temporal precision [START_REF] Michel | EEG source imaging[END_REF]. It allows a reliable reconstruction of EEG/MEG sources responsible for the generation of the scalp-level signals.

In this thesis, we explored the performance of the dynamic analysis of the 'EEG/MEG source connectivity' method followed by decomposition techniques to elucidate elementary BNS that fluctuate over recording time during cognitive tasks. In general, qualitative and quantitative evaluations were executed at both group and subject levels using empirical and simulated electrophysiological data respectively. The major contributions of this thesis are summarized here:

 We performed a comprehensive comparative analysis of the performance, functionality, and limitations among the existing dimensionality reduction methods used to summarize the neural activity in a limited number of dominant brain network states. When applied to three empirical MEG data during cognitive tasks, variability in spatial and temporal precisions was revealed.

Among tested decomposition/clustering methods, it can be seen that (high-order) statistical independence describes the best the temporal brain activity of the most representative BNS at both group and subject levels. Our main message is that researchers should be aware to select the Previously, many studies have shown the alterations in the static functional networks associated with brain disorders. For instance, Hassan et al. associated the cognitive decline in PD disease with the loss of frontotemporal connectivity using HD-EEG (Hassan et al., 2017a). Peled et al. exploits sFC to detect the dysconnectivity related to schizophrenic patients from EEG signals [START_REF] Peled | Functional Connectivity and Working Memory in Schizophrenia: An Eeg Study[END_REF]. Moreover, researchers have performed many studies related to epilepsy to characterize the propagation of the epileptic activity using sFC networks estimated from cortical EEG signals, see review [START_REF] Lemieux | Concepts of Connectivity and Human Epileptic Activity[END_REF].

Recent studies have revealed the potential of the dFC to be more sensitive and informative in identifying brain disease-induced effects than sFC. The additional benefit of dynamic network analysis to static network analysis has recently been demonstrated in a patient cohort undergoing respective neurosurgery, where dynamic hub properties were significantly related to pathology impairments, which is not the case when using static approach [START_REF] Carbo | Dynamic hub load predicts cognitive decline after resective neurosurgery[END_REF]. Using the dynamic approach, several works have attempted to explore the differences between groups of controls and patients, as in Alzheimer diseases [START_REF] Kabbara | Reduced integration and improved segregation of functional brain networks in Alzheimer's disease[END_REF]. The effects of anesthesia were also investigated using sliding window functional connectivity [START_REF] Lee | Diversity of functional connectivity patterns is reduced in propofol-induced unconsciousness[END_REF].

Considering the temporal features of FC along with a meta-state approach, many studies have found that summarizing the dFC into a set of time-varying connectivity patterns can provide an accurate characterization of brain disorders. For instance, several findings show the capability of tICA to identify altered spatiotemporal properties in Alzheimer's Disease (AD) (Koelewijn et al., 2017), epilepsy (Koelewijn et al., 2015), and depression (Knyazev et al., 2016;Nugent et al., 2015) using EEG/MEG data during rest and task. PD differences were detected in modular states using the modularity based parcellation algorithm (Kabbara et al., 2019). PCA and Kmeans methods have also unveiled timevarying topological alterations related to Schizophrenia [START_REF] Damaraju | Dynamic functional connectivity analysis reveals transient states of dysconnectivity in schizophrenia[END_REF]Du et al., 2016;Miller et al., 2016), multiple sclerosis (Leonardi et al., 2013), and Parkinson's disease (Kim et al., 2017), yet, in the context of resting-state fMRI.

In study III, we adopted the EEG source connectivity approach to estimate dFC and applied the tICA method using dense-EEG data recorded during the Simon task. Our objective was to identify significant modulations in the derived dBNS between healthy controls (HC) and PD patients. The reported results

showed that this objective can be achieved and may provide new insights for clinical demands in PD disorder.

On the other hand, it is important to mention that the previous methodology could be applied to provide an objective measure of disease severity (Brookes et al., 2016), or assessing the brain networks