Dynamic Functional Connectivity of EEG: From Identifying Fingerprints to Gender Differences to a General Blueprint for the Brain's Functional Organization

An individual's brain functional organization is unique and can reliably be observed using modalities such as functional magnetic resonance imaging (fMRI). Here we demonstrate that a quantification of the dynamics of functional connectivity (FC) as measured using electroencephalography (EEG) offers an alternative means of observing an individual's brain functional organization. Using data from both healthy individuals as well as from patients with Parkinson's disease (PD) (n = 103 healthy individuals, n = 57 PD patients), we show that “dynamic FC” (DFC) profiles can be used to identify individuals in a large group. Furthermore, we show that DFC profiles predict gender and exhibit characteristics shared both among individuals as well as between both hemispheres. Furthermore, DFC profile characteristics are frequency band specific, indicating that they reflect distinct processes in the brain. Our empirically derived method of DFC demonstrates the potential of studying the dynamics of the functional organization of the brain using EEG.


INTRODUCTION
Due to advances in functional neuroimaging such as resting-state fMRI, there has been an increased focus on studying FC on an individual basis rather than at a group level (Dubois and Adolphs, 2016). It has recently been shown that FC profiles are highly individual, which emphasizes that it is indeed possible to draw inference on individual subjects (Finn et al., 2015). Similarly, the shape of an EEG spectrum qualifies as a statistical signature of a person (Näpflin et al., 2007). Another recent advance in the field of FC research is the recognition of information contained in the temporal dynamics of functional connectivity (Hutchison et al., 2013;Allen et al., 2014;Calhoun et al., 2014;Bassett and Sporns, 2017). This dynamic FC exhibits highly structured spatiotemporal states in which distinct patterns of FC recur across time and across subjects . These states are associated with cognition (Gonzalez-Castillo et al., 2015), consciousness (Barttfeld et al., 2015), neuropsychiatric disorders , personal traits (Shiino et al., 2017), and development (Hutchison and Morton, 2015). Furthermore, these FC states, when measured with fMRI, are associated with distinct EEG spectral signatures (Allen et al., 2018). As illustrated by other studies as well, EEG/MEG can be used to substantiate fMRI FC findings by providing an empirical link to a neural basis underlying BOLD signal changes (Tagliazucchi et al., 2012;Yu et al., 2016). Using a modeling approach, Honey and colleagues showed that the transient synchronization of neural populations is closely correlated with ultra slow fluctuations in the BOLD signal (Honey et al., 2007). Furthermore, they showed that when long time windows are used to calculate FC, the FC pattern is shaped by the underlying topology (Foster et al., 2016). Indeed, as pointed out by Allen and colleagues, the majority of EEG-fMRI integration studies have focused on EEG amplitude/FC modulations in relation to BOLD signal changes (de Pasquale et al., 2012;Chen et al., 2013;Liu et al., 2017;Allen et al., 2018). However, such a comparison is not straightforward because the relation between EEG power and actual spatial integration, as measured through spike-count correlations, is non-linear (Snyder et al., 2015). FC can also be studied using EEG/MEG exclusively and can be calculated over much shorter time-scales compared to fMRI (Van Diessen et al., 2015). This means that information about transient synchronization of neural populations, lasting between 100 and 300 ms (Varela et al., 2001), can be observed directly using EEG. Furthermore, using EEG, synchronization in different frequency bands can be analyzed separately, whereas in the BOLD signal these are lumped together. In this study we propose a novel EEG-based functional neuroimaging method to study the organization of transient functional connectivity. This method quantifies the correlation strength between time-varying FC signals, referred to as dynamic functional connectivity (DFC). Using EEG data from 105 healthy subjects scanned on two occasions up to 1 year apart, we show that, similar to an fMRI-derived FC profile (Finn et al., 2015), a DFC profile obtained from one session can be used to uniquely identify a given individual from a set of profiles obtained in a second session. Furthermore, we show that differences between male and female DFC profiles are substantial, which enables the prediction of a subject's gender based on his or her DFC profile. We also show that EEG-derived DFC profiles are organized according to a highly structured blueprint that is not only shared among people, but also across both hemispheres of the brain. Using additional data from a cohort of 81 Parkinson's disease (PD) patients, we show that our findings not only generalize to a state of health, but also of disease. Finally, our results also indicate that DFC could be used to distinguish PD patients from healthy controls.

