Modeling the Switching behavior of Functional Connectivity Microstates (FCμstates) as a Novel Biomarker for Mild Cognitive Impairment

It is evident the need for designing and validating novel biomarkers for the detection of mild cognitive impairment (MCI). MCI patients have a high risk of developing Alzheimer’s disease (AD), and for that reason the introduction of novel and reliable biomarkers is of significant clinical importance. Motivated by recent findings about the rich information of dynamic functional connectivity graphs (DFCGs) about brain (dys)function, we introduced a novel approach of identifying MCI based on magnetoencephalographic (MEG) resting state recordings. The activity of different brain rhythms {δ, θ, α1, α2, β1, β2, γ1, γ2} was first beamformed with linear constrained minimum norm variance in the MEG data to determine ninety anatomical regions of interest (ROIs). A dynamic functional connectivity graph (DFCG) was then estimated using the imaginary part of phase lag value (iPLV) for both intra-frequency coupling (8) and also cross-frequency coupling pairs (28). We analyzed DFCG profiles of neuromagnetic resting state recordings of 18 Mild Cognitive Impairment (MCI) patients and 20 healthy controls. We followed our model of identifying the dominant intrinsic coupling mode (DICM) across MEG sources and temporal segments that further leads to the construction of an integrated DFCG (iDFCG). We then filtered statistically and topologically every snapshot of the iDFCG with data-driven approaches. Estimation of the normalized Laplacian transformation for every temporal segment of the iDFCG and the related eigenvalues created a 2D map based on the network metric time series of the eigenvalues (NMTSeigs). NMTSeigs preserves the non-stationarity of the fluctuated synchronizability of iDCFG for each subject. Employing the initial set of 20 healthy elders and 20 MCI patients, as training set, we built an overcomplete dictionary set of network microstates (nμstates). Afterward, we tested the whole procedure in an extra blind set of 20 subjects for external validation. We succeeded a high classification accuracy on the blind dataset (85 %) which further supports the proposed Markovian modeling of the evolution of brain states. The adaptation of appropriate neuroinformatic tools that combine advanced signal processing and network neuroscience tools could manipulate properly the non-stationarity of time-resolved FC patterns revealing a robust biomarker for MCI.


Abstract:
It is evident the need for designing and validating novel biomarkers for the detection of mild cognitive impairment (MCI). MCI patients have a high risk of developing Alzheimer's disease (AD), and for that reason the introduction of novel and reliable biomarkers is of significant clinical importance. Motivated by recent findings about the rich information of dynamic functional connectivity graphs (DFCGs) about brain (dys)function, we introduced a novel approach of identifying MCI based on magnetoencephalographic (MEG) resting state recordings.
The activity of different brain rhythms {δ, 2} was first beamformed with linear constrained minimum norm variance in the MEG data to determine ninety anatomical regions of interest (ROIs). A dynamic functional connectivity graph (DFCG) was then estimated using the imaginary part of phase lag value (iPLV) for both intra-frequency coupling (8) and also cross-frequency coupling pairs (28). We analyzed DFCG profiles of neuromagnetic resting state recordings of 18 Mild Cognitive Impairment (MCI) patients and 20 healthy controls. We followed our model of identifying the dominant intrinsic coupling mode (DICM) across MEG sources and temporal segments that further leads to the construction of an integrated DFCG (iDFCG). We then filtered statistically and topologically every snapshot of the iDFCG with data-driven approaches. Estimation of the normalized Laplacian transformation for every temporal segment of the iDFCG and the related eigenvalues created a 2D map based on the network metric time series of the eigenvalues (NMTSeigs). NMTSeigs preserves the non-stationarity of the fluctuated synchronizability of iDCFG for each subject. Employing the initial set of 20 healthy elders and 20 MCI patients, as training set, we built an overcomplete dictionary set of network microstates (nμstates).
Afterward, we tested the whole procedure in an extra blind set of 20 subjects for external validation.
We succeeded a high classification accuracy on the blind dataset (85 %) which further supports the proposed Markovian modeling of the evolution of brain states. The adaptation of appropriate neuroinformatic tools that combine advanced signal processing and network neuroscience tools could manipulate properly the non-stationarity of time-resolved FC patterns revealing a robust biomarker for MCI.

Introduction
The major cause of clinical dementia in the elderly is that of Alzheimer's type (DAT; Qiu et al., 2009), which is mainly characterized by loss of synapses, the accumulation of the Beta amyloid protein (Aβ) and the phosphorylation of the Tau protein. Due to the progressive loss of synapses, which alters the efficient communication within and between various brain subsystems, the DAT may be considered a disconnection syndrome (Delbeuck et al., 2003).
The pathological changes of DAT start decades before the first clinical symptoms appear, thus it is important to design proper analytic pathways for analysing neuroimaging datasets via, e.g., the notion of brain connectivity, which allows an early detecting of such changes (Gomez et al., 2009;Stam et al., 2009;Maestu et al., 2015). It is extremely important to Alzheimer's disease (AD) research to early identify preclinical and prodromal AD as it can assist clinical trials and targeted interventions (Livingston et al., 2017).

