Scalp EEG signals were gathered from the free online database PhysioNet BCI (Database physionet BCI¹) (Schalk et al., 2004). The database consists of N = 94 healthy subjects recorded in two different baseline conditions, i.e., 1-min eyes-open (EO) resting state and 1-min eyes-closed (EC) resting state. In each condition, subjects were comfortably seated on a reclining chair in a dimly lit room. During EO they were asked to avoid ocular blinks in order to reduce signal contamination. The EEG data were recorded with a 64-channel system (BCI2000 system (BCI2000 system) with an original sampling rate of 160 Hz. From the 94 healthy subjects 39 were male with mean_age = 37.69 και std_age = 12.63 and the rest 55 were female with mean_age = 39.4 και std_age = 10.22. The range of age in the training sample was between 19 and 67.

We used a second EEG dataset for external validation of the proposed scheme based on 100 subjects from ages 18—60 from the Human Brain Institute (HBI) Normative Database using 19 EEG channels and the first min from the three recording min of eyes-open and eyes-closed condition.

Ongoing activity was corrected for artifacts through a two-step procedure implemented in Matlab (The MathWorks, Inc., Natick, MA, USA) and Fieldtrip (Oostenveld et al., 2011). Line noise was first removed using a notch filter at 60 Hz and the single-subject data was whitened and reduced in dimensionality by means of Principal Component Analysis (PCA) with a threshold corresponding to 95% of total variance (Delorme and Makeig, 2004; Antonakakis et al., 2013). The resulting signals were submitted to ICA using the extended Infomax algorithm as implemented in EEGLAB (Delorme and Makeig, 2004). A given independent component was considered to reflect ocular or cardiac artifacts if more than 30% of its z-score kurtosis or skewness values, respectively, were outside ±2 of the distribution mean (Antonakakis et al., 2013; Dimitriadis et al., 2015a). The remaining ICs were used to reconstruct a relatively artifact-free signal. The average number of artefactual ICs was 5.5 for eyes-closed and 5.1 for the eyes-open condition (see Figure 1A).

Here, functional connectivity was examined among the following eight brain rhythms of the typical sub-bands of electrophysiological neural signals {δ, θ, α₁, α₂, β₁, β₂, γ1, γ2}, defined respectively within the ranges {0.5–4; 4–8; 8–10; 10–13; 13–20; 20–30; 30–48; 52–70 Hz}. We adopted a 3rd order Butterworth filters applied in a zero-phase mode to get the characteristic brain rhythms (see Figure 1A).

The goal of the analytic procedures described in this section is to understand the repertoire of phase-to-phase interactions and their temporal evolution, while taking into account the quasi-instantaneous 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 1-min continuous EEG recording for both eyes-open and closed condition. The width of the temporal window was set equal to the duration of 10 cycles of δ brain rhythm, as an adequate number to capture the dynamics of every brain rhythm (fast and slows, Dimitriadis et al., 2013a). The center of the stepping window moved forwards every 20 ms and both intra and inter-frequency interactions within and between every possible pair of frequencies were reestimated leading to a quasi-stable in time static iPLV graph. In this manner, a series of 1,856 of iPLV graph estimates were computed per condition, frequencies (8 within frequency + 28 possible pairs of cross-frequency interactions = 36 combinations) and for each participant.

This procedure, the implementation details of which can be found elsewhere (Dimitriadis et al., 2010, 2012, 2013b, 2015a,c), resulted in eight time-varying iPLV graphs per participant (^TViPLV) for within frequency bands (Section 1 in Supplementary Material) and 28 time-varying iPLV graphs per participant (^TViPLV) for each possible cross-frequency pair (Section 2 in Supplementary Material), each serving as an instantaneous snapshot of the surface network. ^TViPLV tabulates iPLV estimates between every possible pair of sensors. For each subject, a 4D tensor [frequencies bands (36) × slides (1,856) × sensors (64) × sensors (64)] was created for each condition integrating subject-specific spatio-temporal phase interactions (Figure 1B).

Afterward, we applied surrogate analysis in order to reveal the dominant type of interaction for each pair of EEG sensors and at each snapshot of the ^TViPLV (Figure 1C—Section 3 in Supplementary Material). Finally, we tabulated both the strength and the type of dominant coupling mode in 2 3D tensor [slides (1,856) × sensors (64) × sensors (64)], one that keeps the strength of the coupling based on iPLV and a second one for keeping the dominant coupling of interaction using integers from 1 up to 36 {1 for δ, 2 for θ,…,8 for γ₂, 9 for δ−θ,…, 36 for β₂ − γ₂} (Figure 1D). The notion of phase-to-amplitude cross-frequency coupling (CFC) estimator has been used also in our previous studies (Dimitriadis et al., 2015b, 2016a,b,c, 2017).

On the previous section, we applied statistical filtering approach in order to explore: (a) if a given iPLV value differed from what would be expected by chance alone, (b) if a non-zero iPLV corresponded to non-spurious coupling, and (c) to integrate dominant intrinsic coupling modes into a single graph called dynamic integrated functional connectivity graph (DIFCG). Apart from 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.

The adopted topological filtering scheme based on orthogonal minimal spanning trees (OMST) (Figure 1D). Usage of the orthogonal MST graph leads to a better sampling of brain network preserving the advantage of MST that connects the whole network with minimum cost without introducing cycles and without differentiated strong from weak connections (Dimitriadis et al., 2017). For further details see also section 4 in Supplementay Material.

The outcome of both statistical and topological filtering scheme is demonstrated in Figure 1E for a static functional connectivity graph (FCG) from δ frequency of the first subject in the database. The resulting functional connectivity graph is sparser compared to the full-weighted (FCG) derived via the bivariate iPLV connectivity estimator. The FCG with dimension sensors x sensors illustrates the strength of the coupling in the surviving set of connections while the comodulogram demonstrates the probability distribution (PD) of the dominant coupling modes. The main diagonal refers to the PD of the within frequency interactions while the off-diagonal to the PD of the cross-frequency coupling namely the phase-to-amplitude couplings. The horizontal x-axis refer to how each the phase of the lower brain rhythms modulates the amplitude of the higher frequencies.

After applying the statistical and topological filtering approach, we estimated the global efficiency (GE) for each node and across time. This approach leads to a 2D matrix of dimensions 64 (sensors) × 1,856 (temporal segments) for each subject. This matrix tabulates the fluctuations of GE for each node and can be seen as network metric time series (NMTS) that keeps the importance of each node for integration of information (Dimitriadis et al., 2013a,b).

The estimation of GE based on the notion of shortest path length. However, we want to favor the selection of strongest functional couplings and for that reason the weights original functional connectivity graphs are inverted (1./FCG). A connection with strong coupling strength means that those two areas are functionally closed even in the case that are distant anatomically. By inverting the original weights, shortest path lengths between every pair of nodes will capture the strongest connections under the constraint of minimizing the overall distance cost. Distance cost (d_ij) of a shortest path length is defined the total sum of weights of the path that connects two nodes, here two EEG sensors.

Network GE reflects the overall efficiency of parallel information transfer within the entire set of 64 sensors and was estimated as the average sensor-specific GE value over all sensors using the following formula:

where d_ij denotes the distance of the shortest path length between target node i and the rest N − 1 j nodes. Here, we estimated the nodal GE which leads to 64 values per time-instant functional connectivity graph.

The derived 2D matrix based on nodal NMTS of GE will be modeled with the proposed method that is described in the following section. The first time instances of the topology of nodal GE is illustrated in Figure 1F.

This subsection serves as a brief introduction to our symbolization scheme, presented in greater details elsewhere (Dimitriadis et al., 2011, 2012, 2013a,b). 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 vector quantization approach (Dimitriadis et al., 2013a), if a preprocessing is added. The algorithm called non-negative matrix factorization (NNMF) (Marimpis et al., 2016). Here, we employed the NNMF as a first step for re-parameterization of the dynamic functional connectivity graph (DFCG) within a reduced space. This reduced space keeps the most informative ingredients for a parsimonious representation of the original connectivity patterns. Within this reduced NNMF-based low-dimensional space, a codebook of k prototypical functional connectivity states (i.e., connectivity microstates) is first designed by applying the neural-gas algorithm (Laskaris and Ioannides, 2001). 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).