Subject Information
The data set used in this study comprises three subsets from different studies conducted at the department of clinical neurophysiology, Basel University hospital and the Mathematics and Informatics department of Basel University. The first is from an ongoing SNFS study, investigating cognitive decline in Parkinson's disease ("Computer aided Methods for Diagnosis and Early Risk Assessment for Parkinson's Disease Dementia, " grant number 159682). Data from both healthy controls as well as PD patients was used. The other two subsets consist of data from the SNSF projects ("Improved prediction and monitoring of CNS disorders with advanced neurophysiological and genetic assessment, " grant numbers 124115 & 140338). All participants provided written informed consent in accordance with a protocol approved by the local ethics committee (Ethikkommission beider Basel). The 3 data sets used are described next. Data set 1. This data set consists of data from 33 healthy subjects (15 females, age 53-76) and 81 PD patients (28 females, age 45-84). EEG data were registered at baseline and again after 4 weeks and 6 months. At the time this study began, follow-up data from 25 healthy subjects and 57 PD patients was available. Baseline and 6 m data from the 33 healthy subjects and baseline and 4 w data from the 57 PD patients was used for the subject identification analysis. In case of the PD patients, we choses to use 4 w data instead of 6 m data because it is not yet known what influence a neurodegenerative disease such as PD has on a patient's DFC profile. Therefore, taking the shortest time-span between scan sessions minimizes the possibility of introducing bias.

Data Set 2
EEG data from 40 healthy subjects (30 females, age 20-49) was registered on 3 occasions over a period of 2 years (baseline, 1y, 2y). Baseline and 1y data from all 40 subjects was used for the subject identification analysis.

Data Set 3
EEG data from 41 healthy subjects (22 females, age 53-76) were registered at baseline and again after 3, 15, and 30 months. Baseline and 15m data were used for the subject identification analysis. For all 41 subjects a baseline scan was available. However, due to some dropouts, 15 m data were only available for 38 subjects. For the remaining 3 subjects 3 m data were used instead.
In all our analyses only data from healthy subjects were used. Data from the PD patients were used as an external validation set. For the analysis of differences between healthy subjects and PD patients, only data from sets 1 and 3 were used because the subjects from set 2 were significantly younger compared to the cohort of PD patients.

Data Acquisition and Preprocessing
EEG was recorded using a 256-electrode Sensor Net R . (Geodesics). First, the correct Sensor Net size was determined by measuring the subject's head circumference. Next, the net was placed over the subject's head such that the central electrode (Cz) was located at the crossing of the midline and lateral line. Subjects were scanned sitting upright and were instructed to sit still. EEG data from all subjects was acquired using the same protocol, which consisted of 4 consecutive sessions. In the first sessions the subject is instructed with an audio cue to alternatingly open and close his or her eyes. This sessions is followed by 15 min of eyes-closed resting state, and subsequently 5 min of eyes-open resting state sessions, after which the first sessions is repeated once more. Total registration duration was around 30 min including instructions. Only data from the 15 minutes eyesclosed were used in this study. Raw EEG signals were recorded with a sampling rate of 1,000 Hz and filtered with a high-order, linear-phase, finite-impulse response filter (MATLAB: Firls, 0.5-70 and 50 Hz notch, filter order: 4.8 × sampling rate). Only 170 out of 256 electrodes were used, midline electrodes and electrodes located in the face and around the neck were excluded (Figure 1).

Functional Connectivity Calculation
The EEG was first band-pass filtered after which phase estimation was performed using the Hilbert transform. The Hilbert transform was implemented using a 50% overlapping sliding window approach with 250 ms long windows. It should be noted that this short window length might approach a critical length for reliable phase estimation, especially for the δ frequency band. However, using longer windows up to 4 s had virtually no influence on the DFC coefficients (results not shown). Next, the continuous phase signals were divided into 0.25 s long nonoverlapping epochs which resulted in about N = 3,600 epochs per EEG recording. For each epoch, functional connectivity between pairs of EEG channels was calculated using the phase-lag-index (PLI, Stam et al., 2007). PLI is based on the asymmetry of the distribution of phase differences between two signals, and is calculated according to: T indexes the number of samples per epoch, and (t) represents the phase difference between two signals, for sample t in radians. This resulted in a symmetric connectivity matrix of size 170 by 170 per epoch per frequency band. A PLI value of 0 indicates no functional connectivity, whereas a value of 1 indicates maximal functional connectivity. Next, the betweenelectrode FC matrices are reduced to between-region FC matrices of size 22 by 22. Each element, or edge, represents the average FC strength between two scalp regions for a given time interval. Based on the lateralization of each scalp region, a distinction can be made between within-hemisphere edges (LL or RR) and between-hemisphere edges (LR or RL).