Mild cognitive impairment (MCI) is considered to be an intermediate clinical stage
between the normal cognitive decline and DAT (Petersen and Negash, 2008). The main parts of the affected brain during the MCI, apart from those involved in action and thought, are those related to memory. For that reason, MCI patients face memory problems on a higher level compared to normal aged population but with no prevalent characteristic symptomatology of dementia-like reasoning or impaired judgement (Petersen et al., 2009).
MCI is a heterogeneous state with different subtypes, which complicates in many cases the prediction to DAT (Portet et al., 2006). Additionally, it is also difficult to accurately discriminate symptomatic predementia (MCI) from healthy aging or dementia (DAT) (Petersen and Negash, 2008). Despite these difficulties to achieve an early diagnosis, an accurate identification of MCI should be attempted. Early diagnosis of MCI even in the absence of a healing strategy is significant for both pharmacological and nonpharmacological interventions. For that reason, new tools based on neuroimaging approaches are needed to increase the sensitivity in the detection of MCI.
Analysis of MEG recordings untangled the association between neural oscillations, functional connectivity assessment and neurophysiological activity (Brookes et al., 2011).
Altered frequency-dependent functional patterns have been linked to the progression of cognitive decline (Poza et al., 2014). Alternative scenarios of analysing MEG recordings include to single channel analysis e.g. power analysis, functional connectivity and brain network analysis in resting state and also in task-based experiments (for a review see Mandal et al., 2018). Analysis of single channel recordings is a less complex approach that identified aberrant oscillations in AD primarily in the left temporal-parietal-occipital brain areas (Gomez et al., 2009). Functional connectivity (FC) and effective (EC) connectivity analysis revealed a loss of connectivity in AD compared to healthy control (HC) subjects found mostly in higher frequency bands (Gomez et al., 2017) while multiplex network analysis of MEG study in AD identified affected regions of hippocampus, posterior default mode network (DMN) and occipital areas (Yu et al., 2017). However, current clinical literature is limited and no strong conclusion can be drawn.
A recent multicenter MEG study addressed this issue using functional connectivity (FC) analysis (Maestu et al., 2015). It revealed a hypersynchronization in MCI as the most discriminative feature of brain connectivity, mainly over the fronto-parietal and interhemispheric links. This pattern was stable across the five different neuroimaging centres that participated in the study (Accuracy ~= 80 %), which might thereby be considered as a preclinical connectomic biomarker for MCI/DAT. Previous MEG studies based on connectivity analysis described a less organized functional brain network a hypersynchrony in fronto-parietal network in MCI subjects (Bajo et al., 2010;Buldú et al., 2011), while patients with DAT demonstrated a less synchronized brain network accompanied with a cognitive decline ). This hypersynchronization might be a compensatory mechanism but it cannot be adaptive since patient's network is closes compared to healthy olds to random (Buldú et al., 2011). In a recent MEG study comparing progressive MCI and stable MCI, authors described hypersynchronization in α band between anterior cingulate and posterior brain areas in the progressive MCI group (López et al., 2014.) Spontaneous fluctuations of functional MRI (fMRI) blood-oxygen-level-dependent (BOLD) signals are temporally coherent between distinct spatial brain areas and not random. Biswal et al., (1995) demonstrated that fluctuations from motor areas were correlated even in the absence of a motor task. FC based on BOLD signal is modulated by cognitive and affective states (Richiardi et al., 2011 andShirer et al., 2012), by learning (Bassett et al., 2011), and also spontaneously , Chang and Glover, 2010and Kitzbichler et al., 2009.
When non-stationarity is taken into account and a dynamic functional connectivity (DFC) approach is adopted for studying FC patterns even in the absence of a task (restingstate), more sophisticated algorithmic analyses should be used. In this line, two studies have been recently published simultaneously that presented a data-driven methodology. In the first one, Allen et al. (2014) proposed a method based on k-means clustering, aimed to detect distinct 'FC states' in the resting brain. These authors clearly showed differences from the stationary static functional brain networks. The second study proposed a data-driven method focused on extracting, out of hundreds of functional connectivity graphs (FCGs) in a multitrial experimental paradigm, distinct brain states called functional connectivity microstates (FCμstates; Dimitriadis et al., 2013a). Both approaches revealed the need of dynamic FC to explore brain dynamics via the notion of brain connectivity, as it is clear that brain FC 'hops' from one state to another (FCμstate) leading to a Markovian chain with characteristic favoured transitions between pairs of FCμstates Dimitriadis et al., 2010aDimitriadis et al., , 2013a.
In parallel, the concept of cross-frequency coupling (CFC) is gaining attention lately in the neuroscience community, as evinced by the increasing number of papers published with the incorporation of this type of interaction in the analysis (Florin and Baillet, 2015;Tewarie et al., 2016;Antonakakis et al., 2016a,b;Dimitriadis et al., 2015cDimitriadis et al., ,d, 2016avan Wijk and Fitzgerald, 2014). Specifically, intrinsic coupling modes and especially CFC bias the task-related response and are sensitive to various brain diseases and disorders such as shown that the dynamics of spontaneously generated neural activity can be informative about the functional organization of large-scale brain networks (Fox et al., 2005;He et al., 2008;Hipp et al., 2012), revealing intrinsically generated "coupling modes" at multiple spatial and temporal scales (Deco and Corbetta, 2011;Engel et al., 2013).
Based on the aforementioned methodological evidences in microscale, it is significant to explore the repertoire of intra and cross frequency interactions across brain rhythms and brain areas under the same integrated graph model (Antonakakis et al., 2017;Dimitriadis et al., , 2017bDimitriadis et al., , 2018a. In a previous study, we demonstrated how to design a connectomic biomarker for MCI based on source-reconstructed MEG activity via static brain network analysis (Dimitriadis et al., 2018a). Here, we extended this work by proposing a scheme to design a dynamic connectomic biomarker under the framework of dynamic functional connectivity (DFC) analysis. Additionally, the proposed scheme will be validated in a second blind dataset (SID did not know anything about the labels).
For this aim, we analysed the MEG activity of healthy controls and MCI patients at resting-state (eyes-open) via DFC analysis. Based on a previous approach Dimitriadis et al., 2017b), we detected the dominant type of interaction per pair of MEG sources and temporal segment (Dimitriadis et al., 2018c). This approach produced a subject-specific dynamic functional connectivity graph (DFCG). This approach created a 2D matrix of size sources x temporal segments that describe the evolution of the eigenvalues across experimental time. Afterward, we used neural-gas to design overcomplete dictionaries for sparse representation of NMTS eigen independently for the two groups (Dimitriadis et al., 2013a). Then, we validated the whole approach in a blind dataset to quantify the generalization of the proposed method.
In Materials and Methods section, we described the data acquisition, preprocessing steps, information about the datasets and the proposed methodological scheme. The Results section is devoted to describe the results including the prototypical network FCμstates, the accuracy of prediction in a blind dataset and network-based information of brain states.
Finally, the Discussion section includes the discussion of the current research results with future extensions.