In this way, the bulk of information contained in the time series of connectivity patterns is represented, in a parsimonious way, 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 code vectors, with a small reconstruction error E. The selection of parameter k reflects the trade-off between fidelity and compression level. In this work, it was set as the smaller value to achieve a reconstruction error lower than 4%. Finally, we need to mention that inside the codebook design there is an algorithmic seriation procedure that “ranks” the symbols. As a consequence, the symbolic time series closely follows the underlying functional connectivity dynamics. The neural-gas algorithm applied to the output of NNMF is called vector-quantization scheme (VQ) and the whole proposed scheme (NNMF-VQ). The derived symbolic times series that keep the information of network FCμstates (nFCμstates) are called hereafter as STS^NNMF-VQ. The combination of NNMF-VQ gives a lower distortion error compared to the VQ alone.

This approach leads to eight brain metastates for eyes-open and eight for eyes-closed based on the minimization of the average distortion error (Figures 1G, 2, 3). Both Figures 2, 3 illustrates the topology of the nodal GE across the eight nFCμstates in both conditions. The nFCμstates no: 4–8 for the eyes-open condition and the nFCμstates no: 6–7 for the eyes-closed demonstrates a topological similarity over parieto-occipital sites but with different activation profile expressed via the nodal GE.

In Figure 4A, we present four state transition vectors STS^NNMF-VQ from three subjects each one from three age groups. These STS^NNMF-VQ expressed the evolution of the eight nFCμstates across the experimental time. The derived transition matrix (TM) that tabulate the probability distributions of the favored transitions is demonstrated also in Figure 4A. Higher values in TM indicate preferred transitions between two states. Figure 4C demonstrates the evolution of dominant intrinsic coupling mode from pair of sensors F3-PO3 from the same subjects as in Figure 4A. The color refers to the strength of coupling namely the iPLV while the changes over the y-axis refer to the alterations of dominant intrinsic coupling modes (DICM). Across time (Figure 4C). On the right of these symbolic time series, we showed the comodulograms that tabulate the probability distribution of each of the 36 DICM for this particular pair of sensors (Figure 4D). The symbolic time series that encapsulate the DICM is a time series with values taken from 1 to 36 and it is called hereafter STS^DICM.