Dynamic Functional Connectivity Calculation
Because the FC matrices are symmetric, we first transformed each 22 by 22 matrix into a 231 long vector by extracting the upper triangle. Then, for each EEG recording and each frequency band, a 231 by N matrix of time-varying FC between all pairs of scalp regions is obtained, where N indicates the number of epochs. Next, DFC is calculated by simply calculating Pearson's correlation coefficients between all pairs of time-varying FC. In this way, a symmetric "dynamic FC (DFC)" matrix of size 231 by 231 is obtained for each EEG recording and each frequency band. Each element in the DFC matrix represents the correlation strength between two time-varying FC values, referred to as "DFC coefficient." Since the DFC matrices are symmetric, each matrix is converted to a 26565 long vector by extracting the upper triangle.

DFC Coefficients Types
To allow a more detailed analysis of identification results we used 2 ways of defining categories of DFC coefficients. One based on the number of scalp regions and another based on the lateralization of the functional connections. Since every DFC coefficient is defined as the Pearson's correlation coefficient between two functional connections, by definition each DFC coefficient involves either 3 or 4 scalp regions. In case of 22 scalp regions shown in Figure 2, this amounts to 4620 3-region and 21945 4-region DFC coefficients. Furthermore, we defined 4 categories of DFC coefficients based on which hemisphere each of the 3 or 4 scalp regions is located: 3025 unilateral within-hemispheric (LL-LL or RR-RR), 2970 unilateral between-hemispheric (LL-RR), 13310 bilateral-unilateral (LR-RR or LR-LL), and 7260 bilateral-bilateral (LR-LR).

Identification Analysis
In our first identification analysis we used the same approach as Finn et al. (Finn et al., 2015). A database was created that consisted of all the subjects' baseline DFC vectors, , consisting of all the subjects' follow-up DFC vectors was created. In case of subject identification based on all five frequency bands, X i and Y i are constructed by concatenating the five 26565 long DFC vectors into a single 132825 long vector. For single-frequency band identification X i and Y i can either be based on the same frequency band (within-frequency band) or on different frequency bands (between-frequency band). To predict a subjects' identity in one database based on the other database, Pearson correlation between pairs of DFC vectors was defined as similarity measures and was calculated between all pairs of DFC vectors from the two databases. This resulted in a 14 by 144 similarity matrix where each element represents the similarity between a DFC vector from the baseline database and a DFC vector from the followup database. Subsequently, each DFC vector from one database was assigned the identity of the DFC vector from the other database that was maximally similar. Identification accuracy was defined as the percentage of subjects for which his or her identity was correctly identified. Following the approach of Finn et al. (2015), we also performed non-parametric permutation testing to assess the statistical significance of identification performance. In each iteration, identification as described above was performed whereby the identity of subjects in one database was permuted. This permutation procedure was repeated 10,000 times.

Gender Classification Analysis
In our gender classification analysis we used baseline data from all 144 healthy subjects. We first investigated the magnitude and extend of DFC profile differences between males and females by assessing how many DFC coefficients were statistically different (t-test) using a significance level of P < 3.76 × 10 −7 (P < 0.05 Bonferroni corrected for 132,825 comparisons). Considering the heterogeneity of our data set with respect to age, male and female subsets were matched for age by randomly excluding 20 FIGURE 1 | Our electrodes were grouped according to the topographical localization of the channels and then divided into 22 non-overlapping regions of interest corresponding to the frontal, central, temporal, parietal and occipital areas (bilaterally) (Hatz et al., 2017). female. Presented results are obtained by averaging 500 Monte Carlo repetitions. To determine with which accuracy a subjects' gender can be predicted based on his or her DFC profile, we applied feature selection and classification using LASSO (least absolute shrinkage and selection operator) (Tibshirani, 1996). Again, presented results are obtained by averaging 500 Monte Carlo cross validation repetitions. In each repetition, 80% of the subjects are assigned to the training data set and the remaining 20% are assigned to the test data set. Furthermore, male and female subsets of the training data sets where matched for age.

Common DFC Blueprint Analysis
To analyze DFC profile characteristics that are shared among individuals we constructed an average gender-matched baseline DFC profile from all 144 healthy subjects. Considering the heterogeneity of our data set with respect to gender and age, an average DFC profile was created 500 times by randomly selected subset of subjects, matched for age and gender. Subsequently the 500 averaged DFC profiles where in turn averaged to obtain a population averaged DFC profile.

Influence of Epoch Length on DFC Information Content
To analyze the influence of epoch length used for PLI calculations on the reliability of DFC coefficients, we repeated the identification and gender classification analysis using epoch lengths ranging from 250 ms up to 32 s. Additionally, we used the EEG recordings from those healthy subjects included in the identification analysis to show that with increasing epoch lengths the test-retest error also increases. The test-retest error for a given DFC coefficient was defined as the absolute difference between its values from both recordings. Subsequently, for each DFC coefficient the errors for all subjects were averaged; the distributions of all 26,565 DFC coefficients are shown in Figure 5.
To illustrate the effect of EEG recording duration, we simulated shorter durations by simply using only the first N minutes of each recording for DFC calculation. N was varied between 1 and 10 min.