Subjects and and Ethics Statement
The training dataset includes eighteen right-handed individuals with MCI (71.89 ± 4.51 years of age) and twenty-two age-and gender-matched neurologically intact controls (70.91 ± 3.85 years of age) were recorded. Table 1 summarizes their demographic characteristics. All participants were recruited from the Neurological Unit of the "The Hospital Universitario San Carlos", Madrid, Spain. They were right handed (Olfield, 1971), and native Spanish speakers. We used also a set of twenty subjects of unknown label (blind author SD) for further validation of the proposed dynamic connectomic biomarker (DCB). for blind classification.
To explore their cognitive and functional status, all participants were screened by means of a variety of standardized diagnostic instruments and underwent an extensive cognitive assessment, as described in López et al. (2016).
MCI diagnosis was established according to the National Institute on Aging-Alzheimer Association (NIA-AA) criteria (Albert et al., 2011), being all of them categorized as "MCI due to AD intermediate likelihood". They met the clinical criteria and also presented hippocampal atrophy, which was measured by magnetic resonance (MRI). According to their cognitive profile, they were also classified as amnestic subtype (Petersen et al., 2001).
The whole sample was free of significant medical, neurological and/or psychiatric diseases (other than MCI). Exclusion criteria included: a modified Hachinski Ischemic score ≥ 4 (Rosen et al., 1980); a geriatric depression scale short-form score participants with medical treatment that could affect the MEG activity (e.g., cholinesterase inhibitors) were required to interrupt it 48 h before the MEG recordings.
The present study was approved by the local ethics committee and all subjects signed an informed consent prior to their MEG recording.
[  (Fischl et al., 2002)) was used to obtain the hippocampal volumes, which were normalized with the overall intracranial volume (ICV) of each subject.

MEG Acquisition and preprocessing
Four minutes of eyes-open resting state data were recorded while the participants were seated in a 306-channel (one magnetometer and two orthogonal planar gradiometers per recording site, sampling frequency of 1kHz ) Vectorview system (ElektaNeuromag) placed in a magnetically shielded room (VacuumSchmelze GmbH, Hanau, Germany) at the "Laboratory of Cognitive and Computational Neuroscience" (Madrid, Spain). Subjects had to fix their look in a cross, which was projected in a screen. The position of the head relative to the sensor array was monitored by four head position indicator (HPI) coils attached to the scalp (two on the mastoids and two on the forehead). These four coils along with the head shape of each subject (referenced to three anatomical fiducials: nasion, and left-right preauricular points) were acquired by using a three-dimensional Fastrak Polhemus system.
Vertical ocular movements were monitored by two bipolar electrodes, which were placed above and below the left eye, and a third one to the earlobe, for electrical grounding.
Four head position indicator (HPI) coils were placed in the head of the subject, two in the forehead and two in the mastoids, in order to online estimate the head position. The HPI coils were fed during the whole acquisition, allowing for offline estimation of the head position.
Maxfilter software (version 2.2 Elekta Neuromag) was used to remove the external noise from the MEG data using the temporal extension of signal space separation (tSSS) with movement compensation (correlation threshold= 0,9m time window= 10s) (Taulu and Simola, 2006). This algorithm removes the signals whose origin is estimated outside the MEG helmet, while keeping intact the signals coming from inside the head. In addition, the continuous HPI acquisition, combined with the tSSS algorithm, allowed for the continuous movement compensation. As result, the signal used in the next steps comes from a set of virtual sensors whose position remains static respect to the head of the subject. Those subjects whose movement along the recording was larger than 25 mm were discarded, following the recommendations of the manufacturer.