Once the DIFCG is formed and it is modeled via the NNMF-VQ scheme, relevant features can be extracted from the data based on the state-transition states and also on the DICM. Our approach adds in the literature also the notion of DICM which is an important characteristic of electroencephalography (Jirsa and Müller, 2013; Dimitriadis et al., 2015b, 2016a). DICM especially based on phase interactions within frequencies and phase-to-amplitude between frequencies are characteristic for EEG signal but not for fMRI. This section will describe the whole repertoire of features that were used here to build the developmental index.

The whole repertoire of features proposed here leads to a total number of 1 (TR) + 8 × 8 (TM) + 64 × 64 (FI) + 36 (comodulograms) + 1 (Flexibility Index) + 8 × 8 (entropy of Markovian trajectories) + 1 (Sample Entropy) = 4.262 features per subject. For that reason, we should apply an appropriate feature selection approach for both tasks: (a) prediction of age and (b) discrimination of young vs. middle age adults.

We estimated the distance correlation of 4262 features across subjects leading to a 4.262 × 4.262 matrix with positive weights. Distance correlation distcor ranges within [0, 1] in contrast to the correlation coefficient and can capture non-linear relationships (Szekely and Rizzo, 2014). Afterwards, we applied a graph clustering algorithm (Massimiliano and Pelillo, 2007) to extract clusters of features that are highly correlated within the cluster and low correlated between clusters. Our analysis revealed 25 clusters for the eyes-open condition and 43 for the eyes-closed. Finally, we extracted one feature per cluster with the criterion of high correlation with the age across the cohort. This procedure leads to 25 features for eyes-open and 43 features for eyes-closed.