RESULTS
In the first part of this study we used data of 105 healthy subjects to show that intra-individual differences in DFC profiles are  Within-subject correlation in α 1 is only significantly higher compared to β. Between-subject correlation in alpha 1 is significantly lower compared to the other 4 frequency bands (one-tailed paired t-test, *p < 0.01, **p ≪ 0.0001). Boxplots represent median (stripe), 25th and 75th percentiles (box), 1.5 times the interquartile range (whiskers), and outliers (crosses). Box colors correspond to the ID rate scale of (A). (C) Number of DFC coefficients per frequency band that are statistically significant different between males and females based on a significance level of α < 3.76 × 10 −7 . substantial and can reliably be observed, for example, they can act as an identifiable "fingerprint." For each subject, a high-density resting-state EEG was obtained twice over a period ranging from 6 to 15 months. Detailed subject information can be found in the online material. For each subject, 12 or 15 min of eyes-closed resting state EEG data were recorded using a 256-electrode EEG recording device (Electrical Geodesics). Signals from a subset of 170 electrodes were used; midline electrodes and those located in the face and around the neck were excluded. By bandpass filtering the EEG recordings, a DFC profile was calculated for each of 5 commonly used frequency bands (δ: 1-4 Hz, θ : 4-8 Hz, α 1 : 8-10 Hz, α 2 : 10-10 Hz, and β: 13-30 Hz). This was done because distinct frequency bands exhibit characteristic changes in response to sensory, motor, and cognitive events (Engel et al., 2001;Varela et al., 2001;Buzsáki and Draguhn, 2004). It is therefore expected that DFC profiles calculated in different frequency bands will differ. Next, from the bandpassfiltered EEG, the instantaneous phases were estimated using the Hilbert transform with a 50% overlapping sliding window. Then, for 250 ms long, non-overlapping epochs, Functional Connectivity (FC) between all possible pairs of the 170 EEG channels was assessed using the phase-lag-index (PLI, Stam et al., 2007). For an EEG recording consisting of N epochs, this resulted in a series of N symmetric FC matrices of size 170 × 170. Because the localization of EEG sources is in the order of centimeters (Cohen, 2017), we reduced the 170 × 170 inter-electrode FC matrices to 22 × 22 inter-region FC matrices. Each element, or edge, represents the time-varying FC strength between two scalp regions. Thus, given 22 scalp regions, there are 231 unique edges. Finally, Pearson's correlation coefficients were calculated for all edge-pairs in order to obtain a 231 × 231 symmetric dynamic FC matrix where each element represents the correlation strength between two edge-pairs, referred to as a "DFC coefficient." A more detailed description of the DFC calculation is provided in the section 2. DFC matrices were calculated for each EEG and for each frequency band such that each subject has 5 matrices reflecting his or her DFC profile during a session. Identification was performed using the same approach as Finn and colleagues (Finn et al., 2015). In brief, 2 datasets were constructed such that each subject had a single DFC profile in each set. Identification implied that each DFC profile in one set ("set1") was assigned the identity of the DFC profile in the second set ("set2") that was maximally similar. Similarity was defined as the Pearson's correlation coefficient between the two vectorized DFC profiles. Once all identities had been predicted, the overall identification rate was calculated as the percentage of scans for which the predicted identity matched the true identity. This process was then repeated with the roles of set1 and set2 reversed. Subsequently, identification performance was averaged over these two attempts.

All Frequency Bands
First, identification was performed using all DFC profiles obtained from the 5 frequency bands. The Identification rate was 89/105 (85%) and 85/105 (81%) based on attempt 1 and 2, respectively. The same permutation test as used by Finn and colleagues was performed which, across 10,000 iterations, yielded a maximal identification rate of 6/105 or less than 5%. Therefore, obtaining 84 correct identifications was associated with p << 0.0001.

Single Frequency Bands
To test the hypothesis that different information, with respect to individual subject discriminability, is contained in different frequency bands, we subsequently tested identification accuracy based on DFC profiles from single frequency bands ( Figure 3A). In case different information is contained in each frequency band, identification performance is expected to be lower when using single frequency bands instead of when using all 5. Furthermore, performance is expected to be reduced even more when different frequency bands are used for each set (see section 2). For within-frequency-band identification, performance indeed dropped significantly (74%, paired t-test α 1 vs. all frequency bands, t 209 = −3.37, P < 0.001). Furthermore, α 1 identification was significantly better compared to δ, θ , and β (paired t-test α 1δ p = 2.5×10 −9 , α 1 -θ p = 0.028, α 1 -α 2 p = 0.078, α 1β p = 0.0017). Furthermore, crossfrequency-band identification performance was significantly lower compared to within-frequency-band identification (ttest within-vs. cross-frequency-band identification, t 13 = 8.5, P < 10 −5 ). To asses whether better α 1 performance can be attributed to higher within-subject correlation and/or lower between-subject correlation, we compared within-and betweensubject correlations between frequency bands (Figure 3B). Compared to α 1 , within-subject correlation was only significantly lower for β (paired t-test, p < 0.01) whereas betweensubject correlation was significantly higher for all 4 frequency bands (paired t-test, p ≪ 0.0001). Superior identification performance for α 1 can therefore be attributed primarily to a relatively low between-subject similarity compared to other frequency bands.