Source Reconstruction
We generated a volumetric grid for the MNI template by adopting a homogenous separation of 1 cm in each direction, with one source placed in (0, 0, 0) in MNI coordinates.
The whole procedure resulted in a source model with 2,459 sources inside the brain surface where each one consisted of three perpendicular dipoles. Every source was then labelled using the automated anatomical labelling (AAL) atlas (Tzourio-Mazoyer et al., 2002). We finally considered 1,467 cortical sources. The computed grid was then transformed to subject specific space employing the original T1 image. The realignment of grid and brain surface was realized manually to Neuromag coordinate system following the three fiducials and the head shape guides. Employing a realistically shaped head, we estimated a lead field (Nolte, 2003). We source reconstructed frequency-dependent brain activity using a Linearly Constrained Minimum Variance (LCMV) beamformer (Van Veen et al., 1997). We ran LCMV beamformer independently for the following eight frequency bands: γ 2 (51-90 Hz). The resulting spatial filters were projected over the maximal radial direction, getting only one spatial filter per source. "Radial direction" means the direction of the segment connecting the dipole location to the center of the sphere best approximating the brain surface. Radial dipoles in a spherical conductor do not produce a magnetic field outside of the conductor (Sarvas, 1987), so this projection avoids the creation of undetectable sources among the target dipoles. Finally, we represent every brain area -region of interest according to the AAL atlas by one source-space time series per frequency band using two alternative solutions: (1) the PCA of all the sources in the area or (2) the source closest to the centroid of the area (CENT).  [ Figure 1 around here]

