Brain network dynamics in the alpha band during a complex postural control task

Objective. To decipher brain network dynamic remodeling from electroencephalography (EEG) during a complex postural control (PC) task combining virtual reality and a moving platform. Approach. EEG (64 electrodes) data from 158 healthy subjects were acquired. The experiment is divided into several phases, and visual and motor stimulation is applied progressively. We combined advanced source-space EEG networks with clustering algorithms to decipher the brain networks states (BNSs) that occurred during the task. Main results. The results show that BNS distribution describes the different phases of the experiment with specific transitions between visual, motor, salience, and default mode networks coherently. We also showed that age is a key factor that affects the dynamic transition of BNSs in a healthy cohort. Significance. This study validates an innovative approach, based on a robust methodology and a consequent cohort, to quantify the brain networks dynamics in the BioVRSea paradigm. This work is an important step toward a quantitative evaluation of brain activities during PC and could lay the foundation for developing brain-based biomarkers of PC-related disorders.


Introduction
The role of postural control (PC) is crucial in daily life. Maintaining upright balance during postural tasks, gait and orientation is demanding and relies on PC (Massion 1994). PC is a complex mechanism based on the interaction of dynamic sensorimotor processes (Horak 2006). It is a central nervous system feedback control system, where neural activity is essential to communicate with muscle activity and physiological as well as visceral autonomic control systems to anticipate and react to gravity, movements, and constantly changing posture (Sherrington 1911, Ivanenko andGurfinkel 2018). Thus, any events provoking a dysfunction of neural activity will disturb PC. Studying PC could help diagnose, or at least differentiate healthy behavior from a pathological one, such as neurodegenerative disease (Buckley et al 2019). This is why the identifying quantitative metrics regarding PC, on the neurophysiological level, is highly interesting. However, even though neural circuitry is prominent in ensuring a steady PC, the underlying neural mechanism is still largely unknown.
Using neuroimaging techniques such as fMRI or functional near-infrared spectroscopy, several studies highlighted the role of brain activity in PC studies (Mirelman et al 2014, Taube et al 2015, to study the impact of brain traumatic injury on balance performance for instance (Caeyenberghs et al 2012). The same approach has been used to study these injuries with electroencephalography (EEG) (Shenoy Handiru et al 2021). EEG proved to be of high interest in PC studies during the control of balance to specific motor stimulation (Solis-Escalante et al 2019) or to decipher adaptation and habituation phenomena (Edmunds et al 2019, Barollo et al 2020). Compared to other imaging techniques, EEG is more flexible: it does not require licensed training, it is inexpensive, but above all, it enables mobile experiment setups and provides a very high temporal accuracy (Cohen 2011), allowing for tracking fast changes of cortical brain activity. For those reasons, EEG is very convenient for PC experiments, where brain activity can be easily recorded while a subject is undergoing different stimulations challenging its balance. Based on EEG analysis, it has been shown that spectral power extracted from the different frequency bands of the brain is an indicator of cerebral cortex involvement during PC tasks (Hülsdünker et al 2015) or to differentiate a healthy behavior from a pathological one (Thompson et al 2005). Other studies used evoked potentials to highlight the correlation between cortical brain activation and balance related tasks (Little andWoollacott 2015, Goel et al 2018). Yet, more advanced techniques using mathematical tools are built to deepen the analysis and identify more accurate metrics.
An innovative methodology based on graph theory has been emerging recently and has shown encouraging results. Brain functional connectivity uses the EEG time series at the scalp level to reconstruct brain sources and compute brain networks. It, therefore, studies the interplay between brain regions as well as the strength of the connections. Its robustness and usefulness have been observed in many studies (Hassan et al 2014, Hassan andWendling 2018) and already showed promising results in terms of PC analysis (Shenoy Handiru et al 2021, Barollo et al 2022). Brain network analysis helps to identify the functional regions involved during a task, thus, the cognitive processes involved. Previous works (Shirer et al 2012) observed specific networks patterns during common tasks (e.g. resting state, visual, auditory tasks). The identification of these patterns and the associated hubs involved during those tasks set the foundation of resting state networks (RSN). Computing brain networks from EEG and comparing their metrics to the known RS, is a concrete application of brain functional connectivity to decipher functional brain response to dynamic tasks. However, for now, the functional connectivity associated with EEG in PC did not fully exploit the high temporal resolution of this measuring technique.
Based on previous studies (Lehmann 1971, Lehmann et al 1998, the analysis of EEG time series showed voltage topographies remain stable for a period of around 80-120 ms. Those stable periods have been mentioned as functional microstates. The study of microstates helped to comprehend brain states and to characterize brain functions in several studies (Khanna et al 2015, Michel andKoenig 2018). More and more tools are available to process and analyze microstates Koenig 2018, Poulsen et al 2018). During dynamic tasks, the review of microstates is a promising tool for decoding brain cortical remodeling. From this aspect, we are interested in studying the dynamicity of brain networks (that are now mostly limited in time) based on the method developed to compute EEG microstates.
The current work will introduce an innovative workflow, combining the two approaches mentioned above. From the EEG time series, sequences of brain connectivity matrices and dynamic brain networks will be computed. Then, we will identify states, by clustering network distributions, mainly due to stable network configurations over time. These resulting states will be called brain network states (BNSs). Studying BNS evolution and characteristics will help to validate this new methodology and open new perspectives on the appreciation of brain function during dynamic tasks. This approach has already been used with a few subjects at rest (Britz et al 2010, Yuan et al 2012 or more recently developed in the study of Parkinson's disease (Duprez et al 2022). Yet, the spatiotemporal dynamics of functional brain networks during complex PC tasks are still unclear.
As a model PC, we have recently developed the 'BioVRSea' (figure 1). This novel measurement setup consists of a moving platform, a virtual reality (VR) simulation mimicking the behavior of a boat in the sea, and several biosignals measurements. BioVRSea paradigm is a dynamic PC task where the stimulation is modulated through time. It already showed promising results in terms of deciphering PC responses, for instance with concussion (Jacob et al 2022), and has been used as a first step toward the identification of markers of motion sickness (Recenti et al 2021). Here, our objective is to decipher the spatiotemporal dynamics of brain networks during the BioVRSea paradigm. We believe that this approach will allow us to observe the dynamic remodeling of the brain during a complex PC task by watching the dynamic network evolution through the experiment. To do so, we use advanced EEG analysis consisting of constructing the time-varying functional brain networks at the source level. The k-means clustering algorithm was then used to segment these time-varying connectivity matrices into a set of BNS. Applied to a relatively large sample size (N = 158) of healthy participants, our results showed that this data-driven segmentation matches very well with different transitions during the task.
The work presented in this study would be the first to show spatiotemporal brain network dynamics in 158 individuals.

Participants
Participants were recruited through snowball sampling. Written information about the study was provided, and all the participants had to sign a written informed consent before starting the acquisition procedure. At the time of the study, 158 individuals (93 women, 65 men, age range = 19-72 years old, mean age = 32.9 years old, std age = 13.2 years old) completed the experiment.

Questionnaire
Before the acquisition, the participants were asked to answer a Motion Sickness Susceptibility Questionnaire (MSSQ) based on the adult questionnaire proposed by Golding (Golding 2006). Subsequently, the participants had to fill out a motion sickness symptoms questionnaire before and after the experiment. This questionnaire evaluates 10 motion sickness symptoms on a score from 0 to 3: General discomfort, Fatigue, Headache, Eye strain, Difficulty focusing or concentrating, Increased salivation, Sweating, Nausea, Blurred vision, Dizziness, or vertigo. From that information, two scores were calculated, the MSSQ score and the symptoms score. The MSSQ score is calculated by using the same formula explicated in Golding. This MSSQ is composed of 9 items (scores from 0 to 3), with a maximum score of 27. The symptoms questionnaire score is calculated by subtracting the average of the symptoms before the experiment from the average of the symptoms Figure 2. Data processing workflow. From the EEG time-series, ROIs are reconstructed. The functional connectivity between those ROIs is computed. Then, BNS are segmented at the group level using clustering algorithm. Graph theory metrics are extracted from the network to quantify and analyze their functions. Then, the segmentation is studied at the subject level to investigate the prevalence of each BNS.
after the experiment, to obtain a final score between 0 and 3.

Acquisition
The acquisition procedure is the same one that has been described in Jacob et al (2022) and Aubonnet et al (2022). Before the acquisition started, the subject was equipped with a wet 64-electrode EEG cap (Sampling frequency 4096 Hz, ANTNeuro, Hengelo, The Netherlands). The participant was then asked to step onto the platform. The subjects had to stand quietly on the platform with their hands by their side observing a mountain view for 1 min 45 s (the last 60 s of this step were used as a baseline). Then, the scene in the VR goggles would change, with a white screen preparing for the sea simulation (15 s). The participants were instructed to stand quietly with their hands by their sides for the first 40 s of the sea simulation. There was no platform movement in this part of the experiment, called task PRE. After 40 s of quiet standing watching the sea simulation, the participant was instructed to hold onto the bars in front of them. The platform then began a synchronized movement with the sea scene in the VR goggles, with 25%, 50%, and 75% of maximal wave amplitude. In this central part, each segment lasted 40 s and the participant held the bars of the platform while continuing to observe the sea simulation. Finally, the platform stopped moving and the participant was asked to remove their hands from the bars and attempt to stand quietly with their hands by their side for the final 40 s of the experiment. The participant still observed the sea scene for the final 40 s. This is called the POST phase of the experiment; it is performed identically to the PRE phase but after the participant has performed movements in the  Figure 2 sums up the data processing workflow.

Data preprocessing
The EEG was recorded using a 64-electrode channel system. Data pre-processing and analysis were performed with Brainstorm (Tadel et al 2011) and Mat-lab2021b (MathWorks, Inc. Natick, 158 Massachusetts, USA), using the Automagic toolbox (Pedroni et al 2019). The preprocessing pipeline was applied per subject, on their whole recording. The data were resampled to 1024 Hz. The data were notch filtered at 50 Hz. A high pass and low pass filter were set respectively to 1 and 45 Hz. Automagic was used to automatically pre-process every dataset, with a manual inspection at the end. The ICA MARA algorithm was used (Winkler et al 2011(Winkler et al , 2014. The residual bad channels were eliminated, identifying the ones with a standard deviation exceeding 20 µV. Finally, bad electrodes were interpolated. Each subject requiring interpolation of more than 15% of the total amount of electrodes was rejected, leading to 158 complete EEG recordings lasting around 318 s. Table 1 sums up the information of the artifactual IC components and bad channels.

Data processing-Dynamic network computation
In the two following sections, several variables are expressed to describe the matrices information. Table 2 sums up the listed variables and their description.
Brain networks were constructed from the pre-processed EEG signals, using an approach called 'dense-EEG source connectivity' , which was developed and used in previous studies (Hassan et al 2014, Kabbara et al 2017, Hassan and Wendling 2018, Duprez et al 2022. This method aims to solve the inverse problem and reconstruct brain sources from EEG time-series. The 45 first seconds of the acquisition have been used to constitute the noise covariance matrix, whereas the sources will be reconstructed from the 273 remaining seconds. Then, the weighted minimum-norm estimation (wMNE) method was used to obtain the source distribution (Hassan et al 2014). Here, we used the Boundary Element Method head model fitted to the ICBM MRI template, composed of three layers (scalp, outer skull and inner skull), using the OpenMEEG plugin from the Brainstorm toolbox. The model's morphology was based on a standard head model template provided by BrainStorm (Tadel et al 2011), and has shown pertinent results in previous studies (Edmunds et al 2019, Kabbara et al 2021, Barollo et al 2022. Based on the Desikan-Killiany atlas (Desikan et al 2006), 56 region of interest (ROI) were estimated. The list of the ROI is provided in the supplementary material. Finally, the connectivity matrices between these 56 ROIs were computed, for the alpha band (8-13 Hz), using the phase-locking value (PLV), a sliding window approach, with a window size of six cycles (Lachaux et al 2000), resulting in a window of 571 ms. It leads to N = 478 matrices per subject, (corresponding to the 273 s of recordings), of dimensions 56 × 56, where each element Ai,j of a matrix represents the weight of the connection from the node i to the node j. This element is called edge.

Data processing-BNS identification
From this final dataset, BNSs segmentation was performed using clustering algorithms to identify stable periods over time. The analysis was done based on some functions developed in the EEGLab Microstate toolbox (Poulsen et al 2018). This toolbox is designed for EEG voltage topography, but we adapted the pipeline to use it with functional brain network as an input. Since the n × n connectivity matrix (n = 56) is symmetric, we vectorized it in a (n * (n − 1)/2) = 1540 elements vector, leading a total of 158 matrices of size 478 × 1540. Those matrices have been averaged across the 158 subjects, to obtain a final averaged matrix of size 478 × 1540, that will be used as the input for the state segmentation. Then, the modified k-means algorithm (Pascual-Marqui et al 1995) was used as a clustering algorithm, and the BNS were sorted according to global explained variance (GEV) (Murray et al 2008, Poulsen et al 2018. We tested the segmentation from 2 to 20 BNS to identify the right number of BNS. The number of random initializations was set at 500, and the max number of iterations was set at 1000, with a default convergence threshold of 10 −6 . To ensure an optimal segmentation, the process has been rerun several times (in our case four times). It has been compared qualitatively to select the most suited amongst the ones presenting good scores in validation criteria, such as GEV (Murray et al 2008), cross-validation (Pascual-Marqui et al 1995), as detailed in Poulsen et al (2018). After the segmentation, a smoothing was performed to remove the small peaks or unstable states. All states being present in less than 20 consecutive matrices were changed to the next most likely BNS, based on global map dissimilarity (GMD) (Murray et al 2008).

Data processing-BNSs and resting-states network correspondence
Each segmented state was then computed from a vector to a 56 × 56 symmetric matrix. This represents the functional network associated with the BNS. Based on the studies published by Kabbara et al (2017), and Shirer et al (2012), we associated each of the 56 nodes with seven RSN: DMN, DAN, SAN, MOT, AUD, VIS, and Other. The affiliation of each node and RSN is summed up in the supplementary material. Then, using Brain Connectivity Toolbox (Rubinov and Sporns 2010), the strength (sum of weights of each link connected to the node) of each node was computed. We then computed the strength of each RSN by summing the values of each node associated with this RSN and dividing it by the number of nodes belonging to the RSN. To identify the dominant RSN for each BNS, we elicited those with a value superior to the sum of the average and standard deviation of the seven RSN strengths. For visualization purposes, only the nodes with a strength >0.1 were displayed, using BrainNetViewer (Xia et al 2013). Each ROI defined from the Desikan-Killiany atlas is affiliated with the corresponding RSN.

Data processing-Quantification
After the BNS segmentation was performed on the averaged matrix, we returned to the subject level by performing a backfitting. Each sample is associated with the most similar BNS, based on GMD (Murray et al 2008, Poulsen et al 2018. As a result, we have a BNS distribution for each subject. To analyze it, the percentage of BNS repartition per subject is computed for each of the 478 matrices. Moreover, two metrics are computed: first, the transition probabilities matrix, revealing the transition probabilities from one state to another, and the occurrence, which indicates the average number of times per second a BNS is dominant (Poulsen et al 2018).

Data processing-Relation with questionnaire or age
The relationship between metrics (global transition matrix, global occurrence) has been investigated with age groups and questionnaire results groups, dividing them with the following criteria: For the age, we chose the 20% top and 20% lowest to observe the age difference impact and have enough subjects in each group to draw a relevant statistical analysis. For the MSSQ and MSS, we chose 10% of the top and lowest scores, to remain coherent in the groups, and not include people that were not feeling sick or prone. Even though the number of individuals is not high, it is enough to perform preliminary statistical analysis, and therefore observe a trend in the results.

Statistics
For the relationship between age or MS related questionnaire, a Wilcoxon rank sum test (α = 0.05) has been performed. If a parameter was observed as statistically significant, it was then corrected for multiple comparisons using a false discovery rate (Benjamini and Hochberg 1995).

BNS segmentation
The functional brain networks were computed over the entire experiment, in the alpha band-due to its significance in previous connectivity and PC studies (Aubonnet et al 2022, Barollo et al 2022-using the EEG source connectivity method (Hassan et al  2014). Briefly, we estimated the dynamics of brain sources, and then we computed the statistical coupling between brain regional time series (N = 56). This was done using a sliding-window approach (window size of six cycles according to Lachaux et al (Lachaux et al 2000). This yields 478 connectivity matrices for the entire task. Then, using the k-means algorithm, we wanted to cluster those matrices into BNSs, a set of networks that dominate over time. The BNS segmentation was performed on all the subjects' averaged dynamic networks. As detailed in the Material and Methods section, the BNS segmentation has been repeated several times. A qualitative inspection has been done in each of them, to result in a final number of clusters k = 5.
The segmentation leads to five BNS during this PC paradigm (figure 3). The correspondence of each BNS to a functional network (according to previous studies (Shirer et al 2012, Kabbara et al 2017 is illustrated in figure 4. The baseline is composed of BNS 4 (default mode network-DMN-), BNS 1 (salience network -SAN-), and BNS 5 (SAN). BNS 3 (auditory and motor networks) appears slightly before the transition phase, until slightly after the PRE phase.
Then BNS 1 appears again until the end of the PRE phase. BNS 3 is present for the whole 25 and 50% tasks, and switches to BNS 1 right before the 75% task. A bit before the POST phase, BNS 1 switches to BNS 2 (Visual network -VIS-), until a third of the POST phase, where it switches to BNS 5.

Quantification
Following this global average segmentation, we investigated the presence of each state at the individual subject's level. This process is called backfitting and consists of associating with each EEG sample the state with which it is the most similar. The backfitting allows us to compute several metrics to quantify, for instance, the transition probability between one state to another and the occurrence of each state per phase of the experiment. Figure 5(A) shows the BNS repartition (in %) after the backfitting. The mentioned metrics completes it, e.g. in figure 5(B), transition matrix probability from one state to another (where the element Ai,j represents the probability of switching from state i to state j), and in figure 5(C), the global occurrence of each BNS which indicates the average number of times per second a BNS is dominant.
The BNS 1 (salience network) is the most present state, being there in around 35% of the individuals during the whole experiment. The second most present states are BNS 4 (default mode network) and BNS 2 with around 15%-20% presence, then BNS 3 and BNS 5, around 10%. The highest transitions ( figure 5(B)) are from all states to BNS1, the state with the highest occurrence in every phase, followed by BNS 4.

Association between network states and questionnaire and age
Here, we investigate the possible correlation between the brain network dynamics metrics (transition matrix and occurrence) as described above and some of the participants' information such as age and motion    sickness assessments (as stated in section 2.9). No significant differences were observed between subjects with high MSSQ scores compared to subjects with low MSSQ scores. The same outcome is observed regarding the MSS scores. However, significant differences were observed between the two age groups, regarding the transition probability from BNS 5 to BNS 1 (figure 6(A)). Figure 6(B) shows a boxplot displaying the significant transition probability value for each group.
It can be seen that most of the young cohort has no transition from state 5 to 1, whereas the transition is more present for the old group.

Discussion
The results from this study highlight a dynamic remodeling of the brain during a complex PC task. A novel methodology has been developed using microstate segmentation at the brain network level.

Group level
A BNS segmentation has been performed from the average dynamic functional brain networks. The first observation is the fact the BNS distribution is representative of the experimental phases. Indeed, for each phase change, a state transition happens slightly before or after. The only exception is for the transition between 25% and 50% phases, which can be easily explained by the fact that there is no difference in stimulation between the tasks, only a higher amplitude of movement. Going deeper into the analysis, we break down each phase and observe the repartition of the BNS. It can be noted that only two states appear only one time in the distribution: BNS 4, at the beginning of the baseline, and BNS 2, between 75% and POST. From the results in figure 4, BNS 4 primary function is linked with DMN. The activation of this network at the baseline is coherent with the paradigm phase: the baseline starts when the subject is already standing on the platform with the same view for around a minute. DMN is a state that is activated when the subject is in the resting-state mode and the focus is on the internal mental state process (Andrews-Hanna et al 2010, Raichle 2015, or mind wandering (Mason et al 2007), and is deactivated when a cognitive task starts.
BNS 2 primary function is linked with VIS. This state occurs at the end of 75% and the beginning of POST, so when the platform is moving at the maximal amplitude and then stops moving, while the subject is always seeing the sea simulation inside the VR headset. The subject must maintain its balance during this phase, but its PC is not challenged by any motion stimulus, only visual. The presence of a visual network state during this period testifies to the cognitive involvement into this task, as its prominence during visual stimulations experiments has been highlighted in previous works (Haynes andRees 2005, Kamitani andTong 2005).
BNS5 appears two times, at the end of baseline and the end of POST. Its primary function is linked with SAN. BNS 1 appears in the middle of the baseline, in the whole PRE phase, and in the whole 75% phase. Its primary function is linked with SAN. Studies showed that SAN interacts with DMN and other executive networks (Menon 2015). This explains the transition from state 4, DMN, to state 1 from state 1 to state 5 in the baseline, as during baseline, the brain is at rest and preparing for what comes next. Moreover, SAN is also more presents when the brain needs to be aware of what it should focus on to move, quickly between information gathering, decision planning, and execution (Uddin 2015). State 1 is present shortly after the PRE starts, which is the first actual PC disturbance: the subject needs to maintain equilibrium even though the platform is not moving, due to the visual stimulation. It needs to adapt and understand where the mismatch between what is seen and what is felt comes from. This also occurs with BNS 5 in the POST phase. Then, BNS 1 is present for the whole 75%, while the platform is still in motion, after already moving for more than a minute. We can hypothesize that the motor cortex is less activated due to the habituation of the phenomenon and that the SAN is dominant in that case to anticipate the upcoming task. To differentiate the 2 BNS, we can look at their second most dominant RSN. For BNS 1, the second RSN is the DAN, which is involved during passive viewing tasks (Bray et al 2015), and in general a focus on a particular task, such as selective attention (Szczepanski et al 2013). This has been experimented during the 75% task, where the subject responds to the platform movements and focuses on a specific point of reference in the VR environment. For BNS 5, the second BNS is Other, which corresponds to the 7 other RSN identified by Shirer et al (Shirer et al 2012), but is not illustrated in this study. However, the strength of the other RSNs is in the same range. This BNS does not have a second dominant RSN. This would explain that it follows BNS1 in the baseline, as they are both salience and prepared to switch to another state. Moreover, this BNS is also at the end of the experiment, in the second half of the POST phase, showing that the subject is focused on information gathering, and decision making to maintain the balance, but also knowing that the experiment is coming to its term. Then the salient network helps with the anticipation of the end of the task.
BNS 3 is present during the whole Transition phase, and the whole 25% and 50% phases. Its primary function is linked to AUD and MOT. During the Transition phase, the screen turns white and the sound of the waves and environment appears, while it is completely silent during the baseline. Since there is no motion or visual stimulation, the presence of a dominant auditory network is coherent, as previous studies identified this network to be more activated during audio stimulation (Laird et al 2011). Secondly, state 3 is also present during the whole 25% and 50% tasks, so when the platform activates and starts moving at different intensities. During this period, we ask the subject to put his hands on the bar. Those phases of the experiment are based on motion stimulation, and the dominance of the motor network is perfectly illustrating this point, as the motor network is primarily solicited during tasks requiring movements (Seidler et al 2015).
Overall, a change in connectivity or BNS is observed when a change in the perturbation is observed. This is particularly striking for the changes from the Transition phase to the PRE phase, which is the start of a visual-only perturbation onset, and then from the PRE phase to the 25% phase, where the perturbation becomes mostly motor. This is in line with results observed in the literature, where visual perturbation showed changes in alpha connectivity during PC tasks, mostly during baseline standing (Peterson and

Subject level
The global BNS segmentation based on the alpha band network shows a clear dynamic brain remodeling during the experimental protocol. Thus, it is interesting to observe the results at a subject-level, and determine if a similar trend can be observed. The backfitting results reveal a clear dominance of BNS 1 during the whole experiment. This BNS is associated with a salience network, which is known to interact between the cognition of DMN and the executive network (Menon 2015). Since our experiment is dynamically composed of different stimulations (visual, motor), the salience network plays a key role in catching the attention of where the brain should focus, and facilitating the transition between information gathering, decision planning, and internal focus (Uddin 2015). This is supported by the average transition matrix. It can indeed be seen that the highest transitions occur from state 1 towards every state, and from states 2 and 4 to state 1, confirming the role of BNS 1 in regulating and orientating brain cortical activity depending on the stimulation.
According to the occurrence metrics, after BNS 1, BNS 2 and 4 have a similar presence during the whole experiment. As the visual stimulation is constant during the experiment, through the VR headset, it is then expected to see an average occurrence of state 2. Regarding state 4, since the DMN testifies to internal focus and cognition (Andrews-Hanna et al 2010, Raichle 2015, it can be hypothesized that a good part of the participants was not challenged a lot during the simulation, and therefore did not need to be involved a big part of their cognition during the experiment.
The analysis of the whole cohort gives promising results. The BNS segmentation reveals an explicit remodeling of the brain during this dynamic PC task. The prevalence of RSN depending on the phase shows the pertinence of this approach.

Subgroups analysis
This study also aimed to study the potential differences in the cognitive behavior between two groups based on different criteria: motion sickness susceptibility and age.
On the age aspect, significance is observed in the transition probabilities. The transition from BNS 5 to BNS 1 is significantly higher in the 30 oldest people than in the 30 youngest. The two states have a salience network as the most dominant RSN, so they principally have the same primary function. However, the second most dominant RSN is DAN for BNS 1 and Other for BNS 5. Thus, since DAN is activated when there is a particular focus on the task, we can speculate that the oldest people, while their salience network is activated, have increased activation of the attentional network, as they need more attention to the task and might suffer more from the experiment, due to the age-related decline in PC, and a need of increasing attention to maintain balance (Maki and McIlroy 1996).

Methodological considerations 4.4.1. Network reconstruction methods
To reconstruct brain connectivity matrices, the combination of wMNE and PLV has been used, as shown in the literature to be very efficient in precisely identifying cortical brain networks from scalp EEG during cognitive activity (Hassan and Wendling 2018). We acknowledge that PLV is a measure that is more sensitive to leakage. However, although this is not a perfect form of control, we used source reconstruction before estimating PLV. Also, although zero-lag methods (such as WPLI) help in this matter, they can also remove true zero-lag connectivity which can be present and even hold an important place in the data (Roelfsema et al 1997, Palva and Palva 2012, Viriyopase et al 2012, Gollo et al 2014, Finger et al 2016. We believe there is no perfect solution for this specific matter. Still, we chose to keep such potential zero-lag connections, especially because recent work from our lab showed good accuracy in source connectivity using the wMNE-PLV combination in a simulation study (Allouch et al 2022). However, we believe that it does not affect the approach presented in this paper, i.e. the combination of functional connectivity and microstate segmentation.

Clustering algorithm
The modified k-means algorithm used to cluster the BNS is stochastic, meaning that it will give different results every time it is run. However, it should converge to an optimal solution, thanks to the settings used. In our study, we are setting the algorithm to 500 restarts and 1000 iterations, with a converging threshold of 10-6, which is sufficient to obtain a reliable segmentation.
Regarding the use of GEV and GMD for BNS sorting and assignation, it is the gold standard measure for EEG voltage time series microstates segmentation (Murray et al 2008). It relies on global field power (GFP), representing the standard deviation across all electrodes. Since our approach is based on dynamic connectivity matrices computation, we believe that the use of GMD is pertinent, as our input data can be assimilated to EEG times series (in our case, 478 matrices of 1540 elements each). Therefore, GFP can be calculated across all the elements, and GMD can be computed for our brain state segmentation.
Concerning the selection of the optimal number of clusters, as stated by Poulsen et al (2018), there is no unique valid answer, because many clusters can describe the data. Even though different criteria exist, a final qualitative decision is still required to define the number of clusters, also based on the physiological coherence of the data. Our work is empirical, and the resulting states segmentation obtained from the modified k-means algorithm clustered our data into five BNSs. This segmentation adequately describes the different phases of the experimental protocol, with a suitable functional network definition.

Variability of the results
In this work, we selected the modified k-means clustering as it is the gold-standard in the segmentation of EEG-based states. Nevertheless, as for other methodological choices, there is a place for high analytical variability. Here we did not aim to analyze this analytical variability related to the clustering algorithms. Some other standard algorithms such as the global kmeans (Likas et al 2003) are not adapted to our data due to the ratio between the higher number of channels compared to the number of samples, leading to no satisfactory clusters.
The use of a stochastic approach used in the paper is a limitation to the method; however, the modified k-means approach has been proven to be an efficient approach, and especially adapted to EEG microstate segmentation.
We are confident that a stochastic approach is still valid to cluster connectivity data, as several studies used k-means algorithm to identify functional connectivity staters based on fMRI data (Allen et al 2014, Calhoun et al 2014. Indeed, the use of a very high number of restarts help to diminish the variability of the results. On top of that, temporal smoothing reduces the variability by removing the 'artefact' clusters with little appearance, and strengthen the resulting segmentation. This is why we believe it is coherent with the global workflow of our approach, which is the application of dynamic functional connectivity to the clustering method derived from the EEG microstates. We are aware about the high analytical variability that can occur by changing one method in the whole pipeline such as changing the clustering algorithm as we have done. We expect different results using different clustering algorithms. Analyzing the analytical variability related to the clustering algorithms can be indeed of interest but it is out of the scope of our paper. This is a key issue indeed in the EEG community (and the neuroimaging in general) and some of the co-authors are investigating the analytical variability in the EEG analysis in their recent studies such as the pre-processing  and the EEG network analysis (Allouch et al 2022). Therefore, we believe that there is, at the moment, no perfectly one validated approach to cluster BNS from EEG signals. In our specific case, this risk was reduced as we have a priori (some kind of ground-truth) about the different phases of the experiment that can may guide us to be more confident about the results.

BNS and RSN assignment
In our case, we reconstructed brain ROIs based on the 68 ROIs of the Desikan-Killiany atlas (Desikan et al 2006). However, we aggregated some regions to get 56 ROIs (aiming at a lower number of ROIs than the number of channels). Regarding their association with RSN, we decided to use the non-overlapping approach which is currently the most used in the network neuroscience community. We acknowledge that other atlases with different modalities (size, anatomy, function, overlap) exist and could be of interest in the frame of RSN attribution interpretation. Concerning the present work, the aim of studying RSN was to understand the dominant function of the resulting BNS.
Regarding the transition between BNS, thanks to the excellent time resolution of the EEG, the dynamics of the brain network can be analyzed at different levels going from sub-second to minutes and hours. The minimum time window here is controlled by using PLV in the alpha band (considering the condition of having a specific number of cycles, in our case, six cycles). Thus, no transition dynamic is expected in a time window lower than 571 ms. In some conditions such as the stimuli-related paradigm, we have shown during a picture naming task that brain network dynamics can be reconfigured quickly (Hassan et al 2015). This was possible because PLV was computed at an ms time scale (in an inter-trial synchronization manner).
Here, our aim was to look at the network reconfiguration and transitions over the entire task and we tried to match these reconfigurations with the reference. This was done by several recent papers using fMRI during a movie (Meer et al 2020) for instance.

Focus on the alpha band
This work is focused on the alpha band, due to its significance in previous studies associated with PC (Barollo et al 2022), and more specifically BioVRSea (Aubonnet et al 2022). Moreover, results from literature highlighted the prevalence of RSN in this frequency band (Kabbara et al 2021). Thus, we aimed to investigate more in depth this band, insisting on the connectivity aspect. Future works will observe how the brain network segmentation and the RSN assignment result in the other frequency bands and on the global frequency spectrum.

Limitations
No significance has been found in global metrics related to motion sickness parameters. This aspect should be pursued to adapt the group definition, as well as the characterization of motion sickness in people.
Lastly, this work has been done on healthy people only. Future studies should compare this work with groups suffering from pathologies affecting their neural activity or PC. This distribution can also be compared with new data from a healthy cohort to confirm that the results obtained in this study are consistent.

Conclusion
This study aimed to decipher brain network dynamics during a complex PC task. 158 individuals underwent the BioVRSea experiment. We demonstrated that, on the alpha band, the BNS segmentation appropriately illustrated the evolution of the different phases of the experiment. Each BNS was studied to determine its dominant RSNs, giving a better insight into functional connectivity and how it translates brain behavior in PC. The BNS distribution has then been monitored on the subject level, showing the prevalence of one state, which comports a dominant salient network.
Finally, a proof of concept has been performed to study differences between groups, based on age, and motion sickness susceptibility. The results showed statistical differences in some metrics and transition probabilities for the age difference.
This work validates an innovative approach, based on a robust methodology and a consequent cohort, to quantify the brain networks dynamics in the BioVRSea paradigm. Further studies will confirm those results by comparing new data to this distribution, and observing the distribution of people presenting pathologies. This is the first step toward defining a reference of brain network behavior in dynamic PC.

Data availability statement
The data cannot be made publicly available upon publication because the cost of preparing, depositing and hosting the data would be prohibitive within the terms of this research project. The data that support the findings of this study are available upon reasonable request from the authors.

Acknowledgments
RA contributed to the data acquisition, processing, analysis and interpretation, and original manuscript writing. MH contributed to the data processing and analysis design, data interpretation, study supervision, and manuscript revision. AM contributed to the data analysis and interpretation and manuscript revision. GDL contributed to the data interpretation, and manuscript revision. HP contributed to the experimental paradigm design, and manuscript revision. PG contributed to the experimental paradigm design, study design, study supervision and manuscript revision.

Ethics statement
The National Bioethics Committee of Iceland (no: VSN-20-101) approved the study. The study followed all relevant guidelines for the study. The participants provided written informed consent for the study.