Gender Differences
In the previous analysis we found that an individual's DFC profile is unique, and can as such act as an identifying fingerprint. In the following analysis we investigated whether the DFC profile of a patient contains other relevant information. An important but often neglected topic in neuroscience is sexrelated influences on brain function (Cahill, 2006). Recent functional and anatomical connectivity studies reported both morphometric (Gong et al., 2009;Chekroud et al., 2016), as well as functional connectivity (Tomasi and Volkow, 2012), differences between the male and the female brain. Our analysis evaluated sex-related differences, and whether DFC profiles can be used to predict an individual's gender. We used the baseline EEGs from the same 105 subjects as in the previous analysis and an additional 8 subjects for whom only a single EEG was available.
As a first step we used a t-test to assess how profound the sex differences of DFC profiles are. Considering the heterogeneity of our data set with respect to age, we performed 500 Monte-Carlo simulations such that in each test the male and female subset was matched according to age. On the basis of a significance level of α < 3.76 × 10 −7 (0.05 Bonferroni-corrected for 132,825 comparisons), we found that the number of DFC coefficients that were significantly different between male and female ranged between 3,531 (13%) and 433 (1.6%) for the β and α 1 frequency band respectively ( Figure 3C). Interestingly, for each frequency band, virtually all of these DFC coefficients were higher for females. The magnitude of gender differences (male -female) ranged from −0.21 to +0.12.
The differences in the number of statistically significantly different DFC coefficients per frequency band indicate that gender differences are not the same for each frequency band. However, it is possible that the frequency bands only differ from each other with respect to gender differences in the magnitude of these differences. Therefore, we tested the hypotheses that gender differences are frequency band specific. To this end, for each frequency band f, we calculated an average male and average female DFC profile and subsequently subtracted one from the other in order to obtain a "gender difference DFC profile" ( DFC f ). To test the hypothesis that gender differences are symmetric with respect to both hemispheres, we correlated FIGURE 4 | Gender differences in the α frequency band. (A) DFC coefficients typically stronger in female brains (left) or male brains (right). For visualization, both sets represent the top 99 percentile of DFC coefficients for which the gender difference is largest. Text on the circle represent all 231 functional connections and we made a distinction between unilateral (LL and RR) and bilateral (LR) connections. Except for connections between cross-hemisphere homologs (LR), connections are ordered such that the left and right half are mirror images. Each line represents a DFC coefficient and its color distinguishes between the 4 lateralization based categories. Additionally, on the outer circle we used a bar for each connection to summarize in how many DFC coefficients it was involved. Furthermore, each bar is divided in a blue and a red part indicating the number of 3-region and 4-region DFC coefficients, respectively. (B) In the barplots the same data is plotted as the incidence of scalp regions in the set of DFC coefficients, averaged over both hemispheres. L, left hemisphere; R, right hemisphere; FP, fronto-polar; FM, frontal-mid; FL, frontal lateral; CM, central mid; CL, central lateral; TA, temporal anterior; TP, temporal posterior; PL, parietal lateral; PM, parietal mid; PO, parietal occipital; OC, occipital.
DFC with its mirrored version (see section 2). Correlation coefficients ranged from 0.9 (α 1 ) to 0.94 (δ) which indicates that gender differences are indeed symmetric. We next determined for which DFC coefficient categories, and for which involved brain regions, the gender differences were most prominent. We restricted our analysis to those DFC δ coefficients for which the magnitude of gender difference were in the bottom 1 and top 99 percentile (Figure 4).
To allow a more detailed analysis of the identification results, we used 2 ways of defining categories of DFC coefficients: one based on the number of regions involved in a DFC coefficient, and the other based on the lateralization of the regions involved. Another distinction can be made based on the number of regions a DFC coefficient is based on, which is either 3 or 4 regions. Lateralization is categorized according to which hemisphere each of the 3 or 4 brain regions is located: unilateral withinhemispheric (LL-LL or RR-RR), unilateral between-hemispheric (LL-RR), bilateral-unilateral (LR-RR or LR-LL), and bilateralbilateral (LR-LR).
In case of the DFC coefficients, which were typically much stronger for females, we found that the majority were either unilateral within-hemispheric (LL-LL or RR-RR) or bilateral-bilateral (LR-LR) 4-region DFC coefficients. These DFC coefficients mainly involved parietal and occipital regions. On the other hand, DFC coefficients that were typically stronger for males were primarily bilateral-unilateral (LR-RR or LR-LL) 3-region DFC coefficients and involved primarily frontal as well as central mid, and occipital regions.