As an appropriate method for build a model for predicted the age based on the selected features is the extension of SVM which is called SVM regression (SVR). Drucker et al. (1997) extended SVM to include also the regression approach for continuous real valued predictions. SVR retains the basic features of SVM classifier but the key difference is that in SVM classification a penalty exists for misclassified data points whereas in SVR a penalty is used for data points far from the regression line in a high dimensional space (Dosenbach et al., 2010).

For the binary classification of the two groups, we followed a different feature selection approach. By adopting the Laplacian Score (He et al., 2005), we estimated the weights of each feature that quantifies the discriminative power of each feature. On this procedure, we used the labels of 1 for young group and 2 for middle-age group. In order to apply a threshold, we re-estimated the LS of the features by shuffling the labels [1, 2] across the subject. We repeated the procedure 200 times and afterward we estimated the mean + 2.5 std (standard deviations) from the 1.000 × 4.263 LS weights. This thresholding criterion was applied to the original LS weights leading to the selection of 109 features.

As a proper classifier for the classification of age groups, we used the famous extreme learning machine (ELM). ELMs appeared as a proper choice due to its ability to manipulate difficult tasks without an extensive training set up (Huang et al., 2006). ELM are feedforward artificial neural networks with a single layer of hidden nodes where the weights that connect the input to the hidden nodes are randomly assigned and never updated (Huang, 2015).

We followed a leave-one out cross-validation scheme (LOOCV) in order to predict the age of a subject. This means that we trained the classifier with the selected features from N − 1 where N denotes the number of subjects and we predicted the age of the Nth. Repeated the whole approach N times, we predicted the individual age of every subject in the cohort. For the purpose of the prediction of age, we adopted the regression version of the Support Vector Machine, the Support Vector Regressor (SVR).

Our results revealed a very high prediction of age for eyes-open (R² = 0.60; y = 0.79x + 8.03) and lower for eyes-closed (R² = 0.48; y = 0.71x + 10.91) where y is a linear regression line applied to the (x, y) points with x being the true age of the subject and y the predicted age (Figures 5, 6).

Figure 7 illustrates the 25 selected features for the regression problem of prediction of actual age based on eyes-open. The 25 selected features referred to 24 FI per pair of sensors (Figure 7A) and 1 for the PD for α₂ (Figure 7B).

Using extreme learning machine (ELM classifier and a LOOCV), we succeeded to correctly classify young vs. middle-age group with 97% accuracy in eyes-open and 87.2% for eyes-closed (Table 1).

Figure 8 illustrates the distribution of selected features for the classification problem across the different sub-groups of features. Our feature selection analysis revealed 108 features referred to FI per pair of sensors (Figure 8A), 2 based on the PD for α₁ and α₂ (Figure 8B), 1 from the transition matrix between state 4 and 7 (Figure 8C), and 2 from the entropy reduction rate on between {4,7} and {8,5} pairs of states (Figure 8D). This makes a total of 103 features from a total of 4.262 initial pool of features.

Following an external validation of the whole approach in a cohort of 100 subjects, we performed reasonably well in continuous age prediction for eyes-open condition (R² = 069, y = 0.78x + 5.696) and eyes-closed condition (R² = 070, y = 0.81x + 3.217). The ELM classifier discriminated between young (18–37) and midde-aged subjects (40–60) with 98% accuracy (p < 10⁻⁷) for eyes-open and with 97% accuracy (p < 10⁻⁹) for eyes-closed (Table 2). Figure 9 illustrates the prediction for eyes-open and Figure 10 the prediction of eyes-closed.

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.