Reliability of Static and Dynamic Network Metrics in the Resting-State: A MEG-Beamformed Connectivity Analysis

The resting activity of the brain can be described by so-called intrinsic connectivity networks (ICNs), which consist of spatially and temporally distributed, but functionally connected, nodes. The coordinated activity of the resting state can be explored via magnetoencephalography (MEG) by studying frequency-dependent functional brain networks at the source level. Although many algorithms for the analysis of brain connectivity have been proposed, the reliability of network metrics derived from both static and dynamic functional connectivity is still unknown. This is a particular problem for studies of associations between ICN metrics and personality variables or other traits, and for studies of differences between patient and control groups, which both depend critically on the reliability of the metrics used. A detailed investigation of the reliability of metrics derived from resting-state MEG repeat scans is therefore a prerequisite for the development of connectomic biomarkers. Here, we first estimated both static (SFC) and dynamic functional connectivity (DFC) after beamforming source reconstruction using the imaginary part of the phase locking index (iPLV) and the correlation of the amplitude envelope (CorEnv). Using our approach, functional network microstates (FCμstates) were derived from the DFC and chronnectomics were computed from the evolution of FCμstates across experimental time. In both temporal scales, the reliability of network metrics (SFC), the FCμstates and the related chronnectomics were evaluated for every frequency band. Chronnectomic statistics and FCμstates were generally more reliable than node-wise static network metrics. CorEnv-based network metrics were more reproducible at the static approach. The reliability of chronnectomics have been evaluated also in a second dataset. This study encourages the analysis of MEG resting-state via DFC.

The resting activity of the brain can be described by so-called intrinsic connectivity networks (ICNs), which consist of spatially and temporally distributed, but functionally connected, nodes. The coordinated activity of the resting state can be explored via magnetoencephalography (MEG) by studying frequency-dependent functional brain networks at the source level. Although many algorithms for the analysis of brain connectivity have been proposed, the reliability of network metrics derived from both static and dynamic functional connectivity is still unknown. This is a particular problem for studies of associations between ICN metrics and personality variables or other traits, and for studies of differences between patient and control groups, which both depend critically on the reliability of the metrics used. A detailed investigation of the reliability of metrics derived from resting-state MEG repeat scans is therefore a prerequisite for the development of connectomic biomarkers. Here, we first estimated both static (SFC) and dynamic functional connectivity (DFC) after beamforming source reconstruction using the imaginary part of the phase locking index (iPLV) and the correlation of the amplitude envelope (CorEnv). Using our approach, functional network microstates (FCµstates) were derived from the DFC and chronnectomics were computed from the evolution of FCµstates across experimental time. In both temporal scales, the reliability of network metrics (SFC), the FCµstates and the related chronnectomics were evaluated for every frequency band. Chronnectomic statistics and FCµstates were generally more reliable than node-wise static network metrics. CorEnv-based network metrics were more reproducible at the static approach. The reliability of chronnectomics have been evaluated also in a second dataset. This study encourages the analysis of MEG resting-state via DFC.