Gender Classification
Based on the sheer number and magnitude of gender differences, it is reasonable to expect that an individual's gender can be determined based on his or her DFC profile. In the next section, we evaluated if this is indeed possible. Since not all DFC coefficients are gender-dependent, we use LASSO (least absolute shrinkage and selection operator) to evaluate if a small subset of DFC coefficients is also sufficient for gender classification. We used Monte-Carlo cross-validation where in each iteration LASSO was applied to 80% of the data in order to fit a regression model between DFC coefficients and gender indicated by 0 (male) or 1 (female). This regression model was then used to predict the gender of the left-out 20%. Prediction performance was assessed as the area under the ROC curve (AUC) of the 500 Monte-Carlo repetitions. On average, 30 (range: 12-103) DFC coefficients were selected to fit the regression model. Median classification performance was AUC = 0.93 (range: 0.6-1). In total, 3152 DFC coefficient were selected by Lasso in at least 1 of 500 repetitions. As could be expected based on the results described in the previous section, DFC coefficients from the δ (27%) and β (23%) frequency bands were chosen more frequently than the other 3 frequency bands (16% each). Shorter epochs improve gender classification performance. Results from 500 Monte Carlo cross-validation repetitions are shown. Boxplots represent median (stripe), 25 th and 75th percentiles (box), 1.5 times the interquartile range (whiskers), and outliers (crosses). (C) Epoch length is inversely related to DFC coefficient estimation accuracy. (D) EEG duration is also inversely related to the duration of the EEG. Boxplots represent the distributions of differences in DFC coefficients between Rest1 and Rest2, averaged over 105 subjects for the α 1 frequency band (26565 DFC coefficients).

A Shared DFC Blueprint Among Individuals and Hemispheres
We further investigated the common properties of DFC organization. To this end we collapsed healthy subjects' DFC matrices into a single average DFC matrix (DFC). Naturally, a DFC was created for each frequency band separately. Similarly to our gender differences analysis, we here also calculated Pearson's correlation coefficients between DFC profiles of the different frequency bands, as well as between DFC and their mirrored versions. Correlation coefficients between frequency bands ranged between 0.90 (δ vs. β) and 0.99 (α 1 vs. α 2 ). Furthermore, the further apart two frequency bands were, the lower the correlation between DFC profiles. This indicates that there is a common blueprint shared among frequency bands on top of which frequency-band-specific characteristics are superposed. Conversely, correlation coefficients between DFC profiles and their mirrored version are above 0.995 for each frequency band, which indicates that both hemispheres share the same blueprint. To study this common blueprint in more detail we first evaluated DFC in terms of the different DFC coefficient categories. Considering the strong similarities between DFC from the different frequency bands, we restrict the following analysis to the α 1 frequency band. We found that 3-region DFC coefficients are generally stronger compared to the 4-region DFC coefficients (Figure 6, t 26563 = 87.8, P ≪ 0.001). Similarly, unilateral withinhemispheric (LL-LL and RR-RR) DFC are stronger compared to the other 3 categories (Figure 6, t 26563 = 46, P ≪ 0.001). Next, we determined which DFC coefficients were in the top we made a distinction between unilateral (LL and RR) and bilateral (LR) connections. Except for connections between cross-hemisphere homologs (LR), connections are ordered such that the left and right half are mirror images. Each line represents a DFC coefficient and its color distinguishes between the 4 lateralization based categories. Additionally, on the outer circle we used a bar for each connection to summarize in how many DFC coefficients it was involved. Furthermore, each bar is divided in a blue and a red part indicating the number of 3-region and 4-region DFC coefficients, respectively. (B) In the barplots the same data is plotted as the incidence of scalp regions in the set of DFC coefficients, averaged over both hemispheres. L, left hemisphere; R, right hemisphere; FP, fronto-polar; FM, frontal-mid; FL, frontal lateral; CM, central mid; CL, central lateral; TA, temporal anterior; TP, temporal posterior; PL, parietal lateral; PM, parietal mid; PO, parietal occipital; OC, occipital. (c) Differences between the several DFC coefficient categories based on the average DFC profile of 7 all healthy individuals. Boxplots represent median (stripe), 25th and 75th percentiles (box), 1.5 times the interquartile range (whiskers), and outliers (crosses). Box colors and legend indicate DFC coefficient category. 99 and bottom percentile. As expected we found that unilateral within-hemispheric (LL-LL and RR-RR) as well as 3-region DFC coefficients were overly represented in the 99 percentile, and virtually absent in the 1 percentile (Figure 6). Conversely, the opposite was found for unilateral between-hemispheric and 4region DFC coefficients (Figure 6). When analyzing the 1 and 99 percentile subsets in more detail, we found that functional connections involving frontal and parietal regions appeared most frequently in the subset of strongest DFC coefficients (Figure 6). On the other hand, in the subset of weakest DFC coefficients, functional connections involving, and especially between, temporal-anterior and parietal-lateral regions appeared most frequently.