Construction of the integrated DFCGs
The DCFG analysis was restricted to the 90 ROIs of the AAL atlas. Adopting a common sliding window of width equal to 1s, to get at least 1 cycle of δ activity and a moving step of 50ms, we estimated the dynamic networks for both intra-frequency (8 frequency bands) and inter-frequency coupling modes (8*7/2=28 cross-frequency pairs) using the following formula of the imaginary part of phase locking value (iPLV).
where φ(t) is the phase of the signal in the corresponding frequency band (intra-frequency modes) and between frequencies (cross-frequency couplings). For further details regarding phase-to-amplitude cross-frequency coupling (see Dimitriadis et al., 2015d and section1 in sup. material).
This procedure, whose implementation details can be found elsewhere (Dimitriadis et al., 2010a(Dimitriadis et al., , 2015d(Dimitriadis et al., ,2017a(Dimitriadis et al., ,2017b(Dimitriadis et al., , 2018b, resulted in a four-dimensional tensor of size  A. Schematic illustration of our approach employed to identify the DICM between two sοurces (Left superior frontal gyrus, Right superior frontal gyrus) for two consecutive sliding time windows (t 1 , t 2 ) during the first 4 secs of resting-state activity from the first healthy control subject. In this example, the functional synchronization between band-passed signals from the two sources was estimated by imaginary Phase Locking (iPLV). In this manner iPLV was computed between the two sources either for samefrequency oscillations (e.g., δ to δ …, γ 2 -γ 2 ; 8 intra-frequency couplings) or between different frequencies (e.g., δ to θ , δ to α 1 ..., γ 1 -γ 2 ; 28 cross-frequency pairs). The sum of 8 + 28 = 36 refers to Potential Intrinsic Coupling Modes (PICM), which are tabulated in a matrix format. In the main diagonal, we inserted the intra-frequency couplings while in the off-diagonal the cross-frequency pairs. Statistical filtering, using surrogate data for reference, was employed to assess whether each iPLV value was significantly different from chance. During t 1 the DICM reflected significant phase locking between α 1 and α 2 oscillations (indicated by red rectangles) in the oscillation list and a '*' in the comodulogram. The DICM remains stable also for the t 2 between α 1 and α 2 oscillations whereas during t 3 the dominant interaction was detected between θ and α 2 oscillations. B. Burst of DICM between Left and Right superior frontal gyrus.
This packeting can be thought to associate the 'letters' contained in the DICM series to form a neural "word.", a possible integration of many DICM. From this burst of DICM, we can quantify the probability distribution (PD) of DICM across experimental time (see C). C. Tabular

Statistical Filtering Scheme
At first place, we have to identify true CFC interactions that are not driven by the changes in signal power. Secondly, following a proper surrogate analysis our DICM model can detect the DICM between every pair of sources and at every temporal segment. The whole procedure of analysis is described in detail in Dimitriadis and Salis, 207b;Dimitriadis, 2018c). and also in see section 2 in sup.material

Topological Filtering Scheme based on OMSTs
Apart from surrogate analysis, which is a statistical filtering procedure of the functional couplings within an FCG akin to a regularization to sparsify the 4D array described above, we adopted a topological filtering to further enhance the network topology and the most significant interactions . To this aim, we applied to each FCG derived from each subject and temporal segment independently, a novel data-driven thresholding scheme, proposed by our group, and termed Orthogonal Minimal Spanning Trees (OMSTs; Dimitriadis et al., 2017aDimitriadis et al., , 2017c. Fig.3.B demonstrates the temporal evolution of the topologically filtered dFCG for the first ten temporal segments.

Graph Signal Processing
After extracting the most significant connections in DCFGs from each individual, we transformed every snapshot of the DFCG into the graph Laplacian variant called the normalized Laplacian matrix. With A being the functional connectivity graph and D, the degree matrix containing the degree of every node in the main diagonal, graph Laplacian L can be defined as L= D -A. The normalized graph Laplacian is defined as L sym = D −1/2 LD −1/2 (Shuman et al., 2013). We estimated the sorted eigenvalues of the Lsym for every snapshot of DFCG resulting to a 2-dimensional matrix of size [source (90) x temporal segments (2.401)] per subject. These 2D matrices were concatenated separately for the healthy control and disease group of the training set. Practically, the concatenation was performance along the temporal direction. Fig.3.C shows the temporal evolution of the normalized Laplacian transformation of the dFCG for the first ten temporal segments while Fig.3.D is dedicated to the temporal evolution of the eigenvalues.

A Vector-Quantization (VQ) Modeling of Group NMTS eigen
This subsection describes briefly our symbolization scheme, presented in greater details elsewhere (Dimitriadis et al., 2011(Dimitriadis et al., , 2012(Dimitriadis et al., , 2013a. The group-specific NMTS eigen patterns can be modelled as prototypical FC microstates (FCμstates). In our previous studies, we demonstrated a better modelling of DFCG based on vector quantization approach (Dimitriadis et al., 2013a(Dimitriadis et al., , 2017b(Dimitriadis et al., , 2018b. A codebook of k prototypical FC states (i.e., functional connectivity microstates-FCμstates) was first designed by applying the neural-gas algorithm (Dimitriadis et al., 2013a). 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 soft-max adaptation rule that minimizes the average distortion error (Martinetz et al., 1993). Neural-gas algorithm has been applied independently to each group by concatenating the 2D matrix of size [2.401 x 90] that describes the fluctuation of Laplacian eigenvalues.
The outcome of the neural-gas algorithm over NMTS eigen is the construction of a symbolic sequence of group-specific prototypical FCμstates, one per subject. An example of such symbolic time series (STS) is a Markovian chain with three FCμstates: {1, 2, 3, 2, 1, 3, 2…} where each integer defines a unique brain state (FCμstates) assigned to every quasi-static temporal segment.

External Validation in a blind dataset
We designed a novel approach of how to classify a blind subject. We reconstructed the subject-specific NMTS eigen with both HC-based prototypical FCμstates and MCI-based prototypical FCμstates Specifically, for every temporal segment expressed via a vector of 90 eigenvalues, we estimated which of the prototypical FCμstates is much closer employing Euclidean distance for an appropriate criterion. Under this scheme, we rebuilt the original NMTS eigen twice, one using prototypical FCμstates of HC and one using prototypical FCμstates of MCI. Then, we estimated the reconstruction mean squared error between the original NMTS eigen and the two rebuilt NMTS eigen based on prototypical FCμstates .Finally, we assign the test sample to the class with the lowest reconstruction error.

Markov chain modelling for synchro state transitions
The temporal sequence of spontaneous activity can be modelled as a Markovian process, which predicts the probabilities of several discrete states recurring or switching among themselves at different time points analysing time-point-based brain activity (Gärtner et al., 2015;. Several studies have investigated transition probabilities between phase-synchronized states on a sub-second temporal scale, untangling the Markovian property and the switching behaviour of finite network-level brain states (Baker et al., 2014;Dimitriadis et al., 2013Dimitriadis et al., ,2015Jamal et al., 2015).

Markovian process of time-sequential FCμstates
Markov model describes the underlying dynamical nature of a system that follows a chain of linked states, where the appearance of a state at any given time instant depends only on the preceding ones (Gagniuc, 2017). In the Markov chain modelling for synchrostate transitions during the deductive reasoning and task-free processes, the first order transition matrices were estimated in a probabilistic framework. According to discrete-time Markov chain theory (Jarvis and Shier, 1999), a finite number (S 1 , S 2 …, S m ) of inferred states, that evolve in discrete time with a time-homogeneous transition structure can be mathematically represented by either its -by-transition probability matrix or its directed graph (digraph).
Here, the inferred states refer to the prototypical FCμstates. A feasible transition is one whose occurring probability is greater than zero. The probability of transition from node (state) i to node j is defined as where Nij is the number of transitions from node i to node j. Obviously, the sum of the transition probabilities along each row of the transition matrix P equals one. The complete digraph for a finite-state Markov process has edges of transition probabilities between every node i and every other node j. Here, nodes refer to FCμstates in the Markov chain. In the digraphs created in this study, if P ij survives a p-value derived by 10,000 shuffled-surrogates of the original symbolic time series.

Assessing the Statistical Significant Level of the Symbolic-based Estimates
To assess the statistical significant level of the aforementioned four estimates (excluding CI), we shuffled the group symbolic time series 10,000 times and we re-estimated the surrogate-based p-values for every estimate per subject. CI is normalized by default with surrogates.

Linking MMSE with Chronnectomics
To investigate the possible relation between MMSE and the chronnectomics derived by the FCμstate symbolic sequence (see section 2.6.2), we used the canonical correlation analysis (CCA) approach to see whether MMSE correlates with seven chronnectomic variables.

Algorithms and MATLAB Code
All the algorithmic steps of constructing the DCFGs were implemented on in house software written in MATLAB, freely available from the first author's github site https://github.com/stdimitr/from_dfcg_to_brain_states.
Same analysis is supported in PYTHON in a novel toolbox that van be found in the following link: https://github.com/makism/dyfunconn LCMV beamformer has been programmed under Fieldtrip's environment (Oostenveld et al., 2002). Fig.4 illustrates the prototypical group-specific FCμstates for each group by assigning the 90 AAL brain areas to five well-known brain networks. The size and colour of every circle decodes the mean degree within every brain network while the colour of each connection defines the mean functional strength between every pair of brain networks. FCμstates can be described based on the most connected brain networks focusing on their degree.

Netμstates eigen
Each test sample with unknown label was classified to one of the two classes using as criterion the minimization of the reconstruction error. The minimum reconstruction error denotes the class label of the sample. In our study, we used twenty samples with distribution of 11 MCIs and 9 controls with 85% accuracy for CENT (17 out of 20) and 70% for PCA representation scheme (14 out of 20). SID received the blind dataset from MEL who evaluated the outcome of this research. Fig.6 illustrates the methodological approach. Fig.6A refers to the temporal resolution of the Laplacian eigenvalues of a blind HC subject while  [ Figure 6 around here]

Group-differences of Temporal measurements derived from FCμstate symbolic sequence
FI, OT and DT were significantly higher than the surrogates based FI values derived from the shuffled symbolic time series (p < 0.001). We detected significant higher FI and CI for HC compared to MCI applying a Wilcoxon Rank-Sum test (Fig.7, p-

Modelling Dynamic Reconfiguration of Functional Connectivity Graphs as a Markovian Chain.
The outcome of the VQ modelling of NMTS eigen is the derived Netμstates eigen called FCμstates (see Fig.4) where its evolution is described via a symbolic time series, a Markovian chain. Fig. 9 illustrates a well-known scheme of the group-averaged transition probabilities (TP) between the three FCμstates for both groups. Our analysis revealed significant group-differences in terms of TP, while the TP were significant different compared to the surrogates' symbolic time series. Self-loops defined the 'staying' TP of brain dynamics to the same brain state.

Figure 9. A finite-state diagram showing group-averaged transition probability matrix
(TP) of the symbolic time series, which describes the temporal evolution of the brain,

states (FCμstate). A) For healthy control and B) for mild cognitive impairment.
The symbolic time series illustrated in Fig. 5 is a Markovian chain of order 1 and it is shown schematically with a diagram of three nodes defining the three FCμstates (Fig.4) while the arrows from one state to the other show the TP. Our results revealed significant group differences between every possible brain state transition (Wilcoxon Rank-Sum test, p < 0.0001/9).

Comodulograms of Dominant Intrinsic Coupling Modes (DICM)
Probability distributions (PD) of prominent intrinsic coupling modes across all sensors, sensor pairs, and time windows were summarized for each group in the form of an 8 x 8 matrix. The horizontal axis refers to the phase modulating frequency (Hz) while the y-axis to the amplitude modulated frequency (Hz). The main diagonal of the comodulograms keep the PD of intra-frequency phase-to-phase coupling. Group-averaged comodulograms in Figure 9 demonstrate the empirical PD of DICM revealing a significant role of α 1 as phase modulator of the whole studying spectrum up to high-gamma (γ 2 ) activity, which covers almost the 50% of pair-wise sensor connections and time-windows. No significant trend was detected regarding the PD of each pair of frequencies between the two groups (p < 0.05, Wilcoxon rank-sum test, Bonferroni corrected). Moreover, no significant difference was found regarding the PD of the groups for every possible pair of sensors (p < 0.05, Wilcoxon ranksum test, Bonferroni corrected). Finally, transition dynamics of DICM between consecutive time-windows at every sensor-pair did not uncover any group-difference (for further details see Dimitriadis et al., 2016b).

Figure 10. Group-averaged empirical Probability Distribution values of dominant intrinsic coupling modes for MCI (B) compared to control group (A).
[ Figure 10 around here]

Discussion
We have demonstrated here a novel framework for designing a proper dynamic connectomic biomarker (DCB) for the detection of MCI subjects from spontaneous neuromagnetic activity. The whole approach exhibits novel, data-driven, algorithmic steps that can be summarized as follows: • The construction of a dynamic functional connectivity graph (DFCG) that incorporates dominant types of interactions, either intra-(e.g. θ -θ) or inter-frequency (phase-to-amplitude coupling (PAC) (e.g. θ -γ)) coupling.
• The application of a new thresholding scheme termed Orthogonal Minimal Spanning Trees (OMSTs) as a topological filtering applied to DFCG to extract a 'true' network topology.
• The VQ modelling of network metric time series (NMTS) based on nodal Laplacian eigenvalues for prototyping the spatiotemporal dynamics of both control and MCI subjects.
• Modelling of the switching behavior of brain states as a Markovian chain • The validation of the whole approach to a second blind dataset achieving 85% of classification accuracy for the CENT ROI representation scheme compared to 70% for the PCA scheme • ROI representation scheme matters on the designing of connectomic biomarker in general and also for MCI • CCA analysis between chronnectomics and MMSE revealed that the DT of brain states associates strongly with the CCA mode of MMSE variability.
We proved that the VQ modelling of NMTS eigen is an effective approach to extract an overcomplete dictionary for the representation of dynamic functional connectivity that can accurately classify subjects as either control or MCI based on their resting state MEG activity. Adopting a static network analysis, the classification accuracy was 12 out of 20, demonstrating the need of a dynamic functional connectivity approach for studying resting brain dynamics Damarazu et al., 2014;Dimitriadis et al., 2010aDimitriadis et al., , 2012cDimitriadis et al., ,d, 2013aDimitriadis et al., ,b, 2015bDimitriadis et al., ,d, 2016bDimitriadis et al., , 2017aDimitriadis et al., ,b, 2018aKopell et al., 2014).
The capture of time-varying coupling between variables is a topic that has been heavily studied in other fields and in communications for signal processing in particular.
However, the specific application to whole-brain functional connectivity is relatively new (Sakoglou et al., 2010;Dimitriadis et al., 2013a;Calhoun et al., 2014), and its application to brain-imaging data poses particular challenges, which are the topic of active current research.
One important challenge is how to best identify relevant features from the high-dimensional brain imaging data. The main algorithms highly used for manipulating functional brain network dynamics in fMRI are group ICA (Calhoun and Adali, 2012) or spatial-constrained ICA (Lin et al., 2010) and tensor decompositions (Acar and Yener, 2009). Complementary, to characterize the dynamics of time-varying connectivity brain patterns, the basic approach is the metastate analysis based on the sliding window or more adaptive approach (Damaraju et al., 2015;Dimitriadis et al., 2013a;2015a;Nomi et al., 2016). From the dynamic connectivity patterns, functional connectivity microstates (FCμstates) are extracted that are 'quasi-stable' distinct brain states. Then, the state vectors can be modelled via a Markovian chain (Calhoun et al., 2015;Damaraju et al., 2015;Dimitriadis et al., 2013a;. CFC mechanisms support the brain interactions across space over multiple temporal scales (Canolty & Knight, 2010;Fell & Axmacher, 2011). Computational models have explored the theoretical advantages of the existence of cross frequency coupling (Lisman & Idiart, 1995;Neymotin et al., 2011). These models untangled the major mechanisms of the importance of CFC, which may serve as brain's neural syntax. Segmentation of spike trains into cell clusters ('letters') and sequences of assemblies (neural 'words') is supported by the existing syntactic rules (Buzsaki, 2010).
In the present study, we demonstrated a methodology whose main scope is to provide a framework for modelling DFCG into a repertoire of distinct 'quasi-static' brain states called FCμstates. Here, we modelled the network metric time series (NMTS) derived by the dynamic functional connectivity patterns expressed via the Laplacian eigenvalues (Dimitriadis et al., 2015a(Dimitriadis et al., , 2017b(Dimitriadis et al., , 2018b. After extracting the virtual source time series, we followed an algorithmic approach with main aim to minimize the effect of a priori selection of variables that can minimize the reproducibility of the results. The main steps of the proposed methodology are: 1) the construction of one integrated dynamic functional connectivity per subject. which incorporates the dominant intrinsic coupling mode per each pair of brain areas and at every temporal segment, 2) the application of a data-driven topological filtering scheme to reveal the backbone of the network topology at every temporal segment, 3) the estimation of Laplacian eigenvalues to extract the so-called network metric time series (NMTS eigen ; Dimitriadis et al., 2010Dimitriadis et al., ,2015aDimitriadis et al., ,2017bDimitriadis et al., ,2018b, 4) the modelling of these NMTS eigen via a vector-quantization approach and 5) the validation of the whole approach to a second blind dataset.
The analysis of spatio-temporal evolution of Laplacian eigenvalues during the training phase revealed three prototypical brain states (FCμstates). For a better illustration of the FCGs linked to the prototypical eigenvalues, we assigned the 90 AAL brain areas to five well-known brain networks. In Fig. 4, we mapped the average functional strength between ROIs belonging to every pair of brain networks while the size and colour of every node defines the within brain network degree. The most connected brain networks in FCμstates are the DMN and CO. CO plays a key role in working memory mechanisms (Wallis et al., 2015) while cognitive complaints related to AD are linked to alterations of resting-state brain networks and mostly FPN and DMN (Contreras et al., 2017). Functional coupling strength between FPN and DMN was significantly higher for HC compared to MCI for FCμstates 1 and 3 (Fig. 4). The functional strength between FPN-DMN was positively correlated with a better episodic memory performance (Contreras et al., 2017).
Well-known and novel chronnectomics were estimated from the Markovian (symbolic) Chain that describes the evolution of brain states. We detected significant higher flexibility and complexity for HC as compared to MCI described from FI and CI, correspondingly ( Fig.7). A summarization derived from OT and DT revealed a significant trend: HC subjects spent significantly higher time compared to MCI to FCμstates 1 and 3 while MCI spent significantly more time to the second FCμstate (Fig. 8). Following a CCA analysis between the extracted chronnectomics and the MMSE score; we found a significant contribution of the DT for the three NMTS eigen . OC related to the 2 nd NMTS eigen didn't associate with the CCA mode of MMSE variability (Fig.11).
In the era of data-sharing and aggregating large datasets from different research groups worldwide contributed to big consortiums, it is important to test the reproducibility of the proposed biomarkers (Abraham et al., 2016). Our study is a first step in this direction to diminish the effect of any arbitrary selection of algorithmic steps till the extraction of biomarkers. The next step is to reproduce the results on the source space and to extend the analysis in larger populations from different sites and MEG scanners. A recent study showed that a 70 percent of the preclinical research from academic labs could not be replicated (Collins and Tabak, 2014). Abraham's work is one of the very first neuroimaging studies that lays the ground on the reliability and reproducibility of biomarkers extracted from neuroimaging data.
There is a large body of research based on different imaging methods covering various temporal and spatial scales that document the association of electrophysiological rhythms with distinct cognitive processes within narrowly or broadly anatomical areas (for review see Engel et al., 2001;Buszaki, 2006;Siegel et al., 2012;Bazar et al., 2013). For example, low frequency δ rhythms (1-4 Hz) are known to coordinate large portions of the brain (Fujisawa and Buzsaki, 2011;Nacher et al., 2013) while γ oscillations play a dominant role in stimulus processing and detection locally anatomically constrained (Engel et al., 2001). Recently, an extension of Brodmann's areas was suggested in order to associate distinct anatomical areas with preferable connectivity estimator and cognitive functions in both normal and brain disease/disorder populations as an initial step to summarize the large body of current brain connectivity research (Basar et al., 2016).
In the last few years, an increasing number of studies appeared studying crossfrequency coupling (CFC) at resting-state (Antonakakis et al., 2016a), during cognitive tasks (Dimitriadis et al., 2015c(Dimitriadis et al., ,d, 2016a and in various brain diseases and disorders such as mild traumatic brain injury (Antonakakis et al., 2016a), amnestic MCI (Dimitriadis et al., 2015d), dyslexia , schizophrenia (Kirihara et al., 2012), etc. It has been suggested that CFC is the key mechanism for the integration of information between anatomically distribution subsystems that function on a dominant frequency (Canolty and Knight, 2010;Jirsa and Muller, 2013;Florin and Baillet, 2015). However, only a few MEG studies have explored CFC at resting-state (Antonakakis et al., 2015a(Antonakakis et al., , 2016Florin and Baillet, 2015) and especially in a more dynamic fashion (Antonakakis et al., 2016b;Dimitriadis et al., 2015dDimitriadis et al., , 2016b. Our next goal is to work on the MEG source level after first improving sensitivity of beamformers in amplitude, phase and specifically to phase-to-amplitude coupling (PAC).
MEG source connectivity is in a mature level compared to a decade ago (Ioannides et al., 2012), and it is an active research area aimed at improving many aspects of a 'true' brain connectivity (Schoffelen and Gross, 2009;Colclough et al., 2016). The most significant issue is the parcellation of cerebral cortex. In many cases, the AAL template (90 ROIs) is feasible for the detection of those changes induced by a specific task or obtained after compare different groups. But in others such as the design of a reliable connectomic biomarker, there is a need to oversample to more than 90 areas. The FC which is directly linked to functional parcellation of cerebral cortex is an active area, which will further improve both the interpretation and the predictive power of source connectivity of many brain diseases such as MCI. The solution of a functional parcellation template for MEG source connectivity will improve the classification performance on the source level with the additional advantage, compared to sensor level, of facilitating the anatomical interpretation of the results.
Adopting the same framework and including also stable and progressive MCI groups, we will attempt to connect DCB with neuropsychological measures and cognitive scores (Cuesta et al., 2015a,b). It is evident that a multifactorial model that includes cognition, neuropsychological measures and anatomical information can reliably predict the conversion from MCI to DAT, while genetic variation of risk genes like APOE-e4 allele or cognitive reserve might play a secondary role (Lopez et al., 2016).
Going one step further from our previous studies demonstrating the significance of a DCB (Dimitriadis et al., 2013b(Dimitriadis et al., , 2015b, where we used network microstates extracted from DFCG patterns, in the present study we introduced a sparse modelling approach of NMTS eigen estimated over DFCGs that preserve the dominant type of coupling (intra or inter-frequency intrinsic coupling mode). Our study demonstrates the effectiveness of the data-driven analytic pipeline tailored to DFCG to the correct classification of a blind dataset based on control and MCI subjects compared to a static connectivity approach. Given these outcomes, it is evident the need, in the next years, to adopt data-driven techniques that will not introduce bias, subjectiveness and assumptions in neuroimaging datasets and also to improve the reproducibility of the outcome in large databases.
In magnetoencephalography (MEG) the conventional approach to source reconstruction is to solve the underdetermined inverse problem independently over time and space. Different algorithms have been proposed so far with alternative regularization procedures of the space and time like with a Gaussian random field model (Solin et al., 2016).
Commonly used techniques include the minimum-norm estimate (MNE) (Hamalainen & Ilmoniemi, 1994) and Linearly Constrained Minimum-Variance (LCMV) beamformer (Van Veen et al., 1997). It is on the right direction to compare the consistency of the outcome of current study with alternative inverse solution algorithms to measure their consistency and sensitivity to the design of connectomic biomarkers tailored to MCI.

Conclusions
In this study, we presented a novel DCM for the prediction of MCI from age-matched control group validated over a blind dataset. The novelties of the proposed analytic scheme is the incorporation in the DFCGs the dominant intrinsic coupling mode (DICM, either intra or inter-frequency coupling based on PAC), the adaptation of a novel data-driven thresholding scheme based on OMSTs, the estimation of laplacian eigenvalues across time and the extraction of prototypical network microstates (FCμstates) for both control and MCI group.
It is important for the near future to work in source space on MCI subjects that convert to AD after a following up study to further validate the proposed scheme as a potential tool of clinical importance. Complementary, it would be interesting to explore how Apoe-e4 allele can induce changes to the DFC of spontaneous activity. Moreover, multimodal neuroimaging biomarkers is a novel trend that will further been validated (Jack et al., 2016).

Author Contributions
SID: conception of the research analysis, methods and design, data analysis, and drafting the manuscript; MEL: data acquisition; MEL, FM, EP: critical revision of the manuscript. Every author read and approved the final version of the manuscript.

Conflict of Interest Statement
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.  Table 2. Mean ± standard deviation of the demographic characteristics of the blind sample of controls and MCIs. M= males; F= females; Educational level was grouped into five levels: 1) Illiterate, 2) Primary studies, 3) Elemental studies, 4) High school studies, and 5) University studies; MMSE= Mini-Mental State Examination; LH ICV, Left hippocampus normalized by total intracranial volume (ICV); RH ICV, Right hippocampus normalized by ICV