INTRODUCTION
The coordination of spontaneous activity can be characterized with functional connectivity (FC), which refers to statistical dependencies between the activity of distinct brain areas (Pereda et al., 2005) and has been linked to the efficiency of an individual's brain functioning (Baldassarre et al., 2012;Yamashita et al., 2015).
A functional connectivity graph (FCG) can be constructed by estimating the statistical dependencies between the brain activity of all the areas in a pair-wise fashion. An FCG represents statistical or causal relationships measured as cross-correlations, coherence, or information flow (Dimitriadis et al., 2009(Dimitriadis et al., , 2015d. Neuroscientists first examined resting-state FC with functional magnetic resonance imaging (fMRI) by correlating blood oxygenation level-dependent (BOLD) signals (Biswal et al., 1995;van den Heuvel et al, 2009;Biswal, 2011Biswal, , 2012. After 20 years of using fMRI as a dominant neuroimaging tool, the community has succeeded in mapping brain areas to specific brain functions, creating an anatomical-functional atlas (Bandettini, 2012). Although fMRI is of high interest and a key modality to explore human brain function, ultra-slow activity described via BOLD signals is only an indirect measure of brain activity (Logothetis, 2008).
In the last few years, greater attention has been given to explore FC via electro-magneto-encephalography. Even though the spatial resolution of magnetoencephalography (MEG) is lower when compared to fMRI, MEG can capture the multiplexity of human brain activity by providing insight into the spectro-temporo-spatial dynamics of human brain activity. MEG-based FC provides us with a direct measure of neuromagnetic activity with a high temporal resolution (Deco et al., 2011).
Resting-state networks (RSNs) have been successfully extracted with MEG over the past few years using source-space FC (de Pasquale et al., 2010;Brookes et al., 2011a,b;Hipp et al., 2012;Luckhoo et al., 2012;Hall et al., 2013;Wens et al., 2015). Moreover, resting-state MEG FC has been proven to detect abnormal brain functioning in a variety of diseases, including Alzheimer's disease (López et al., 2014(López et al., , 2017Engels et al., 2015), multiple sclerosis (Tewarie et al., 2015), in schizophrenia (Bowyer et al., 2015), in dyslexia (Dimitriadis et al., 2013b(Dimitriadis et al., , 2015c, in mild cognitive impairment (Dimitriadis et al., 2015b) and in mild traumatic brain injury (Dimitriadis et al., 2015c;Dunkley et al., 2015;Antonakakis et al., 2016Antonakakis et al., , 2017. Several studies have thus captured alterations of MEG parameters in the resting state in order to estimate FC in disease groups compared to controls. However, FC estimates at restingstate could be affected by subject's cognitive, emotional state and other scanning-related systematic differences. For that reason, it is unclear up to which level FC estimates are repeatable for an individual. Moreover, in large studies of hundreds of participants, there is a significant cost, both in financial resources and time, to scan all the subjects two or more times. To establish MEG as a clinically reliable neuroimaging tool that can distinguish disease from healthy populations, the reliability of FC patterns should be explored from repeat scans. Up to date, only a few studies accessed the test-retest reliability of MEG/electroencephalography (EEG) FC (Jin et al., 2011;Hardmeier et al., 2014;Garcés et al., 2016) while only one study has quantified the test-retest reliability of FC estimates in the source-space MEG (Garcés et al., 2016). Colclough et al. (2016) attempted to report the reliability of every edge-weighted connections with a high number of connectivity estimators but using a split-half strategy from a large pool of subjects. Practically, the results cannot be adopted as reliability of static network metrics since the analysis involved single MEG scan recordings. However, no study has ever explored the reliability of both static and dynamic networks at the source space in MEG.
In the present study, we investigated the test-retest reliability of both static and dynamic FC measures derived from MEG resting-state data. For that purpose, we computed whole-brain FC for 40 subjects who were scanned twice with a 1-week testretest interval. For each subject and session, MEG-beamformed source activity was estimated and FC was computed between 90 brain areas. FC was estimated with the imaginary part of the phase locking index (iPLV) and the correlation of the amplitude envelope (CorEnv) in both static (SFC) and dynamic models (DFC) by adopting a sliding window approach (de Pasquale et al., 2010;Dimitriadis et al., 2010aDimitriadis et al., , 2012aDimitriadis et al., , 2013aDimitriadis et al., , 2015aDimitriadis et al., , 2016aDimitriadis et al., , 2017aDimitriadis and Salis, 2017c). Afterwards, statistical and topological filtering schemes were applied to both SFC and DFC to reveal the true topology (Dimitriadis et al., 2017a,b). For the SFC approach, we estimated well-known network metrics in a node-wise fashion and the reliability was accessed via correlation values between the two measurements and across the cohort. Graph-based reliability was assessed with a novel graph diffusion distance metric.
For the DFC approach, node-wise network metrics were estimated across experimental time. To explore spatiotemporally the derived network activity, we first designed a codebook of prototypical network microstates and then assigned each of the instantaneous connectivity patterns to the most similar code symbol (e.g., functional connectivity graph-FCG) (Dimitriadis et al., 2010a(Dimitriadis et al., , 2012a(Dimitriadis et al., , 2013a(Dimitriadis et al., ,b, 2015a(Dimitriadis et al., , 2016aDimitriadis and Salis, 2017c). A codebook is a set of prototypical functional connectivity graphs (FCGs). In this way, we derive a unique symbolic time series from each individual where each symbol corresponds to one of the predefined prototypical functional connectivity microstates (FCµstates). The evolution of these symbol-patterns encapsulates significant state transitions. Furthermore, the evolution of these FCµstates can be seen as a first order Markovian Chain (MC) that can be modeled representing an individualized state transition model of resting-state FCµstates. Fractional occupancy of each FCµstate, transition rates of FCµstates and MC models are the key features to explore the reliability of chronnectome in MEG source space. The group-consistency of subject-specific FCµstates was further explored. The whole analysis of dynamic functional connectivity graphs and the definition of FCµstates have been described in previous paper (Dimitriadis et al., 2013a).
Many techniques have already been proposed to summarize brain activity into short-lived transient brain states using the spectrum of neuromagnetic recordings (Vidaurre et al., 2016) and also the band-limited amplitude envelopes of source reconstructed MEG data (Baker et al., 2014;O'Neill et al., 2015). In detail, Vidaurre et al. (2016) proposed a combination of multivariate autoregressive model with hidden markovian modeling (MAR-HMM) in order to model the temporal, spectral and spatial properties of MEG reconstructed activity into very short-lived brain states. Similarly, (Baker et al., 2014) modeled resting-state source-reconstructed MEG activity with HMM into distinct spatio-temporal activation profiles called brain states. These brain states were linked to well-known anatomical brain areas. O'Neill et al. (2015) mined MEG source activity from two tasks, a self-paced motor and a Sternberg working memory task. He used a sliding window canonical correlation analysis (CCA) to estimate the functional connectivity at each time-window and a k-means clustering to detect repeatable spatial patterns of connectivity that form transiently synchronizing sub-networks (TSNs) or functional connectivity microstates. Here, we must underline the distinction of summarizing brain activity using the raw time series (ROIs × sliding windows; Baker et al., 2014;Vidaurre et al., 2016) which is a 2D matrix and the dynamic functional brain networks (ROIs × ROIs × sliding windows) which is a 3D matrix (O'Neill et al., 2015). Currently the mapping and the relationship between raw activity and brain connectivity and also the relationship of microstates (raw activity) with functional connectivity microstates (dynamic graphs; Allen et al., 2012;Dimitriadis et al., 2013aDimitriadis et al., ,b, 2015aDimitriadis et al., , 2016aDimitriadis and Salis, 2017c) is still unknown. Further research is needed to explore their mapping at resting-state and during tasks.
The proposed methodological scheme entails two distinct ways of analyzing dynamic functional connectivity patterns. These patterns are representative brain network topologies across subjects and brain rhythms and are directly linked to a brain state (Buzsáki and Draguhn, 2004). The very first approaches in fMRI constitutes novel contributions to an emerging neuroimaging field called chronnectomics (Allen et al., 2012;Calhoun et al., 2014). Previously, we reported the notion of FCµstates (Dimitriadis et al., 2013a) and the developmental trends in cognition (Dimitriadis et al., 2015a) using electroencephalographic recordings. The concept of chronnectome is the incorporation of a dynamic view of functional brain connectivity networks and the evolution of revisiting distinct spatio-temporal brain states (functional connectivity microstates-FCµstates). To the best of our knowledge, this study constitutes the first attempt to assess the test-retest reliability of Dynamic Functional Connectivity at the MEG source level.
Despite growing enthusiasm in the neuroscience community about the potential contribution of neuroimaging and especially brain networks in the designing of connectomic biomarkers for various brain diseases/disorders, many challenges remain open (Stam, 2014). At first level, it is more than significant to explore how reliable are network metrics at both temporal scales (static and dynamic) by analyzing a group of control subjects with repeat scans (e.g., diffusion MRI: Dimitriadis et al., 2017d). Here, we assess evidence of the reliability of neuromagnetic (MEG) based functional connectomics to lead to potential clinically meaningful biomarker identification in target populations through the lens of the criteria used to evaluate clinical tests.

MATERIALS AND METHODS
Subjects 40 healthy subjects (age 22.85 ± 3.74 years, 15 women and 25 men) underwent two resting-state MEG sessions (eyes open) with a 1-week test-retest interval. For each participant, scans were scheduled at the same day of the week and same time of the day. The duration of MEG resting-state was 5 min for every participant. The study was approved by the Ethics Committee of the School of Psychology at Cardiff University, and participants provided informed and written consent.

MEG-MRI Recordings
Whole-head MEG recordings were made using a 275-channel CTF radial gradiometer system. An additional 29 reference channels were recorded for noise cancelation purposes and the primary sensors were analyzed as synthetic third-order gradiometers (Vrba and Robinson, 2001). Two or three of the 275 channels were turned off due to excessive sensor noise (depending on time of acquisition). Subjects were seated upright in the magnetically shielded room. To achieve MRI/MEG coregistration, fiduciary markers were placed at fixed distances from three anatomical landmarks identifiable in the subject's anatomical MRI, and their locations were verified afterwards using high-resolution digital photographs. Head localization was performed before and after each recording, and a trigger was sent to the acquisition computer at relevant stimulus events.
All datasets were either acquired at or down-sampled to 600 Hz, and filtered with a 1-Hz high-pass and a 200-Hz lowpass filter. The data were first whitened and reduced in dimensionality using principal component analysis with a threshold set to 95% of the total variance (Delorme and Makeig, 2004). The statistical values of kurtosis, Rényi entropy and skewness of each independent component were used to eliminate ocular and cardiac artifacts. Specifically, a component was deemed artifactual if more than 20% of its values after normalization to zero-mean and unit-variance were outside the range of [−2, +2] (Delorme and Makeig, 2004;Escudero et al., 2011;Antonakakis et al., 2016). The artifact-free multichannel MEG resting-state recordings were then entered in the beamforming analysis (see next section).
Subjects further underwent an MRI session in which a 3T GE scanner with an eight-channel receive-only head RF coil T1weighted 1-mm anatomical scan was acquired, using an inversion recovery spoiled gradient echo acquisition.

Beamforming
An atlas-based beamformer approach was adopted to project MEG data from the sensor level to source space independently for each brain rhythm. The frequency bands studied were: δ (0.5-4 Hz), θ (4-8 Hz), α 1 (8-10 Hz), α 2 (10-13 Hz), β 1 (13-20 Hz), β 2 (20-30 Hz), γ 1 (30-45 Hz), γ 2 (55-90 Hz). First, the coregistered MRI was spatially normalized to a template MRI using SPM8 (Weiskopf et al., 2011). The AAL atlas was used to anatomically label the voxels, for each participant and session, in this template space (Tzourio-Mazoyer et al., 2002). The 90 cortical regions of interest (ROIs) were used for further analysis, as is common in recent studies (Hillebrand and Barnes, 2002;Hillebrand et al., 2016;Hunt et al., 2016). Next, neuronal activity in the atlas-labeled voxels was reconstructed using the LCMV source localization algorithm as implemented in Fieldtrip (Oostenveld et al., 2011). The beamformer sequentially reconstructs the activity for each voxel in a predefined grid covering the entire brain (spacing 6 mm) by weighting the contribution of each MEG sensor to a voxel's time series-a procedure that creates the spatial filters that can then project sensor activity to the cortical activity. Each ROI in the atlas contains many voxels, and the number of voxels per ROI differs. To obtain a single representative time series for every ROI, we defined a functional-centroid ROI representative by functionally interpolating activity from the voxel time series, within each ROI, in a weighted fashion. Specifically, we estimated a functional connectivity map between every pair of source time series within each of the AALs ROIs (Equation 1) using correlation (Equation 2). We then estimated the connectivity strength of each voxel within the ROI by summing its connectivity values to other voxels within the same ROI (Equation 3) and finally we normalized each strength by the sum of strengths (Equation 4) to estimate a set of weights within the ROI that sum to a value of 1. Finally, we multiplied each voxel time series with their respective weights and we summed across them in order to get a representative time series for each ROI (Equation 5). The whole procedure was applied independently to every quasi-stable temporal segment derived by the settings of temporal window and stepping criterion.
The following Equations 1-5 demonstrated the steps for this functional interpolation.
The outline of the methodology is described in Figure 1. An exemplar of the representative bandpass filtered ROI based time series is given in Figure 1. Figure 2 illustrates the preprocessing steps described in Equations (1-5).

Intra-Frequency Connectivity Estimators
Among the available connectivity estimators, we adopted one based on the imaginary part of phase-locking value (iPLV) (Lachaux et al., 1999) and adjusted properly so as to extract time-resolved profiles of intra-frequency coupling from MEG multichannel recordings at resting state. The original PLV is defined as follows: where k, l denote a pair of MEG sources and the imaginary part of PLV is equal to: The imaginary part of PLV (iPLV) investigates intra-frequency interactions without putative contributions from volume conductance. In general, the iPLV is mainly sensitive to nonzero-phase lags and for that reason is resistant to instantaneous self-interactions from volume conductance (Nolte et al., 2004). In contrast, it could be sensitive to phase changes that not necessarily imply a PLV oriented coupling. Correlation of the Envelope coupling (CorEnv) is based upon correlation between the oscillatory envelopes of two frequency band limited sources . See Figure 1 for a schematic diagram of phase and envelope based connectivity analyses based upon neural oscillations. Correlation of the Envelope coupling (CorEnv) is based upon correlation between the oscillatory envelopes of two band limited sources (A) while phase coupling searches for a constant phase lag between signals, FIGURE 1 | Outline of the methodology for accessing the reliability of network metrics derived from functional connectivity graphs (FCGs). (SF-Statistical Filtering, TF-Topological Filtering, FCE-Functional Connectivity Estimator). Statistical and topologically filtering of the FCGs will be described in sections Surrogate Data Analysis of iPLV/CorEnv Estimates-Statistical Filtering of Brain Networks and A data-driven Topological Filtering Scheme based on Orthogonal Minimal Spanning Trees (OMSTs), correspondingly. One can understand how from a full-weighted FCG, a more sparse version is derived via the statistical and topological filtering approaches. in the example a difference of π (B). The time series for the estimation of CorEnv were orthogonalized between each other using the bivariate version of this correction for signal leakage effects (Colclough et al., 2015).

Static Functional Connectivity Analysis
Using both connectivity estimators, we estimated the fullyweighted (90 × 90) anatomical oriented FCG, one for each subject, recording session and frequency band. To construct the static FCG (SFCG), we incorporated in the analysis the whole 5 min of the recording session.

Dynamic iPLV Estimates: The Time-Varying Integrated iPLV Graph ( TVI iPLV Graph)
The goal of the analytic procedures described in this section was to understand the repertoire of phase-to-phase interactions and their temporal evolution, while taking into account the quasiinstantaneous spatiotemporal distribution of iPLV estimates. This was achieved by computing one set of iPLV estimates within each of a series of sliding overlapping windows spanning the entire 5-min continuous MEG recording for eyes-open condition. The width of the temporal window and the stepping criterion were optimized for each frequency band separately using as objective criterion the reliability of transition dynamics between scan session 1 and 2 for each brain rhythm (Dimitriadis et al., 2013a; see sections State Transition Rate and Optimizing the Width of the Time-Window and the Stepping Criterion). The center of the stepping window moved forwards every frequencydependent time-window (see sections Optimizing the Width of the Time-Window and the Stepping Criterion and Tuning Parameters for Dynamic Functional Connectivity Analysis) for the optimization of the parameters) for every intra-frequency interactions and a new functional brain network is re-estimated between every pair of "swifting" temporal segments of MEG activity, from two sources, leading to a "quasi-stable in time" static iPLV graph. In this manner, a series of 598 (for δ) to 2,140 (for γ 2 ) iPLV graph estimates were computed for each frequency (8 within frequency), for each participant and for both repeat scans.
For each subject, a 4D tensor (frequencies bands (8) × slides (598-2,140) × sources (90) × sources (90); see sections Optimizing the Width of the Time-Window and the Stepping Criterion and Tuning Parameters for Dynamic Functional Connectivity Analysis) was created for each condition integrating subject-specific spatio-temporal phase interactions ( Figure 3A).

Estimates-Statistical Filtering of Brain Networks
To identify significant iPLV/CorEnv-interactions which were estimated for every pair of frequencies within and between all 90 sources, and at each successive sliding window (i.e., temporal segment), we employed a surrogate data analysis. Accordingly, we could determine (a) if a given iPLV/CorEnv value differed from what would be expected by chance alone, and (b) if a non-zero iPLV/CorEnv corresponded to non-spurious coupling.
For every temporal segment, sensor-pair and frequency, we tested the null hypothesis H 0 : "the observed iPLV/CorEnv value comes from the same distribution as the distribution of surrogate iPLV/CorEnv-values". One thousand surrogate time-series were generated by cutting at a single point at a random location the original time series and exchanging the two resulting time courses (Aru et al., 2015). We restricted the range of the selected cutting point in a temporal window of width up to 10 s in apart from the middle of the recording session (between 140 and 160 s). This surrogate scheme was applied to the original whole time series and not to the signal-segment at every slide. Repeating this procedure leads to a set of surrogates with a minimal distortion of the original phase dynamics, while the non-stationarity of the brain activity is less destroyed compared to shuffling the time series or cutting and rebuilding it in more than one time points.
This procedure ensures that the real and surrogate indices both have the same statistical properties. For each data set, the surrogate iPLV/CorEnv ( s iPLV/ s CorEnv) was then computed .
We then determined a one-sided p-value for each iPLV/CorEnv value that corresponded to the likelihood that the observed value could belong to the surrogate distribution. This was done by directly estimating the proportion of "surrogate" s iPLV/ s CorEnv that was higher than the observed iPLV/CorEnv. The p-value reflected the statistical significance of the observed iPLV/CorEnvlevel (a very low value revealed that it could not have appeared from processes with no iPLV coupling).
At a second level, we applied the FDR method (Benjamini and Hochberg, 1995) to control for multiple comparisons within each snapshot of the dynamic graph (FCG-a 90 × 90 matrix with tabulated p-values) with the expected fraction of false positives set to q ≤ 0.01. Finally, for each subject the resulting TV iPLV/ TV CorEnv profiles constituted of two 3D arrays of size [598 to 2,140 for δ to γ 2 (time windows) × 90 (sources) × 90 (sources)] with a value of 0 indicated a non-significant iPLV/CorEnv value.
The aforementioned statistical filtering approach was applied independently for each frequency band, session, subject, and connectivity estimator for both static and dynamic functional connectivity graphs.

A Data-Driven Topological Filtering Scheme Based on Orthogonal Minimal Spanning Trees (OMSTs)
As well as the statistical filtering approach, it is important to adopt a data-driven topological filtering approach in order to reveal the backbone of the network topology over the increment of information flow.
Recently, it was proved that MST is an unbiased method that yields reliable network metrics (Tewarie et al., 2015). In this study, we adopt a variant of this topological filtering scheme called orthogonal minimal spanning trees (OMST), which leads to a better sampling of brain networks, preserving the advantage of MST, that connects the whole network with minimum cost without introducing cycles and without differentiated strong from weak connections compared to the absolute threshold or the density threshold (Dimitriadis et al., 2017a,b). MST is too sparse to capture the "true" network and for that reason leading to the selection of N-1 connections where N denotes the number of nodes. We introduced OMST which samples the weights of a brain network via the notion of MST and under the optimization of the global information flow under the constraint of the total Cost of preserving the functional connections (Dimitriadis et al., 2017b,d).
Our criterion for topologically filtering a given brain network is based on the maximum value of the following quality formula: We applied the data-driven topological filtering scheme based on OMST at every static and quasi-instantaneous FCG from the dynamic DFCG. After statistical and topological filtering approaches applied to both SFCG and the DFCG, we estimated network metrics at the node/source level. One can understand that human brain demonstrates a preferred transition from FCµstates 2 to FCµstates 1 compared to the opposite direction (see 2D colormap). The chronnectomics were derived from this symbolic time series. ED, Euclidean distance; LEG, Laplacian EiGenvalues; ED, Euclide and distance; AAL, automated anatomical labeling; LEG, Laplacian EiGenvalues. Figure 1 demonstrates an example of a full-weighted FCG after applying both statistical and topological filtering approach. Our algorithm was validated over all the existing thresholding schemes with a large EEG dataset over brain fingerprinting and with a multi-scan fMRI dataset over reliability of nodal network metrics (Dimitriadis et al., 2017a). Additionally, we demonstrated the importance of a data-driven topological filtering technique in functional neuroimaging by using OMST in a multi-group study with MEG resting-state recordings (Dimitriadis et al., 2017b). The MATLAB code of the OMST method and also of the majority of existing filtering methods can be downloaded from the: https://github.com/stdimitr/ topological_filtering_networks & researchgate https://www. researchgate.net/profile/Stavros_Dimitriadis.

Graph Diffusion Distance Metric for Brain Networks
In order to assess group and scan sessions differences in the topologically filtered FCG at the single-case level, we computed the Graph Diffusion Distance as a distance metric (Fouss et al., 2012;Hammond et al., 2013) from the OMST-derived final Functional Connectivity Graphs (FCG). The graph laplacian operator of each subject-specific FCG was defined as L = D − FCG, where D is a diagonal degree matrix related to FCG. This method entails modeling hypothetical patterns of information flow among sources based on each observed (static) SFCG. The diffusion process on the person-specific FCG was allowed for a set time t; the quantity that underwent diffusion at each time point is represented by the time-varying vectoru(t) ∈ ℜ N . Thus, for a pair of sources i and j, the quantity FCG ij (u i (t) − u j (t)) represents the hypothetical flow of information from i to j via the edges that connect them (both directly and indirectly). Summing all these hypothetical interactions for each sensor leads to u ′ j (t) = i FCG ij (u i (t) − u j (t)), which can be written as: where L is the graph laplacian of FCG. At time t = 0 Equation 9 has the analytic solution: u(t) = exp(−tL)u (0) . Here exp(-tL) is a N × N matrix function of t, known as Laplacian exponential diffusion kernel (Fouss et al., 2012), and u (0) = e j , where e j ∈ ℜ N is the unit vector with all zeros except in the jth component. Running the diffusion process through time t produced the diffusion pattern exp(-tL)e j which corresponds to the jth column of exp(-tL). Next, a metric of dissimilarity between every possible pair of person-specific diffusion-kernelized FCGs (FCG 1 , FCG 2) was computed in the form of the graph diffusion distance d gdd (t). The higher the value of d gdd (t) between two graphs, the more distinct is their network topology as well as the corresponding, hypothetical information flow. The columns of the Laplacian exponential kernels, exp(-tL1) and exp(-tL2), describe distinct diffusion patterns, centered at two corresponding sources within each FCG. The d gdd (t) function is searching for a diffusion time t that maximizes the Frobenius norm of the sum of squared differences between these patterns, summed over all sources, and is computed as: where . F is the Frobenius norm.
Given the spectral decomposition L = V V, the laplacian exponential can be estimated via where for , exp(-t ) is diagonal to the ith entry given by e −t i,i . We computed d gdd (FCG 1 ,FCG 2 ) by first diagonalizing L1 and L2 and then applying Equations (9, 10) to estimate d gdd (t) for each time point t of the diffusion process. In this manner, a single dissimilarity value was computed for each pair of participants based on their individual characteristic FCGs. For further details see (Hammond et al., 2013). The GDD metric can be downloaded from: https://github.com/stdimitr/multi-group-analysis-OMST-GDD.

Static Network Metrics
After applying the statistical and topological filtering approach, we estimated the global efficiency for each node in static approach. The static approach leads to 90 (sources) values for each network, metric, frequency band and session per subject. We adopted complementary features that measure the importance of each node in segregation, integration and the information flow within a weighted functional brain network (Dimitriadis et al., 2010a(Dimitriadis et al., ,b, 2013a(Dimitriadis et al., ,b, 2015a. In this study, we estimated four basic network metrics, the global and local efficiency, the strength of each node and the mean first passage time based on random walks (Goñi et al., 2013).
Network global efficiency (GE) reflects the overall efficiency of parallel information transfer within the entire set of 90 sources and was estimated as the average source specific GE value over all sources using the following formula (Latora and Marchiori, 2001): where d denotes the shortest path length from i to j. Local efficiency (LE) indicates how well the subgraphs exchange information when a particular node is eliminated (Achard and Bullmore, 2007). Specifically, each node is assigned the shortest path length within its subgraph G i G i where k i corresponds to the total number of spatial (first level neighbors) neighbors of the ii-th node, while d denotes shortest path length.
The strength is equal to the total sum of the weights of the connections of each node.
As a fourth candidate network metric, we adopted the mean first passage time (MFPT). Starting a random walk process on a brain network, an analytic expression can give the probability that a single particle departing from a node i arrives at node j for the first time within exactly L steps (Wang and Pei, 2008). This criterion can be applied for each MEG source pair by setting L to their shortest-path-length. We denote with G = π ij the n × n symmetric matrix containing, for each pair of nodes, the probability of a single particle going from node i to node j via the shortest path. Each entry π ij can be computed as where matrix B j is the transition matrix P introduced above, but with all zeros in the j-th column, i.e., with j acting as an absorbing state (Wang and Pei, 2008). Evaluating shortest-pathlengths ensures that ∀i, j π ij > 0. By considering one particle here, the average shortest-path probability of a graph is defined as The derived 2D matrix based on nodal NMTS of GE, LE, MFPT and strength will be modeled with the proposed method that is described in the following section.

Modeling of Dynamic Functional Connectivity Graphs (DFCG) as a 3D Tensor
This subsection serves as a brief introduction to our symbolization scheme, presented in greater details elsewhere (Dimitriadis et al., 2012a(Dimitriadis et al., ,b, 2013a. The dynamic functional connectivity patterns can be modeled as prototypical functional connectivity microstates (FCµstates). In a recent study, we demonstrated a better modeling of dynamic functional connectivity graphs (DFCG) based on a vector quantization approach (Dimitriadis et al., 2013a). In our previous work (Dimitriadis et al., 2013a(Dimitriadis et al., ,b, 2015a, we used the neural-gas Each subplot illustrates the (dis)similarities of static FCGs across scanning sessions and subjects. The 2D matrix demonstrates the (dis)similarities of the static FCGs across the subjects and both repeat scans. Scanning sessions were coded with blue and red circles correspondingly and a black line connects the FCG of each subject between the two scanning sessions. With this representation one can read out the similarity of a static FCG between two scanning sessions and participants.
Stress expresses the loss of information expressed in the projected Frequency-Dependent Static Functional Connectivity Graphs in 2D feature space from an original 80D space. The low stress values mean that the relationship of the 80 FCGs in the original 80 × 80 matrix is preserved in the projected 2D space. R 1,2 refer to the 2D projected space of the 80 FCGs. FCG, functional connectivity graph; gDD, graph diffusion distance.
algorithm (Martinetz et al., 1993) to learn the 2D matrix (vectorized version of 2D matrix × time) leading to a codebook of k prototypical functional connectivity states (i.e., FCµstates). This algorithm is an artificial neural network model, which converges efficiently to a small number k of codebook vectors, using a stochastic gradient descent procedure with a softmax adaptation rule that minimizes the average distortion error (Martinetz et al., 1993). In a recent study, we adopted non-negative matrix factorization (NNMF) as an appropriate learning algorithm of the 2D vectorized version of a dynamic functional brain network (Marimpis et al., 2016). In our previous study where we first demonstrated how to model dynamic functional connectivity graph (dFCG) (Dimitriadis et al., 2013a), we vectorised the upper triangular of each of the quasi-static FCGs building a 2D matrix where the 1st dimension is the number of temporal segments and the 2nd the vectorised version of a static FCG. The final outcome of this approach is to define the so-called functional connectivity microstates (FCµstates). In a next study, we moved one step further by estimating node-wise global efficiency as the best descriptor to characterize the brain activity. The final outcome of the modeling using the same methodology of neural-gas algorithm was task-based network microstates (Dimitriadis et al., 2015a). Here, the vectorized version of a 90 × 90 FCG produces a long vector of 4,005 values while the number of temporal segments ranged from 598 to 2,140 which caused the so-called curse of dimensionality where the number of number of the temporal segment over which the modeling will learn the brain states is much smaller compared to the vectorized snapshot of FCG. Simultaneously, the vectorized notion of a brain network didn't maintain the inherent format of a functional brain network which is a 2D matrix, a tensor.
The outline of this procedure is illustrated in Figure 3. In Figure 3A, a characteristic bandpass filtered time series for each of the studying frequency band is estimated from each ROI. Here, instead of vectorising the upper triangular of an undirected FCG, we used the statistical and topological filtering FCG on its inherent format which is a 2D tensor. In the case of dynamic networks, the dimension is a 3D tensor where the 3rd dimension is the time. Figure 3B illustrates a few snapshots of the dFCG for δ frequency of the first subject. At the next level, we estimated the laplacian matrix of each quasi-static FCG. Given a FCG, the laplacian matrix is given by: where D is the degree matrix and A is the FCG.  Figure 3C demonstrates the laplacian matrix of the FCG in Figure 3B. In the main diagonal, the degree of each node is tabulated. Afterward, we applied an eigenanalysis for each of these laplacian matrixes and the eigenvalues of this procedure describes the synchronizability of the original FCG. Figure 3D illustrates for the 1st min the eigenvalues for each quasi-static FCG. One can easily detect the abrupt transition between the brain states. Here, neural-gas algorithm was applied on the 2D matrix presented in Figure 3D after first concatenated across subjects independently for each frequency band and scan session. The main scope of this codebook learning algorithm is to define FCµstates.
By estimating the reconstruction error E between the original 2D matrices (90 × [slides × subjects]) and the one reconstructed via the k FCµstates assigned to each snapshot of the DFCG for each predefined threshold, we can detect the optimal threshold T for each case. In this work, the criterion of the reconstruction error E was set less than 4%. Practically for all the frequency bands and in both connectivity estimators, the reconstruction error E was less than 2%. The selected threshold was detected based on the plateau by plotting of reconstruction error E vs. the threshold T.
In this way, the richness of information contained in the dynamic connectivity patterns is represented, by a partition matrix U, with elements u ij indicating the assignment of input connectivity patterns to code vectors. Following the inverse procedure, we can rebuild a given time series from the k FCµstates, with a small reconstruction error E. The selection of parameter k reflects the trade-off between fidelity and compression level. As a consequence, the symbolic time series closely follows the underlying functional connectivity dynamics. The derived symbolic times series that keep the information of network FCµstates (nFCµstates) are called hereafter as STS L−EIGEN (L:Laplacian − Eigen:Eigenalysis). Figure 3E tabulates the correlation of the eigenvalues between every pair of temporal segments while in Figure 3F, the matrix in Figure 3E was reordered such as the FCµstates to be revealed via the neural-gas algorithm. The network topology of the extracted FCµstates is illustrated in Figure 3G. From, Figure 3F, one can understand that the two FCµstates describe the DCFG of this subject.
An exemplar of prototypical FCµstates is illustrated in Figure 3G. The outcome of this clustering procedure is also to extract a symbolic time series per subject, repeat scan and frequency that describes the transition of brain activity between the extracted brain states (FCµstates; Figure 3D). The transition probability P for this example and for the two FCµstates is illustrated with a classical figure for first order Markovian Chain. The self-arrows refer to the percentage of sliding windows where the brain stays stable in a FCµstate without any transition while the directed arrow gives the percentage of transition from one FCµstate to the other.This symbolic time series can be seen as FIGURE 6 | Reliability of node-wise network metrics derived from static brain networks with iPLV connectivity estimator. Each subplot demonstrates the correlation coefficient (CC) of each network metric at every studying frequency band of each brain area between the two scanning sessions. CC, the correlation coefficient; AAL, automated anatomical labeling. a first order Markovian chain where these switching between "quasi-static" FCµstates can be modeled as a finite Markov chain (Dimitriadis et al., 2013a(Dimitriadis et al., ,b, 2015aO'Neill et al., 2015;Vidaurre et al., 2016). One can clearly understand that human brain demonstrates a preferred transition from FCµstates 2 to FCµstates 1 (off-diagonal lines of the TP) compared to the opposite direction ( Figure 3H). The sketch of the markovian chain and the colored TP matrix can reveal the aforementioned trend of preferred direction FCµstates 2 to FCµstates 1 .
From the symbolic timeseries, specific metrics tailored to the dynamic evolution of FCµstates were estimated (see next section) and their reliability was assessed via the correlation coefficient between scan session 1 and scan session 2. The whole approach was repeated independently for each frequency band and connectivity estimator by integrating subject and scan-based DFCG.

Characterization of Time-Varying Connectivity
Once the integrated DFCG is formed and it is modeled via the combination of neural-gas and laplacian eigenanalysis scheme [N-GAS L−EIGEN (L:Laplacian -Eigen:Eigenalysis)], relevant features can be extracted from the data based on the statetransition states. There features are called chronnectomics (chronos-Greek word for time and connectomics for network metrics) which are described in the following section.

Chronnectomics
The following chronnectomics (dynamic network metrics) will be estimated on the STS L−EIGEN which expresses the fluctuation of the FCµstates.

State Transition Rate
Based on the state transition vectors STS L−EIGEN as demonstrated in Figure 3A, we estimated the transition rate (TR) for every pair of states as followed: where slides denote the number of temporal segments using the sliding window approach. TR yields higher values for increased numbers of "jumps" of the brain between the derived brain states over consecutive time windows. This approach leads to one feature per participant.

Occupancy Times of the nFCµstates
Complementary to the aforementioned chronnectomics, we estimated also the occupancy time (OC) of each FCµstates as the percentage of its occurrence across the experimental time. OC was estimated from STS L−EIGEN as follows: where k denotes the FCµstates.
FIGURE 7 | Reliability of node-wise network metrics derived from static brain networks with CorEnv connectivity estimator. Each subplot demonstrates the correlation coefficient (CC) of each network metric at every studying frequency band of each brain area between the two scanning sessions. CC, the correlation coefficient; AAL, automated anatomical labeling.

Reliability of Static Network Metrics and Chronnectomics
The reliability of static node-wise network metrics and the chronnectomics was assessed with the correlation coefficient between forty values derived from scan session 1 and forty values from scan session 2 for each frequency band, condition and connectivity estimator (see Figures 1-3).

Optimizing the Width of the Time-Window and the Stepping Criterion
We optimized both the width of the time-window and the step criteria for the sliding-window approach based on the maximization of the reliability of TR. The reliability was estimated based on the correlation coefficient of the TR across the whole group between scan session 1 and 2. The whole procedure was followed independently for each brain rhythm. The settings for the width of temporal window and the step were defined as a percentage of the cycles of the studying frequencies: {from 1 up to 10 cycles with step equals to 0.5 cycle} for the width of the temporal window and {from 0.1 cycles to 2 cycles with step equals to 0.1 cycle} for the step (see Table 1).
To avoid overfitting of both TR and OT since, we used TR for both the optimization of the width of the temporal window and the stepping criterion, we used the optimized parameters in an external second repeat scan dataset for further evaluation (see Supplementary Material).

Tuning Parameters for Dynamic Functional Connectivity Analysis
The optimization of the temporal window and the stepping criterion for each brain rhythm reveals a nice trend for dynamic functional connectivity analysis. The width of the temporal window increased from δ to γ 2 while the stepping criterion decreased in both connectivity estimators (see Table 1).
FIGURE 8 | Frequency-dependent FCGµstates iPLV . Network topologies of the FCGµstates iPLV for each of the studying frequency band. To enhance the visualization and contrast of FCGµstates across frequency bands, we adopted a network-wise representation instead of plotting the brain network in a 90 nodes layout. The 90 ROIs of the AAL template were assigned to each of the five networks: default-mode (DMN), fronto-parietal (FPN), occipital (O), sensorimotor (SM), and cingulo-opercular (CO). The color of each node denotes total strength of within network connections while the color of each line the total strength of between network connections. Both strength values were normalized across both frequencies and FCGµstates.

Common Projection Space of Frequency-Dependent Static FCG
To demonstrate the (dis)similarities between sessions and subjects of the frequency-dependent static FCG, we constructed a distance matrix of dimensions 80 × 80 (subjects × sessions) using the graph diffusion distance metric. Then, we applied multidimensional scaling (MDS) to project the distance matrix into a common 2D feature space. Using a different colored circle for each scanning session (blue for session 1 and red for session 2) and connecting both of them with a black line for each subject, we further enhanced the (dis)similarities of the static FCGs. Figures 4, 5 illustrate these FCG-based projections for static FCG Iplv and FCG CorEnv correspondingly. In Figure 4G one can detect a few subjects with high reliable static FCG between the two scan sessions and also subject-specific network topologies that occupied an isolated subarea in the common projection FCG space. The stress index estimated via the MDS approach was low and the relationship of the 80 FCGs in the original 80 × 80 matrix is preserved in the projected 2D space.

Reliability of Static Network Metrics
Figures 6, 7 demonstrate the correlation coefficients for each node-wise network metric between the two scanning sessions for every frequency-dependent static FCG. From the visual comparison of both figures one can clearly reveal that the correlation values are higher for CorEnv compared to iPLV.
Applying Wilcoxon Rank-Sum Test for every frequency and network metric between the 90 correlation values, we detected statistical significant differences in every case (p < 0.01, p ′ < p/32, Bonferroni Corrected). However, the averaged correlation values did not reach high reliability (e.g., >0.9) even for the CorEnv. It is obvious from the correlation plots that the reliability of nodewise static network metrics has high spatial variability in both connectivity estimators.

Frequency-Dependent FCGµstates and Reliability Chronnectomics for iPLV
Our analysis of DFCG based on iPLV revealed two FCGµstates iPLV for each frequency band. The topology of these frequency-dependent FCGµstates iPLV is illustrated in Figure 8. We integrated the nodes into five well-known brain networks: default-mode (DMN), fronto-parietal (FPN), occipital (O), sensorimotor (SM), and cingulo-opercular (CO). The mapping between the 90 ROIs of AAL and the five brain networks can be retrieved from section Results in Supplementary Material. One can clearly detect that the functional coupling between the default mode network and the cingulo-opercular dominates the coupling strength across the frequency bands and FCGµstates with less pronounced effect in both γ bands. Complementary, the coupling strength between and within the networks is diminished after α 2 frequency. This behavior can be interpreted as a reduction of the connections up to the defined threshold following the increment of the frequency. The two FCGµstates iPLV showed also a different distribution of strength globally and locally. Both types of chronnectomics, transition rates (TR) (Figure 9) and occupancy times (OC) (Figure 10) demonstrated high reliability (Corr. > 0.9, p < 10 −7 ) across frequency bands. Similar results, we obtained also for the second external dataset (see section 2 in Suplementary Material).

Frequency-Dependent FCGµstates and Reliability Chronnectomics for CorEnv
Our analysis of DFCG based on the correlation of the envelope connectivity estimators revealed two FCGµstates CorEnv for each frequency band. The topology of these frequency-dependent FCGµstates CorEnv is illustrated in Figure 11. The mapping between the 90 ROIs of AAL and the five brain networks can be retrieved from section 2 in Supplementary Material. One can clearly detect that the functional coupling between the default mode network and the cingulo-opercular dominates the coupling strength across the frequency bands and FCµstates. Complementary, the coupling strength between and within the networks is diminished after α 2 frequency as it was observed for FCGµstates iPLV . Complementarily, the network topologies of FCGµstates CorEnv between low and high frequencies based on the strength coupling are more common than the FCGµstates iPLV .This common substrate across the FCGµstates CorEnv is consistent with the general notion that correlation of the envelope is more stacked to the structural connectome compared to the phase-based connectivity patterns which demonstrate higher degrees of freedom (Engel et al., 2013; compare Figure 8 with Figure 11).
Only transition rates (TR) showed high reliability for CorEnv (Corr. > 0.8, p < 10 −4 ) in all the frequency bands with the only exception of β 1 (Figure 12). Occupancy times (OT) showed low reliability across the frequency bands (p > 0.4) (Figure 13). TR of FCGµstates CorEnv increased with the increment of frequency reaching a plateau in γ 1 . In contrast, TR of FCGµstates iPLV did not show such a frequency-dependent behavior. Similar results, we obtained also for the second external dataset (see section 2 in Supplementary Material).

DISCUSSION
In the present study, we assessed the reliability of both static and dynamic functional connectivity network descriptors using resting-state MEG data from 40 subjects with repeat scan sessions. This is the first time that the reliability of chronnectomics, at least for the MEG modality, has been taken into account. Source time series were first beamformed  independently for each frequency band (Hillebrand et al., 2005;Schoffelen and Gross, 2009;Brookes et al., 2011b), and then representative voxel time series based on the AAL atlas were extracted using a novel linear interpolation analysis. This procedure produces informative timeseries with a characteristic carrier frequency compared to the noisy time series derived by PCA or by selecting the voxel time series within a ROI that encapsulates the maximum power. Then, both static and dynamic frequency-dependent functional connectivity graphs were computed for each subject and scan session using the imaginary part of phase locking value (iPLV) and the correlation of the amplitude envelope (CorEnv). Both static and dynamic FCG (SFCG-DFCG) were filtered both statistically (surrogates) and topologically (OMST; Dimitriadis et al., 2017a,b).
Here, we adopted a data-driven pipeline of how to estimate both static and dynamic FCG statistically and topologically filtered using an algorithm previously applied to EEG recordings. We explored the reliability of both static network metrics and chronnectomics (dynamic network metrics) by employing two representative connectivity estimators for the construction of static and dynamic brain networks. Using this pipeline, prototypical FCµstates were derived which were highly reproducible across subjects and scan sessions in both connectivity estimators and in all frequencies. The reliability of node-wise static network metrics based on four network metrics was low and spatially variable with both connectivity estimators while the CorEnv demonstrates higher ICC values compared to iPLV. The reliability of chronnectomics (TR, OT) for iPLV was high while for CorEnv the reliability of only the TR reaches high acceptable levels. Our results were reproduced also in a second external dataset (see Supplementary Material). Our study strongly encouraging the use of DFCG with neuromagnetic recordings that takes the advantage of the nature of MEG modality, its high temporal resolution.
To our knowledge, this is the very first study that explored the reliability of both static and dynamic FCG and the related network metrics and chronnectomics, respectively in neuromagnetic source space at. In static FCG, node-wise network metrics demonstrated poor reliability for iPLV and poor to medium for CorEnv. The node-wise reliability was highly spatial variable and static FCG have also demonstrated low repeatability in both connectivity estimators and especially in CorEnv. In contrast, prototypical FCµstates were high reproducible across subjects and scan sessions in both connectivity estimators and in all frequencies supporting by the low reconstruction error (<2%) of our brain network learning algorithm. Complementary, the reliability of chronnectomics (TR,OT) for iPLV was high while for CorEnv the reliability of only the TR reaches high acceptable levels. These results strongly encourages the neuroscientists to adopt the notion of DFCG with neuromagnetic recordings that takes the advantage of its high temporal resolution.
Two main studies explored the reliability of static FCG on the source level using MEG-beamformed resting-state connectivity analysis. Garcés et al. (2016) studied the reliability of restingstate networks using four connectivity estimators: phase-locking value (PLV), phase lag index (PLI), direct envelope correlation (d-ecor), and envelope correlation with leakage correction (lcecor). They adopted intra-class correlation coefficient (ICC) and Kendall's W for assessing within and between-subjects agreement respectively. Higher test-retest reliability was found for PLV from θ to γ, and for lc-ecor and d-ecor in β. They commented that high ICC in PLV and d-ecor could be artifactual due to volume conduction effects. Colclough et al. (2016) investigated the reliability of static FCG at resting-state using beamformed source static connectivity analysis. They reported high reliability mostly for the partial correlation analysis and the correlation of the envelope among 12 connectivity estimators. Two more studies, Deuker et al. (2009) estimated the reliability of restingstate network metrics derived from MEG in sensor space using mutual information. They obtained high ICC for clustering, global efficiency and strength at a network level. Jin et al. (2011) found medium ICC for nodal global efficiency, nodal degree and betweenness centrality in α and β bands.
Our results revealed that nodal network metrics derived from static FCG are less reproducible then their dynamic counterparts. In contrast, chronnectomics are highly reproducible with both adopted connectivity estimators. These results complemented with the results presented in (Colclough et al., 2016) where they adopted multiple connectivity estimators for the construction of static brain networks on the source level using MEGbeamformed resting-state activity. Colclough et al. (2016) showed that the static full-weighted FCG are high repeatable within the group-level mostly for the correlation of the envelope adopting a split-half strategy on a dataset with single scans. Here, we accessed the reliability of any network metric using a two scan session design per subject. We should state here that edge-weights are significant for the construction of network topology and the reliability of connectomic biomarkers .
One of the key findings of our analysis are the frequencydependent FCµstates for each connectivity estimator. Figures 8, 11 illustrate the strength of the coupling within and between brain networks for the prototypical FCµstates at every frequency band. It is obvious that the highest strength within a network is observed within the DMN in both connectivity estimators. The strength between the brain networks is mainly distributed between DMN and the rest of the networks demonstrating the highest value till α 2 and dropped from β 1 to γ 2 (Figures 8, 11). DMN reignited a high interest the last years for the description of intrinsic ongoing activity in studies of the human brain in health and disease (Raichle, 2015). Disruptions of functional connections within the DMN and between DMN and the rest of brain networks have been linked to various neurological and neuropsychiatric disorders (Mohan et al., 2016). Studies in healthy aging and Alzheimer's disease have revealed the significant role of DMN (Mevel et al., 2011).
Flexible hub theory based on clustering analysis of functional networks gave an explanation of how temporal functional modes exist where one neural region may be switched from a certain network at one time to a different network at another time (Smith et al., 2012). It remains still unclear how the different brain networks are connected together during spontaneous and task-related activity. Dosenbach et al. (2008) proposed that the FPN may serve to initiate and adjust cognitive control, whereas another control-type network, the CO network (CON), provides stable set-maintenance. Cole and colleagues (Cole et al., 2013) helped to untangle the flexible role of the FPN, many questions remain regarding the interaction between the FPN and the CON and also with other networks such as the DMN,SM and O. In the present study, we characterized the dynamic relationships of the brain networks across time at resting-state in various frequency bands and using representative connectivity estimators. We found that these functional patterns are high reproducible which will help multi groups worldwide to explore these interactions and build reproducible connectomic biomarkers in various diseases and disorders. Understanding the neural basis of intrinsic activity, cognition and structure-function relationships, will further enhance the prognostic/diagnostic abilities in clinical populations.
The interactions of large-scale brain networks at resting-state and during tasks is characterized by the studying frequency. Frequency-specific functional interactions between large-scale brain networks may be individual fingerprints of the idle activity and cognition . It will be interesting in the future to explore how the FCµstates from a dynamic integrated functional connectivity graph (Dimitriadis and Salis, 2017c) that incorporates different intrinsic coupling modes both intra and cross-frequency coupling can be used for brain fingerprinting (Engel et al., 2013).
It is critically important to take advantage of new imaging modalities to untangle the mechanisms that produce circuit dysfunctions in many brain diseases and disorders. The development of biomarkers is very important and for that reason the proposed experimental paradigm and analytics of the metadata derived from the analysis of human brain activity must be highly reliable and reproducible. Magnetoencephalography (MEG) allows us to measure neuronal events noninvasively with millisecond resolution and recent advanced methods opens new avenues of exploring and answering fundamental key research questions tailored to each brain disease/disorder. MEG can become a pioneering clinical research tool for mental disorders (Bowyer et al., 2015;Grent-'T-Jong et al., 2016;Uhlhaas et al, 2017), Alzheimer's disease (López et al., 2014(López et al., , 2017Koelewijn et al., 2017), dyslexia (Dimitriadis et al., 2010b(Dimitriadis et al., , 2013b, traumatic brain injury (Dimitriadis et al., 2015c;Antonakakis et al., 2016Antonakakis et al., , 2017, multiple sclerosis (Tewarie et al., 2015), and other brain diseases. To establish MEG-based biomarkers that can be used for daily clinical practice and clinical evaluation, their reproducibility should be further explored. Complementary, the transition rate and also the occupancy times could be personalized biomarkers of a subject's resting-state condition where more task-related FCµstates and the related markers derived from them could build a subject specific database for longitudinal studies. Transition rates could be also correlated with IQ scores and also with behavioral performance during execution of cognitive tasks.
In the present study, we proposed a data-driven analytic pathway to assess the reliability of connectomics using MEGbeamformed connectivity analysis. Our results clearly support the notion of dynamic functional connectivity on the source level, the derived prototypical FCµstates and the related chronnectomics. Last years, many studies explored the dynamic functional connectivity graphs in many modalities (EEG/MEG/fMRI) and in both resting-state and during tasks (Dimitriadis et al., 2010a(Dimitriadis et al., , 2012a(Dimitriadis et al., ,c, 2013a(Dimitriadis et al., ,b, 2015a(Dimitriadis et al., ,b,c,d, 2016aBassett et al., 2011;Rosazza and Minati, 2011;Allen et al., 2012;Handwerker et al., 2012;Ioannides et al., 2012;Hutchison et al., 2013;Liu and Duyn, 2013;Braun et al., 2015b;Mylonas et al., 2015;Toppi et al., 2015;Yang and Lin, 2015;Calhoun and Adali, 2016, for reviews see Calhoun et al., 2014). This is the very first study according to authors' knowledge that the reliability of chronnectomics was explored. The outcome of this study opens new avenues in the exploration of human brain dynamics with MEG-beamformed source activity and under the notion of dynamic functional connectivity.
We addressed the key question of the readiness of neuromagnetic-based based functional connectomics to lead to clinically meaningful biomarker identification through the reliability approach that offers a repeat scan study in healthy controls. It is more than significant to customize stable approaches for analyzing neuromagnetic recordings and present reproducible brain connectomics across scans in healthy control populations without sacrificing the individual characteristics that can be used for personalized intervention neuroscience (Gratton et al., 2018). It is highly recommend to access the reliability of any metric derived from any neuroimaging modality in a repeat scan protocol in healthy control population before applying it to a larger disease group where the cost of scanning is too high (diffusion MRI: Dimitriadis et al., 2017d). Additionally, we will expand this analysis in future efforts to identify disease status alone including clinical variables related to genetic risk (Lancaster et al., 2018), expected treatment response and prognosis.

CONCLUSIONS
In conclusion, we provided the first source-space test-retest reliability of dynamic functional connectivity of neuromagnetic recordings at resting-state. We computed both static and dynamic functional connectivity based on 90 ROIs according to AAL templated and using two connectivity estimators, the iPLV and the CorEnv. Nodal network metrics were unreliable in both connectivity estimators but with higher reliability demonstrated for CorEnv. Moreover, their reliability demonstrates highly spatial variability. Static FCG were also unreliable and especially for CorEnv. In contrast, prototypical FCµstates were reliable in both connectivity estimators and across frequency bands. The derived chronnectomics (TR, OT) were highly reproducible for iPLV while only TR was reliable for CorEnv within acceptable levels. Our results strongly encourages future studies with main scope to explore resting-state networks in both healthy control and disease populations to apply a data-driven dynamic functional connectivity analysis using MEG-beamformed source reconstructed brain activity.

AUTHOR CONTRIBUTIONS
SD: conception of the research analysis, methods and design, data analysis, and drafting the manuscript; BR: data acquisition; BR, DL, and KS: critical revision of the manuscript. Every author read and approved the final version of the manuscript.

ACKNOWLEDGMENTS
SD and DL were supported by MRC grant MR/K004360/1 (Behavioral and Neurophysiological Effects of Schizophrenia Risk Genes: A Multi-locus, Pathway Based Approach). SD is also supported by a MARIE-CURIE COFUND EU-UK Research Fellowship. BR and the CUBRIC MEG lab are supported by an MRC UK MEG Partnership Grant, MR/K005464/1 and an MRC Doctoral Training Grant, MR/K501086/1. We would like to acknowledge RCUK of Cardiff University and Wellcome Trust for covering the publication fee.