Influence of Epoch Length on DFC Information Content
For very long time-windows, functional connectivity patterns closely resemble structural connectivity patterns (Foster et al., 2016). However, when functional connectivity is calculated for shorter time windows the dynamics of functional connectivity is revealed which enables analyses such as the one presented in this study. A question arises about what the effect of increasing epoch length is on the amount of information contained in the DFC coefficients. To this end, we repeated the subject identification and gender classification analyses using 7 additional epoch lengths to calculate functional connectivity: 500 ms, 1, 2, 4, 8, 16, and 32 s. Increasing the epoch length resulted in lower identification rates (Figure 7) as well as lower gender classification performance (Figure 7). This means that with decreasing epoch lengths the DFC coefficients can be estimated with higher accuracy, as expressed by the smaller differences (RMSE) in DFC coefficient values between Rest1 and Rest2 (Figure 7). Another factor that is expected to play a role in accuracy of which the DFC coefficients can be identified is the duration of the EEG recording. To test the hypotheses that shorter EEGs also result in reduced accuracy, we calculated DFC profiles with 250 ms long epochs, using truncated EEGs ranging from 1 to 10 min. We observed that shorter EEGs, and hence less data, indeed result in decreased accuracy (Figure 7).

External Validation: Parkinson's Disease Patients
As a stronger proof of generalizability we applied our subject identification and gender classification models on a completely independent validation set consisting of resting state EEGs from 81 patients diagnosed with Parkinson's disease (PD). Of these 81 patients, 58 had 2 EEGs recorded over a period of 4 weeks. Identification based on all 5 frequency bands resulted in an identification rate of 79%, which is in accordance with the 82% that was achieved using the 105 HC subjects. Additionally, an identification rate of 80% was achieved when identification was based on the joint set of 163 healthy subjects and PD patients.
We next applied the gender classification models from the analysis described in the "Identification and gender classification using all frequency bands" section to the 81 PD patients as well. The models trained on the data from healthy subjects could also successfully discriminate between males and females: AUC = 0.89 ± 0.02.
Finally we investigated if DFC could also distinguish between PD patients and healthy subjects. Considering the sex-related differences in DFC coefficients, and the fact that we have almost twice as many male as female PD patients, we restricted our analysis to those DFC coefficients that were not significantly different in our gender analysis (p > 0.05). This yielded a total of 41616 DFC coefficients for which we assessed differences between PD patients and healthy subjects using a t-test. P-values below 1.2 × 10 −6 (0.05 Bonferroni-corrected for 41,616 comparisons) were considered statistically significant. We also excluded a subset of healthy subjects from this analysis because they were much younger than the PD patients. Only 13 DFC coefficients survived Bonferroni correction (5, 4, 3, and 1 from the δ, θ , 1, and β frequency band respectively). All of these DFC coefficients were higher for healthy subjects, and all except 1 involved either a left or right occipital region (Figure 7).

DISCUSSION
In practice PLI and other functional connectivity measures are typically estimated by averaging over trials and/or time (Aydore et al., 2013). Here we propose an alternative method of functional connectivity analysis. We demonstrate that correlation coefficients between pairs of functional connectivity timeseries contain meaningful information about the functional organization of the brain. First, we demonstrated that an individual's DFC profile is unique and reliably observable which makes it possible to use it to identify an individual from a large group of subjects. Second, we also demonstrate that an individual's DFC profile is highly gender-specific, which makes it possible to determine a subject's gender solely based on his or her DFC profile. Third, we show that there is a general blueprint for DFC organization that is not only shared among individuals, but also between left and right hemispheres. These observations are in accordance with our current understanding of functional brain organization, namely, that there is an intrinsic standard architecture for functional brain organization on top of which task- (Cole et al., 2014) and individual-specific characteristics (Finn et al., 2015) are superposed. Furthermore, this architecture is also shared between left and right hemispheres which agrees with observations that similar bilateral functional networks develop independently (Tyszka et al., 2011). Fourth, our results also generalize to a population of PD patients, and furthermore indicate that there are specific changes in DFC profiles from PD patients. This observation in particular is promising because it suggests that DFC can be used to study neurodegenerative diseases. Last, we showed that DFC coefficients dependent on the frequency band used to calculate PLI values. A large body of research links specific patterns of oscillations to perceptual, cognitive, motor, and emotional processes (Cohen, 2017). It is therefore reasonable to expect that DFC profiles calculated based on different frequency bands are associated with different neuronal processes.

DFC Coefficient Interpretation
A DFC coefficient is defined as the correlation strength between two PLI time-series. First, it is important to recognize that PLI calculated over a period of several hundreds of milliseconds cannot reflect an intrinsic property of the brain's functional organization. This is because the magnitude of epoch length and the typical duration of phase-locking are the same (Varela et al., 2001). Consequently, during most epochs the instantaneous phases of two signals only cross once, twice, or not at all. Assuming that only 1 phase-crossing occurs, the resulting PLI value only depends on when during this epoch the phase-crossing took place. Naturally, in such a scenario a single PLI value is determined by chance. However, if two PLI time-series are positively correlated, it means that the phase crossing in both sets of signals repeatedly occur around the same time. Hence, a DFC coefficient can best be understood as a measure of the timing relationship of the phase crossings between sets of EEG signals. Contrary to single PLI values calculated over short epochs, it is reasonable to assume the timing relationship between phasecrossing does reflect an intrinsic property of the brain.
Our observation that the majority of DFC coefficients are stronger in females fits to other studies demonstrating stronger anatomical (Gong et al., 2009) as well as functional connectivity (Tomasi and Volkow, 2012). However, we also found a subset of DFC coefficients which were stronger in males. This is in line with observations that some functional connections are stronger in females while others are stronger in males (Biswal et al., 2010). If our hypothesis of higher DFC coefficients reflecting shorter time-intervals between phase-crossings is correct, differences in brain size might also partly explain differences in DFC coefficients between the sexes (Hänggi et al., 2014). Why gender differences are most profound in the δ (1-4 Hz) and β (13-30 Hz) frequency bands, and least profound in the α 1 (8-10 Hz) frequency band, might be explained by the opposite pattern observed in our subject identification analysis, where between-subject DFC correlations were lowest in the α 1 (8-10 Hz) frequency band. However, caution is warranted when trying to find a physiological interpretation for a measure of statistical dependence between EEG signals. A method such as DFC might indicate which EEG features are relevant for deeper investigation. For example, it would be interesting to use an effective connectivity approach, such as direct causal modeling (Friston, 2011), to better understand the mechanism of brain connectivity.

Additional Considerations
The development of our DFC method is an empirical extension of existing FC methods, and therefore based on several assumptions such as the use of epochs, averaging over 22 scalp regions, PLI as functional connectivity measure, and choice of frequency bands. However, it is not straightforward that these assumptions are valid, let alone optimal, when studying the dynamics of FC. For example, one could argue that for 250 ms-long epochs, it does not make sense to use PLI as a FC measure. Nevertheless, we showed that our DFC method provide meaningful information about an individual. Therefore, our DFC method can best be seen as the starting point of developing theoretically sounder methods to study the dynamics of phase interactions.
We deliberately used minimal pre-processing, in the form of band-pass filtering, and also used the complete resting state EEG recording. A common approach in quantitative EEG analysis is to let one or more EEG experts select parts of an EEG recording which are then used for subsequent analysis, instead of using the complete recording. Although this might prevent the results being negatively influenced by various artifacts, it also obscures the process of how the results are obtained. Furthermore, it enables cherrypicking of EEG data such that better results are obtained. However, it goes without saying that careful selection of highquality EEG data should remain part of any quantitative EEG analysis, especially when used to draw inference on individual subjects.

CONCLUSION
These results underline the potential of EEG-based DFC as a tool to investigate brain functioning on both a population and individual level. This could in turn lead to the discovery of new "neuromarkers" for cognitive behavior, psychiatric conditions, and cognitive decline in neurodegenerative diseases. Furthermore these results underline the validity of phasebased functional connectivity methods such as PLI, be it as an intermediate step for a more advanced quantification of functional brain organization.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by Ethikkommission Beider Basel, Nordwestschweiz. The patients/participants provided their written informed consent to participate in this study.

AUTHOR CONTRIBUTIONS
GB, UG, PF, and VR conceptualized the study. GB, MC, FH, and AM collected and processed the data. GB designed and performed the analyses, wrote the data visualization tools, and wrote the manuscript, with comments from all authors. GB and FH wrote the analysis tools. PF and VR provided guidance with data interpretation. The study was performed under the supervision of PF and VR.

FUNDING
This work was supported by the Schweizerischer Nationalfonds (SNF grant 159682).