Prediction and classification of sleep quality based on phase synchronization related whole-brain dynamic connectivity using resting state fMRI

Recently, functional network connectivity (FNC) has been extended from static to dynamic analysis to explore the time-varying functional organization of brain networks. Nowadays, a majority of dynamic FNC (dFNC) analysis frameworks identified recurring FNC patterns with linear correlations based on the amplitude of fMRI time series. However, the brain is a complex dynamical system and phase synchronization provides more informative measures. This paper proposes a novel framework for the prediction/classification of behaviors and cognitions based on the dFNCs derived from phase locking value. When applying to the analysis of fMRI data from Human Connectome Project (HCP), four dFNC states are identified for the study of sleep quality. State 1 exhibits most intense phase synchronization across the whole brain. States 2 and 3 have low and weak connections, respectively. State 4 exhibits strong phase synchronization in intra and inter-connections of somatomotor, visual and cognitive control networks. Through the two-sample t-test, we reveal that for the group with bad sleep quality, state 4 shows decreased phase synchronization within and between networks such as subcortical, auditory, somatomotor and visual, but increased phase synchronization within cognitive control network, and between this network and somatomotor/visual/default-mode/cerebellar networks. The networks with increased phase synchronization in state 4 behave oppositely in state 2. Group differences are absent in state 3, and weak in state 1. We establish a prediction model by linear regression of FNC against sleep quality, and adopt a support vector machine approach for the classification. We compare the performance between conventional FNC and PLV-based dFNC with cross-validation. Results show that the PLV-based dFNC significantly outperforms the conventional FNC in terms of both predictive power and classification accuracy. We also observe that combining static and dynamic features does not significantly improve the classification over using dFNC features alone. Overall, the proposed approach provides a novel means to assess dFNC, which can be used as brain fingerprints to facilitate prediction and classification.


Introduction
Functional connectivity (FC) profiles, which can be derived from the functional magnetic resonance imaging (fMRI) time series, have provided a promising measure to characterize the connectivity patterns within the human brain Smith et al., 2009 ;van den Heuvel and Hulshoff Pol 2010 ;Hutchison and Morton 2015 ;Smitha et al., 2017 ;Cohen 2018 ). Due to the unique feature of FC pat-this type of FC is referred to as functional network connectivity (FNC) ( Jafri et al., 2008 ).
A majority of FC/FNC studies assume that the FC/FNC is static throughout the entire fMRI session, resulting in a single measure averaging all time points. Although static FC/FNC can characterize behaviorally relevant individual variations, it loses information at a temporal scale. A growing consensus in the neuroscience field is that individuals engage in different mental activities at different time points ( Cohen 2018 ;Lurie et al., 2020 ). Accumulated evidence has shown that the spontaneous dynamic fluctuations in brain activity are consistent with dynamic changes in FC/FNC during the course of the fMRI session ( Chang and Glover 2010 ;Hutchison et al., 2013 ;Rashid et al., 2014 ). The dynamic patterns in FC/FNC provide additional valuable information about individual differences in human behavior and cognition. Due to the unconstrained mental activity, dynamic patterns of functional brain connectivity are potentially more prominent in the resting state ( Shehzad et al., 2009 ;Allen et al., 2014 ;Betzel et al., 2016 ). DFNC methods have been used in the analysis of electroencephalography (EEG)/magnetoencephalography (MEG) ( O'Neill et al., 2018 ;Coquelet et al., 2020 ), which provide better temporal information on brain activity in comparison with fMRI.
Considering the fluctuating information within the brain, dynamic functional network connectivity (dFNC) analysis has been proposed ( Sakoglu et al., 2010 ;Allen et al., 2014 ). In dFNC analysis, a data-driven approach, namely group independent component analysis (GICA), is used to decompose the brain into independent component networks (ICNs) ( Calhoun et al., 2001 ). Next, the dFNC among brain networks can be captured by a sliding window technique, whereby the time series of ICNs are segmented into windows and Pearson correlation coefficients (PCCs) between pairs of time series are computed within each window ( Sakoglu et al., 2010 ;Handwerker et al., 2012 ;Allen et al., 2014 ;Barttfeld et al., 2015 ;Choe et al., 2017 ). Following this procedure, k-means clustering is employed on all the dynamic correlation matrices to identify a finite number of discrete dFNC patterns or states within subjects. This dFNC analysis framework is often used as chronnectome, focusing on identifying time-varying but recurring patterns of coupling among brain regions .
From the dFNC analysis, both direct and indirect measures can be extracted. Due to the simplicity, connectivity strength has been widely used as a direct measurement. Researchers have successfully identified between-group differences of connectivity strengths for some dFNC states in terms of age, sex, intelligence, diseases and etc., which cannot be reflected by the static FNC de Lacy et al. 2019 ;Espinoza et al., 2019a ;Xia et al., 2019 ). However, there is a challenge in identifying individual difference and predicting behaviors and cognitions based on the connectivity strengths of each dFNC state. More specifically, a considerable number of research has shown that one or more dFNC states are absent for many individuals after the k-means clustering of estimated dFNCs Liu et al., 2017 ;Cai et al., 2018 ;Vergara et al., 2018 ;Zhi et al., 2018 ;Espinoza et al., 2019b ). Although a fuzzy meta-state approach can alleviate the stateabsent problem by allowing a subject's state to be represented by varying degrees of multiple states , the challenge is still not resolved. Summary measures such as the variance of dFNC matrices ( Liu et al., 2017 ) and amplitude of low-frequency fluctuation of FNC ( Hare et al., 2018 ) are indirect measures, which have been used to classify patients with schizophrenia or idiopathic generalized epilepsy from health controls at individual level. But these summary measures are less interpretable for brain functional organization. Another indirect measure namely regression coefficients against dFNC states is proposed recently by Rashid et al. ( Rashid et al., 2016 ), where they extracted 5 dFNC states for schizophrenia, bipolar and healthy subjects respectively. However, the regression coefficients cannot provide direct measure of functional connectivity.
The above described dFNC analysis approaches identify recurring functional connectivity patterns, which are mainly based on the linear correlation. However, the brain is a complex nonlinear dynamical system, which exhibits intriguing features such as oscillations and synchronization of neural activity ( Aydore et al., 2013 ). Hence, as a nonlinear measure, phase synchronization can be incorporated into the FNC analysis ( Pereda et al., 2005 ). Currently, a few studies have examined the static ( Cordova-Palomera et al., 2016 ;Gravel et al., 2018 ;Wang et al., 2019 ) or dynamic phase synchronization ( Glerean et al., 2012 ;Ponce-Alvarez et al., 2015 ;Demirtas et al., 2016 ;Omidvarnia et al., 2016 ;Pedersen et al., 2017 ) for fMRI. Among these studies, only one study belongs to chronnectome research ( Ponce-Alvarez et al., 2015 ), but the identified recurring patterns are binarized ones from phase difference, which will result in a loss of information. Moreover, Ponce-Alvarez et al. examined the recurring patterns for all scanning fMRI sessions, but not for each subject. It should be noted that all these phase synchronization studies for fMRI involve no behavior and cognitive performance prediction/classification. As stated in Zhang et al's work , amplitude-based static and dynamic FNC/FC have been used for predicting behavior and cognitive performance, but it is still an open problem on whether neural synchronization characterized by analytical signal (e.g., instantaneous phase) can serve as an alternative.
In order to extract phase synchronization based dFNC patterns, the fMRI time series should be transformed into narrowband signals. The empirical mode decomposition (EMD) has been identified as a selfadaptive filtering method to decompose the acquired signal into a collection of intrinsic mode functions (IMFs) with narrowbands ( Huang et al., 1998 ). The EMD has been used to analyze static functional connectivity in resting-state fMRI ( Moradi et al., 2019 ;Qian et al., 2018 ;Zhang et al., 2015 ). In a recent study, Goldhacker et al. proposed an approach for identifying dynamic connectivity states in rs-fMRI based on IMFs of ICA time courses and K-means clustering ( Goldhacker et al., 2018 ). But the functional connectivity is calculated based on the amplitude of IMF, not the phase of IMF.
To this end, we introduce a new chronnectome analysis framework in order to estimate phase synchronization based dFNC for group/individual difference identification and behavior and cognitive performance prediction/classification. In this framework, multivariate empirical mode decomposition (MEMD) is employed on time series of ICNs to obtain the characteristic IMFs with physiologically meaningful frequency range. The dFNCs are estimated by computing phase locking values (PLVs) between all pairs of characteristic IMFs using a series of sliding windows. Then k-means clustering is employed on the Fisher transformed dFNC matrices to identify the recurring dFNC states. According to the identified dFNC states, the median dFNC matrix per state is calculated at individual level for investigating the sleep quality relevant group/individual differences. To verify the feasibility of the proposed framework, the dynamic features selected from the median matrices are used for prediction and classification of sleep quality.
In the following sections, we first review the conventional framework for estimating static and dynamic FNC based on Pearson correlation between brain networks, and then introduce a new framework for estimating PLV-based dFNC in Section 2 . Two methods introduced including multivariate empirical mode decomposition and phase locking value estimation are demonstrated in this section. We then validate our approach by assessing the dynamic FNC applying to HCP dataset in Section 3 . Some discussions and concluding remarks are given in Sections 4 and 5 .

Data acquisition and preprocessing
The resting state fMRI (rs-fMRI) data used in this work are collected as part of the Human Connectome Project (HCP) 1200 subject release. The data acquisitions are performed on a customized 3T Siemens scanner (Connectome Skyra) using a standard 32-channel head coil and a body transmission coil. Resting state data are obtained across two ses-sions on consecutive days with a total of four 14.4-minute runs. All data are sampled at TR = 720 ms resulting in 1200 timepoints per run. Detailed descriptions of inclusion and exclusion criteria on the dataset can be found elsewhere Glasser et al., 2016 ).
In this work, we use the time series released with the extensively preprocessed HCP "PTN " (Parcellation + Timeseries + Netmats) dataset to estimate the functional network connectivity at individual level. The HCP "PTN " data are publicly available ( https://db.humanconnectome.org/data/projects/HCP_1200 ).
The details of preprocessing steps can be found in . In brief, a minimal preprocessed pipeline and ICA + FIX ( Griffanti et al., 2014 ) are adopted to remove the potential artifacts in each run of rs-fMRI data. Thereafter, Group-ICA parcellation is performed on the preprocessed fMRI data temporally concatenated across subjects. For each given parcellation (Group-ICA decomposition), the rs-fMRI time series per ICA component at the individual level is derived by mapping the ICA spatial maps onto each subject's rs-fMRI data. Among the HCP PTN dataset, the data from the earliest 184 subjects are reconstructed by an initial version of data reconstruction software, while the data from the latest 812 healthy subjects use an improved reconstruction software to reduce spatial blurring in fMRI. The reconstruction software makes a notable signature on the data that can cause a large difference in fMRI data analysis. Therefore, our study limits the analyses to the latest 812 healthy subjects (ages 22-35 years-old, 410 females) whose rs-fMRI data are reconstructed by the later improved version of software. In addition to the rs-fMRI data of the 812 subjects, we adopt Pittsburgh Sleep Quality Index (PSQI) score as the selected behavioral measure to investigate its relationship with human functional connectome. The scores of PSQI of the 812 healthy subjects range from 0 to 19, where a score of 5 or more is considered as poor sleep quality ( Buysse et al., 1989 ). These data are derived from self-report answers to the Pittsburgh Sleep Quality Index, a 24-item questionnaire comprising 7 component scores.
In the HCP PTN dataset, Group-ICA is applied to the preprocessed resting state fMRI images to obtain group-average parcellations at several different dimensionalities (15,25,50,100,200,300). Time series corresponding to each ICN are extracted for each subject. Previous publication has reported the full details regarding the Group-ICA parcellation, and time series acquisition ( Smith et al., 2015 ). We choose the time series of the first 14.4-minute run from 100 independent components (ICs) for analysis. The number of ICs is chosen as 100 for this study as it represented a balance between a necessary parcellation of major neural networks and the avoidance of the extreme parcellation of visual and cerebellar brain areas ( Nomi et al., 2017 ). Out of the 100 ICs, 51 ICs are characterized as components of resting state networks (RSNs) after discarding those related to movement, imaging artifact or those exhibited peak activations in white matter ( Cordes et al., 2000 ). Based on the anatomical and presumed functional properties, previous rs-fMRI studies Rashid et al., 2016 ), and the information on the Neurosynth ( https://neurosynth.org/ ), the selected 51 RSNs are grouped into 7 functional networks : subcortical (SC), auditory (AUD), somatomotor (SM), visual (VIS), cognitive control (CC), default-mode (DM), and cerebellar (CB). The spatial maps of the 51 selected RSNs grouped by functional domains along with their IC number and peak activation coordinates (x, y and z) can be seen in Inline Supplementary Fig.S1.

FNC estimation based on Pearson correlation coefficient
In conventional analysis, Pearson correlation coefficients estimated over entire time series are adopted as the measure of static FNC. In this work, by using a 5th order Butterworth filter, subject-specific time series are filtered into a frequency range of 0.01 − 0.1 Hz, consistent with other studies reported so far ( Fox and Raichle 2007 ;Wang et al., 2015 ;Wehrle et al., 2018 ). Pairwise correlations are computed between each pair of filtered RSN time series, resulting in a symmetric 51 × 51 correlation matrix for each subject. Notably, Butterworth filter is a type of infinite impulse response (IIR) filter, which has a nonlinearphase distortion, especially in the transition regions of the filters. Despite the inevitable phase distortion in Butterworth filter, extensive researches show that the Pearson correlation calculated from the Butterworth filtered time series is a solid FNC estimation ( Jafri et al., 2008 ;Damaraju et al., 2014 ;Zhi et al., 2018 ).
For the time-varying FNC, previous works used windowed Pearson correlation moving along the time axis. Pre-processing for dFNC analysis requires the same bandpass filtering as static FNC analysis. We use a tapered window, created by convolving a rectangle (window size = 30 TRs = 21.6 s) with a Gaussian ( = 3 TRs) to localize the dataset at each time point. We slide the window in steps of 1 TR, resulting in total W = 1170 windows. For each subject = 1 , 2 , … , 812 , we calculate the Pearson covariance matrices ∑ ( ) , = 1 , 2 , …,1170, from windowed segments of time series to measure the dFNC. Because the windowed data may not have sufficient information to characterize the full covariance matrix, we use graphical LASSO method to estimate the covari- ( Fu et al., 2019 ). The graphical LASSO method induces a shrunken estimator of the precision matrix by maximizing the penalized log-likelihood with an L1 penalty on the precision matrix for sparse estimation. In this method, if the edge linking nodes (brain regions) j and k is absent, then nodes j and k are conditionally independent given the other nodes, and the corresponding entry of the precision matrix is zero. The regularization parameter lambda ( ) is optimized for each subject by using a crossvalidation framework. The dFNC estimates of each window for each subject, 1 ∑ ( ) , = 1 , 2 , …,1170, are concatenated to a 51 × 51 × 1170 array, representing the changes in brain connectivity between components as a function of time. Both static and dynamic FNC estimates are Fisher transformed to stabilize variance prior to further analysis.

FNC estimation based on phase locking value
In this section, we adopt the phase locking value as the index of phase synchronization to measure the dFNC from windowed segments of IC time series. The bandpass filtering is not performed on the IC time series. We first use multivariate empirical mode decomposition (MEMD) to decompose the IC time series into different Intrinsic Mode Functions (IMFs) for every subject, with the specified frequency band for each IMF. The center frequency bands of the first 6 IMFs (i.e., IMF1, IMF2, IMF3, IMF4, IMF5, and IMF6) are 0.35-0.69 Hz, 0.17-0.35 Hz, 0.08-0.17 Hz, 0.03-0.08 Hz, 0.02-0.04 Hz, and 0.01-0.02 Hz, respectively. Then, we identify the IMF4 components as the characteristic IMFs for each IC time series. These IMFs are related to the coherent fluctuations in the lowfrequency range ( < 0.1 Hz) of the blood oxygenation level-dependent (BOLD) signal. Thereafter, phase locking dynamics between each pair of characteristic IMFs are calculated from their windowed segments. Although there existed relatively few MEMD based phase synchrony studies for fMRI analysis, they have been widely employed in other neurophysiological modalities ( Farahmand et al., 2018 ;Mutlu et al. 2011 ;Rutkowski et al., 2008 ;Sweeney-Reed et al. 2007.

Multivariate empirical mode decomposition
The empirical mode decomposition (EMD) is firstly proposed by Huang et al. as an efficient data-driven approach without the requisite of predefined basis ( Huang et al., 1998 ). It has been demonstrated that EMD is preferable for extracting physiological dynamics without either harmonics or limitations with respect to the time and frequency. The EMD approach decomposes a time series into a finite number of components named as intrinsic mode functions (IMFs) by iteratively conducting the sifting process. The IMFs resulted from the sifting process must satisfy the following conditions: (1) The sample-wise mean of the envelope defined by local maxima and local minima must be zero. (2) The number of zero crossings and extrema of the IMF does not differ by more than 1. Through EMD decomposition, the time series data x ( t ) is separated into IMFs as where ( ) , = 1 , 2 , … , are the IMFs and r ( t ) is the remaining residue. The IMFs v 1 ( t ), v 2 ( t ),…, v n ( t ) include different narrow frequency bands ranging from high to low.
EMD can only handle single channel time series in its original formulation (single channel EMD), which means that multiple time series must be decomposed separately. There are several obstacles limiting its application to processing multichannel data, most significantly with the problem of uniqueness associated with single channel EMD, i.e., IMFs differ between time series in number and identified frequency . Specifically, if we use single channel EMD in rs-fMRI time series analysis, EMD is repeatedly applied to analyze each time series. In this way, it is difficult to establish a consistent pattern of EMD decomposition for all rs-fMRI time series of each subject. In addition, singlechannel EMD decomposition technique not only loses the cross-channel information, but also suffers from the mode mixing problem wherein specific frequency information is not localized .
To extend the algorithm for the computation of IMFs of multiple time series, MEMD algorithm is proposed to align the decomposed components of each time series at the same IMF level ( Rehman and Mandic, 2010 ). This decomposition technique utilizes the cross-channel information and provides high localization of the specific frequency component. A noise-assisted method can be introduced in MEMD algorithm. Noise-assisted MEMD introduced extra channels of multivariate noise. The extra noise channels can serve as references to enable more accurate and stable IMFs, which can also alleviate mode mixing problem . MEMD insures the formation of IMFs of multiple time series in the same frequency range for consistent results. Same frequency range of each IMF level is important to estimate the FNC characterized by coherent fluctuations in specific frequency range (for example, < 0.1 Hz) of the BOLD signal. The Matlab script for the MEMD algorithm is provided by Rehman and Mandic, and can be found at http://www.commsp.ee.ic.ac.uk/~mandic/research/emd.htm . The script is improved so that the input data can have any number of time series, instead of the 36 established in the original Matlab code ( https: //github.com/mariogrune/MEMD-Python-). After applying MEMD to the 51 IC time series of each subject, the Fourier spectrum of each IMF level for every subject is estimated to detect the characteristic IMFs with dominant frequency range of 0-0.1 Hz. The extracted characteristic IMFs ( ) , = 1 , 2 , … , 51 corresponding to the 51 IC time series of each subject are used for PLV estimation at individual level.

Phase locking value estimation
The characteristic IMFs of each subject obtained by MEMD satisfy the narrowband requirement of Hilbert transform. Thus, every characteristic IMF ( ) , = 1 , 2 , … , 51 for each subject can have its Hilbert analytic signal as where Here ̃ ( ) represents the Hilbert transform of c k ( t ), the symbol P in Eq. (3) indicates the Cauchy principal value, a k ( t ) is the envelope of the instantaneous amplitude and k ( t ) is the instantaneous phase.
Then, given c m ( t ) and c n ( t ) the characteristic IMFs over two IC time series numbered m and n, m ( t ) and n ( t ) are their corresponding instantaneous phases. The phase differences between the characteristic IMFs are found to fluctuate around a constant value along time. An elegant way of significance testing, known as the phase locking value (PLV), provides a reliable measurement of phase synchrony. The PLV within the time window is defined as where T is the window size or data length for PLV calculation. In the fully synchronized time series, the phase difference is a constant and PLV = 1. If the time series are unsynchronized, the phase difference follows a uniform distribution and PLV = 0. A simple way to calculate the PLV is using normalized Hilbert analytic signal. According to Eq. (2) , the normalized Hilbert analytic signal of c m ( t ) is So the PLV can be calculated as where the symbol ( · ) * stands for complex conjugate, ̄ is the vector version of ̄ ( ) , and ̄ is the conjugate transpose of ̄ .
In this study, we use a rectangle window (window size = 30 TRs = 21.6 s) to calculate the PLV at each time point. We also slide the window in steps of 1 TR, resulting in total W = 1170 windows. For each subject = 1 … 812 , we calculate the PLV matrices ( PLV ) i ( w ), = 1 , 2 , …,1170, from windowed segments of time series to measure the dFNC. The PLV-based dFNC estimates of each window for each subject are concatenated to form a 51 × 51 × 1170 array, representing the changes in brain connectivity between components as a function of time.
We also estimate the PLV-based static FNCs by computing PLVs between pairs of IMFs over the total duration of the recording for each subject. Both static and dynamic FC estimates are Fisher transformed to stabilize variance prior to further analysis.

Clustering analysis
To assess the frequency and recurring patterns of FNC, k-means clustering is applied to the windowed FNC matrices of all subjects. The kmeans algorithm uses the L1 distance function (Manhattan distance) based on our previous dFNC work . The optimal number of dFNC states is determined by the elbow criterion of the cluster validity index, computed as the ratio between within-cluster distance to between-cluster distance.
To ensure that all the identified states are present in all subjects, we add an additional k-means clustering step at subject level : the resulting centroids (cluster medians) from group level k-means clustering are used to initialize the clustering of the windowed FNC matrices for each subject. A similar idea on using clustering but with Principal Component Analysis at fMRI scanning session level can be found in Ponce-Alvarez et al's work ( Ponce-Alvarez et al., 2015 ).

Prediction ability for sleep quality
To test whether PLV-based dFNC features can better represent sleep quality associated with sleep quality, we conduct two studies: continuous phenotype prediction and classification analysis. Note that sleep quality in the following experiments is represented by the PSQI scores.

Continuous phenotype prediction
In our experiments, we compare the prediction performance of the proposed PLV-based dFNCs and conventional FNCs. We use a 10-fold cross-validation (CV) technique to evaluate the prediction performance of all these FNCs. That is, the whole set of subjects is first randomly partitioned into 10 disjoint subsets of about equal size; then each subset is successively selected as the test set and the other 9 subsets are used for training the predictive model; across all subjects in the training set, a linear regression of Pearson correlation is conducted to relate each edge in the FNC matrices to the PSQI score, and the most relevant edges is selected for use in the predictive model; and finally the trained linear regression model is applied to predict PSQI scores of the subjects in the test set. This process is repeated for 100 times independently to reduce the effect of sampling bias in the CV. The prediction procedure closely follows the similar metrics used elsewhere ( Finn et al., 2015 ;Shen et al., 2017 ). The predictive power is calculated as the correlation values between the predictive and observed PSQI scores. The significance of the predictive power is assessed by the permutation testing, which is used to generate an empirical null distribution of the test statistic. Test statistics are calculated under random rearrangements of the labels for the data. Specifically, permutation is done by preserving the structure of the FNC matrices but randomly reassigning PSQI scores (e.g., subject 1 ′ s FNC matrix is paired with subject 2 ′ s PSQI score). For predictions based on dFNC, the features selected from one state or multiple states are separately used for predicting the PSQI scores. As each dFNC state may contain different information for prediction, incorporating functional connectivity information from a spectrum of dFNC states into a single predictive model may yield best performance. Finally, predicted PSQI scores are generated for each FNC profile (static or dynamic) or their combination.

Classification analysis
We investigated whether we can classify two subsets according to sleep quality. To perform the experiment, we define bad and good sleep quality groups according to PSQI scores. The 151 subjects (80 females) with PSQI scores ranging from 0 to 2 are adopted for the group of good sleep quality, while the 162 subjects (84 females) with PSQI scores ranging from 7 to 19 are selected for the group of bad sleep quality. Then we implement the feature selection procedure as demonstrated in the section of continuous phenotype prediction. Due to the nature of crossvalidation, a slightly different set of edges in FNC is selected as features in each iteration. In this work, the edges that appear in at least 90% across the iterations are selected as the features for classification.
A support vector machine (SVM) approach is used for classification ( Cortes and Vapnik, 1995 ). Additionally, we intend to determine whether applying PLV-based dFNC can improve the ability of discrimination relative to using static FNC profiles. Thus, features extracted from static and PLV-based dFNC are separately used for classification. Each classification experiment is repeated 100 times. For each run, all subjects are randomly partitioned into 10 equal size subsamples. Of the 10 subsamples, a single subsample is retained as the validation data for testing, and the remaining 9 subsamples are used as training data. This cross-validation process is repeated 10 times, with each of the 10 subsamples used exactly once as the validation data. A standard SVM function built in Matlab with Gaussian kernel is applied, and a grid search is performed to optimize parameters within the SVM (e.g., the radius of Gaussian kernel, the weight of the soft margin cost function).

dFNC based on Pearson correlations between the band-passed time series
Static FNC is estimated as the Pearson correlation coefficient (PCC) between each pair of IC time series for intrinsic networks. We separately average static FNC matrices for subjects in the group of bad and good sleep quality to obtain the group-averaged FNC matrix. The results are presented in Fig. 1 (a). A two-sample t -test is performed to examine significant ( q < 0.05, corrected for false discovery rate (FDR) and multiple comparisons) sleep quality-related differences in the strength of static connections among intrinsic networks. The between-group differences are visualized through plotting the log of p value with the sign of t statistics, -sign ( t ) log 10( p ). In general, we observe that the relationships between SC and AUD/VIS networks are stronger in subjects with poor sleep quality, while the relationships between SC and CB, VIS and AUD are weaker in subjects with good sleep quality ( Fig. 1 (b)).
The PCC-based dFNCs are estimated using graphical LASSO method, wherein we compute Pearson correlation coefficient matrices from a series of windowed segments. The optimal number of dFNC states is determined to be 4. Dynamic states for the two groups, subjects with bad and good sleep quality, are shown in the first two rows of Fig. 2 . Among the four states, state 3 can be characterized as weakly connected due to the weak connectivity intensity among functional networks, while the others can be characterized as connected. As shown by the percentage on the top of Fig. 2 , the weakly connected state is detected more frequently than the connected ones. Specifically, state 3 accounts for 47.35% and 44.62% of the total time period for subjects with bad and good sleep quality, respectively. Note that not all subjects have dynamic windows that are assigned to every state, thus the count of observations varies with state. For instance, the number of subjects in states 2 and 4 is below the total number of each group (162 and 151 for subjects with bad and good sleep quality, respectively).
To investigate the between-group differences in PCC-based dFNC, an element-wise subject median is computed for every state using each subject's windowed PCC matrices. The median matrix has been considered as the representative pattern of the subject for a particular state ( Jafri et al., 2008 ). Similar to the between-group difference study of static FNC, a two-sample t -test with the significance level q < 0.05 (corrected for FDR and multiple comparisons) is performed on the median matrices for each state across the two groups. The significant betweengroup differences in functional networks among these four identified states are depicted in the last row of Fig. 2 . Overall, the general betweengroup differences in static FNC are found in state 2 and 4 of the dFNC. But the dFNCs show considerable variation in the significant sleep quality-related connectivity differences among the four brain states. As shown in Fig. 2 , no sleep quality-related differences are observed in the weakly connected state (state 3), and only a few differences in functional networks are obtained in state 1 across the two groups. By contrast, a large number of significant sleep quality-related differences are detected in states 2 and 4. The details are discussed as follows. In general, for subjects with bad sleep quality, the inter-network relationships between SM/VIS and SC/CC/DM/CB are reported to have increased connectivity in state 4 in contrast to those with good sleep quality. At the same time, there are weaker connections within individual network systems of SC, SM, and VIS. Besides, weaker connections can also be found between networks SM and VIS in state 4. State 2 shares some similar between-group differences as state 4: increased connectivity between SM and SC/CB; decreased connectivity within individual network systems of SC and SM. But the connectivity between SM and CC/DM, and between VIS and DM/CB is weaker in state 2 for subjects with bad sleep quality. Another difference in state 2 is the increased connectivity in individual network system of VIS. Other stronger relationships can be found in inter-network relationships between CB and AUD/CC/DM, and in intra-network of DM in state 2.
According to the above results, PCC-based dFNC provides more informative measures to detect sleep quality-related group differences than static FNC. However, the states including most distinguishing features (states 2 and 4) are not available for all subjects, which prevents the application of using PCC-based dFNC features per state for prediction or classification of individual's sleep quality. Previous researches on capturing dynamic brain states based on PCC-based dFNC features encountered the same problem, i.e., not all subjects have dynamic windows  Fig. 2. Sleep-quality-related differences in PCC-based dynamic FNC. PCC-based dFNC across 4 brain states was estimated with k-means clustering of windowed FNC matrices for the subjects with bad sleep quality (top) and good sleep quality (middle). FNC is inverse Fisher transformed for display. Total percentage of occurrence for each state and total number of subjects in each state are listed above each picture. Significant sleep quality-related differences in FNC between all networks are computed for each of the 4 states using with 2-sample t-tests at a significance level of q < 0.05, corrected for false discovery rate and multiple comparisons. Effects are color-coded for bad sleep quality -good sleep quality computations, where the hotter the color is, the worse the sleep quality is. that are assigned to every state Liu et al., 2017 ;Cai et al., 2018 ;Vergara et al., 2018 ;Zhi et al., 2018 ;Espinoza et al., 2019c ).
With subject-level clustering, four identified states and the betweengroup differences among these states can be seen in Inline Supplementary Fig. S2. These estimated states are similar to those estimated from the original framework as shown in Fig. 2 . But the total percentage of occurrence for each state and total number of subjects in each state vary. Similar but fewer between-group differences in functional networks are detected in states 2 and 4 in comparison with the results in Fig. 2 .

PLV-based dFNC estimation
The MEMD method is adopted to decompose the ICA time series of 51 networks into different IMF levels which have common frequency components. For each subject, the ICA time series from a 14.4-minute run period (1200 timepoints) is composed of the multivariate time series (51 × 1200). Then, the MEMD method is used to decompose the multivariate signal (51 × 1200) of each subject into a new multivariate time series (51 × 1200 × M ) with M = 10 or 11 IMF levels. The different amount of IMF levels between subjects is reflected by the lowest frequency IMFs and the residual, which is related to the linear and curvilinear trend in the time series. The Fourier frequency spectrum of IMFs with 10 or 11 levels across all subjects can be found in Inline Supplementary Fig. S3. The IMF4 components are consistent with the physiologically meaningful frequency range of spontaneous oscillations of the brain. Hence, in this section, the time-varying FNCs are estimated by computing PLVs between pairs of IMF4 components of the 51 networks using a series of sliding windows. We also estimate the PLV-based static FNCs for further comparison, and the results can be found in Inline Supplementary Fig. S4. We apply the k-means clustering algorithm to the windowed PLV matrices and identified 4 dynamic patterns that reoccurred in time and across subjects. PLV-based dynamic states for the two groups, subjects with bad and good sleep quality, are shown in the first two rows of Fig. 3 . Among the four states, state 1 shows the most intense phase synchronization across the global brain, while state 3 is a weakly connected state due to the low phase synchronization between the functional networks. State 2 can be characterized as low connected because most of the PLV values in this state are relatively low. State 4 exhibits strong phase synchronization within and between somatomotor, visual and cognitive control networks. As shown by the percentage on the top of Fig. 3 , the weakly connected and low connected states are detected more frequently than the other states. As the count of subjects in each state listed above each picture of Fig. 3 , all subjects have dynamic windows that are assigned to the 4 states.
To investigate the between-group differences in PLV-based dFNC, subject median matrices are computed from every subject's windowed PLV matrices per state. Then, we examine significant ( q < 0.05, corrected for FDR and multiple comparisons) sleep quality-related differences in the strength of phase synchronization among intrinsic networks for each state. The results are depicted in the last row of Fig. 3 . It is evident that no sleep quality-related differences survived the FDR correction of q < 0.05 for state 3. Decreased phase synchronization in state 1 can be observed only in inter-network relationships between SM and VIS, and in intra-network of SC for subjects with bad sleep quality. By contrast, more sleep quality-related differences are significant in state 2 and 4. As shown in state 4, for subjects with bad sleep quality, decreased phase synchronization can be found in inter-networks between SC and AUD/SM/VIS/DM/CB, AUD and SM/VIS/DM/CB, and between SM and VIS/CB. The weaker phase synchronization can also be found in the intra-networks of SC, SM, VIS, and CB. Increased connectivity can be observed within individual network system of CC, and in inter-networks between CC and SM/VIS/DM/CB, DM and SM. State 2 displays some opposite results. For subjects with bad sleep quality, decreased phase synchronization can be observed in intra-network of CC, and in internetworks between CC and SM/VIS/DM/CB. Increased connectivity can be found within individual network system of SC. We also attempt to add an additional subject level k-means clustering step to the framework for PLV-based dFNC estimation. By comparing the results between the modified and previous approaches, we find that PLV-based clusters resulted from the modified framework are almost identical to the original ones.
At the condition of k = 4 clustered states, the detailed information about the number of dynamic windows for each state of the two groups is presented in Fig. 4 . The minimum numbers of dynamic windows for both groups are in the state with the most intense phase synchronization (state 1), which are 32 and 56 for subjects with bad and good sleep quality, respectively. We also try to cluster the windowed PLV matrices across subjects in each group into 5 and 6 states, the results can be found in Inline Supplementary Fig. S5. It is found that the minimum numbers of dynamic windows are still in the state with the most intense phase synchronization. Under the conditional of k = 5 clustered states, the minimum number of dynamic windows is 11 and 4 for subjects with bad and good sleep quality, respectively. The minimum count of dynamic windows reduces to 8 and 1 for subjects with bad and good sleep quality at the condition of k = 6 clustered states.
For a fair comparison, we also estimate the dFNC based on Pearson correlation derived from the amplitude of IMF4 components, instead of the band-passed time series. We consider two cases: with and without including subject level k-means clustering into the framework of dFNC analysis. For the first case, only k-means clustering at group level is applied to the windowed FNC matrices to assess the recurring patterns of FNC. Four identified states and the between-group differences in functional networks among these states can be seen in Inline Supplementary Fig. S6. The estimated states are similar to those based on Pearson correlations between the band-passed time series. But the clustering results are different, as reflected by the total percentage of occurrence for each state and total number of subjects in each state. Significant sleep quality-related differences are detected in states 2 and 4. But with only group level clustering, not all the identified states are present in all subjects. For the second case, k-means clustering at subject level is added after the group level clustering in the framework. Four identified states and the between-group differences in functional networks among these states can be seen in Inline Supplementary Fig. S7. The estimated states are also similar to those based on Pearson correlations between the band-passed time series. Significant sleep quality-related differences are detected in states 2 and 4. In this case, all the identified 4 states are present in all subjects. Hence the dynamic features of each state can be used for prediction and classification of sleep quality.

Prediction of connectivity for sleep quality
In this section, we explore whether PLV-based dFNC profiles can better describe the characteristics of sleep represented by the total score across all items on the Pittsburgh Sleep Quality Index (PSQI). Specifically, we firstly test the prediction ability of the PLV-based dFNC profiles for sleep quality. Secondly, we investigate whether we could train a classifier to discriminate between two subsets based on differing sleep quality by applying PLV-based dFNC.

Continuous phenotype
We first teste whether PLV-based dFNC within the whole brain can better predict sleep quality for novel subjects. A 10-fold cross-validation approach is repeated 100 times to evaluate the sleep quality prediction performance. Prediction scores applying both PCC-based static FNC and PLV-based dFNC are displayed in Fig. 5 . As shown in the figure, the correlations between predicted and observed PSQI scores estimated by PLV-based dFNC of both state 2 and 4 are significant better than that estimated by PCC and PLV-based static FNC, wherein the prediction using PLV-based dFNC of state 4 achieves the best correlation coefficient of 0.312 ± 0.016 among the four states. PLV-based dFNC of state 1 and 3 has similar prediction performance as the PCC-based static FNC. These  Fig. 3. Sleep-quality-related differences in PLV-based dynamic FNC. PLV-based dynamic functional connectivity (FC) states for the subjects with good sleep quality (top) and bad sleep quality(middle). FNC is inverse Fisher transformed for display. Total percentage of occurrence for each state and total number of subjects in each state are listed above each picture. Significant sleep-quality-related differences in FNC between all networks are computed for each of the 4 states using 2-sample t-tests at a significance level of q < 0.05, corrected for false discovery rate and multiple comparisons. Effects are color-coded for bad sleep quality -good sleep quality computations, a hotter color indicates worse sleep quality. results are consistent with the t -test results of Fig. 3 , i.e., the state providing more informative features to detect between-group differences achieved better prediction performance of sleep quality. By combining the PLV-based dFNCs of state 2 and 4, the correlation coefficient between predicted and observed PSQI scores is further improved to 0.364 ± 0.019. However, the prediction performance is reduced when combing more brain states for PSQI score prediction. The statistical significance for the correlation between predicted and observed PSQI scores is assessed by performing 1000 nonparametric permutation tests for every scenario. The results demonstrate that the prediction of PSQI scores is above chance for scenarios including predictions in terms of PLV-based dFNC according to state 2, state 4, the combination of state 2 and 4, the combination of state 1, 2, and 4, and all four states. The other scenarios including predictions in terms of PCC and PLV-based static FNC, and PLV-based dFNC state 1 and state 3 yield weaker predictions across results that are not statistically significant under permutation testing.
We then evaluate the sleep quality prediction performance according to the dFNCs based on Pearson correlations between the band-passed time series. The correlation coefficients between the predicted and observed PSQI scores estimated by dFNC state 1, state 2, state 3 and state 4 are 0.067 ± 0.021, 0.143 ± 0.022, 0.018 ± 0.023, and 0.131 ± 0.023, respectively. The dFNCs derived from Pearson correlations between the IMF4 components are also used for sleep quality prediction. The corre-lation coefficients between the predicted and observed PSQI scores estimated by dFNC state 1, state 2, state 3 and state 4 are 0.096 ± 0.023, 0.151 ± 0.021, 0.040 ± 0.024, and 0.193 ± 0.018, respectively. The sleep quality prediction results for these two types of PCC-based dFNCs can be found in Inline Supplementary Fig. S8.

Classification analysis
For further analysis, we investigate whether we could discriminate between two subsets according to the sleep quality. The features (edges) selected from static or dynamic FNC for classification are shown in Fig. 6 . For each scenario of Fig. 6 , the selected features are those most relevant to the sleep quality. As we can see from Fig. 6 , more features in both static and dynamic FNCs are selected for classification than those showing significant between-group difference across all subjects. This is particularly noticeable for PCC and PLV-based static FNC and PLVbased dFNC state 1 and 3, in which few or no features survived the FDR correction in estimating the between-group differences associated with sleep quality.
The comparisons of the classification are illustrated in Fig. 7 . The PCC and PLV-based static FNC approach show an average classification accuracy of 61.19% and 61.93%, respectively. Based on the selected features in the PLV-based dFNC, the accuracy increases to 64.50, 76.88%, 70.35%, and 85.07% for state 1, state 2, state 3 and state 4, respectively. This is consistent with the results in Fig. 3 . State 4 has the largest   5. The prediction performance based on the selected features in PCC-based static FNC, PLV-based static FNC, and PLV-based dFNC profiles including state 1, state 2, state 3, state 4, the combination of state 2 and 4, the combination of state1, 2, and 4, and all states. The white dot in the middle is the median value and the thick black bar in the center represents the interquartile range.
between-group differences, as detected by the significant t -tests, providing evidence to prefer this state for classification. State 2 followed State 4. A significant increase is observed, boosting accuracy from 85.07% to 92.61%when combining states 2 and 4 for classification. A marginal but significant increase in accuracy is achieved when combining states 1, 2 and 4 for classification (93.87%). The accuracy can be further improved when including all states in dynamic FNCs for classification (94.17%). Combining the PLV-based static and dynamic FNC features results in an average accuracy of 94.53%, but it does not significantly improve the classification over using the PLV-based dynamic FNC features alone.
We also evaluate the classification results between bad and good sleep quality groups according to the two types of PCC-based dFNCs. A 10-fold cross-validation classification is repeated 100 times for both types of dFNCs. The results can be seen in Inline Supplementary Fig. S9. Based on the selected features in the PCC-based dFNC between pairs of the band-passed time series, the classification accuracies are 67.23%, 70.33%, 63.14%, and 69.72% for state 1, state 2, state 3 and state 4, respectively. The accuracy is improved when including all states of dynamic FNCs for classification (76.97%). Based on the selected features in the PCC-based dFNC between pairs of IMF4 components, the classification accuracies are 69.93%, 74.20%, 65.32%, and 78.16% for state 1, state 2, state 3 and state 4, respectively. The accuracy is improved to 86.11% when including all states of dynamic FNCs for classification.

Discussion
In this work, we present a new framework to evaluate the dFNC based on phase locking value instead of conventional Pearson correlation coefficient. FC/FNC estimation based on Pearson correlation have three drawbacks: First, if two signals are 90°out of phase, they will result in a Pearson correlation of 0, even though these signals are highly locked ( Gravel et al., 2018 ). PLV does not suffer from these drawbacks and the PLV measure is convenient to interpret. Second, there is an ongoing debate on the underlying mechanisms of negative functional connectivity and how to deal with the negative correlations ( Chen et al., 2011 ;Buckner et al., 2013 ;Gravel et al., 2018 ;Parente et al., 2018 ;Pedersen et al., 2018 ;Zhu et al., 2018 ). A few studies excluded negative functional connectivity-related results to avoid uncertainty. In addition, the negative functional connectivity can lead to non-stationarity for PCC-based dFNCs. This kind of nonstationarity is detected when strong positive Pearson correlations of edges temporally transit towards strong negative correlations, or vice  versa ( Pedersen et al., 2018 ). Such temporal transitions are not detected in PLV-based dFNCs because strong negative and positive Pearson correlations will both have high phase synchrony. The above drawbacks of PCC will lead to a large variance of PCC-based dFNC states, which probably results in the inconsistency of recurring dFNC states across individuals.
According to the strength of brain functional integration and segregation, brain connectivity states for a subject can be roughly divided into three types . When the integration is low, there is high segregation, and nodes (brain regions) are weakly connected and behave more independently. By contrast, when integration is high, segregation is low. The optimal working point of a brain network occurs when segregation and integration are balanced . Con-verging evidence indicates that phase synchronization is a basic mechanism for brain functional integration across the whole brain over time ( Roskies, 1999 ;Varela et al., 2001 ). There is a large field of research of phase synchrony in MEG/EEG/ECoG analysis, where PLV has been used to much larger extent than in fMRI research ( Aydore et al., 2013 ;Fries 2015 ;Palva et al., 2012 ;Singer, 1999 ). PLV-based dFNC states are consistent with the three types of brain functional integrations. In our work, PLV-based dFNC state 1 is the state with strongest integration exhibiting most intense phase synchronization across the global brain, while state 3 is the state with the weakest integration. The functional integration and segregation strengths of states 2 and 4 are balanced, so states 2 and 4 are most likely the states having the optimal working points of the brain networks. The group difference between subjects with bad and good sleep quality also exists in these two states. On the contrary, it is difficult to map the PCC-based dNFC states to the three types of brain functional integrations. A large part is because of the negative functional connectivity in PCC-based dNFC, the interpretation of which has encountered with criticism. We find that even when the PCC-based dFNC states are clustered into 3 states, not all the identified states are present in all subjects. The 3 PCC-based dynamic states for subjects with bad and good sleep quality can be seen in Inline Supplementary Fig. S10. Therefore, PLV-based dFNCs are probably more appropriate for identifying the recurring states. In addition to the above described characteristics of PLV-based dFNC, there are other two important requirements to ensure that all the identified PLV-based dFNC states are present in all subjects, i.e., higher temporal resolution and enough scanning time of fMRI. These two requirements can provide the basic guarantee that all recurring states can be collected during fMRI signal acquisition. One of the big drawbacks of PLV measures is the bias discussed in detail in Aydore et al. ( Aydore et al., 2013 ). The bias is inversely proportional to the number of samples (N) and is negligible for N > 50. It also reduces as true PLV increases for fixed N. In this work, PLVs are estimated by using the average within a time-window with 30 samples. The PLV values are higher than 0.5 for states 1, 2, and 4, and higher than 0.4 for the state 3 (see Fig. 3 ). According to the work of Aydore et al., the bias of the estimated PLV is very low for states 1, 2, and 4, and low for state 3 in our study.
Until recently, there is few research on characterizing individual difference associated with dFNC state. Rashid et al. proposed a method using regression coefficients against PCC-based dFNC states for the classification of schizophrenia and bipolar patients, which is the first trial of using individual difference from dFNC states to improve the classification performance ( Rashid et al., 2016 ). But as the authors stated, the differences in pairwise correlation across different states and different groups are not used as dynamic features, rather the states are used as regression matrix while obtaining regression coefficients. The regression coefficients are secondary measures and not good for straightforward interpretation of brain functional organization associated with behaviors. It is hard to know the contribution of each specific FNC feature (edge) to behavior prediction from the regression coefficients. Followed by Finn et al. and Shen et al.'s work ( Finn et al., 2015 ;Shen et al., 2017 ), the prediction model with PLV-based dFNC can project the features back into the brain space, which can facilitate the interpretation of relationships between brain functional organization and behaviors, in addition to the comparison with existing literatures.
We have followed Rashid et al.'s method ( Rashid et al., 2016 ) by regressing out the dFNC states. Three types of dFNC is considered: (a). dFNC based on Pearson correlations between pairs of band-passed time series; (b). dFNC based on Pearson correlations between pairs of IMF4 components; (c). dFNC based on Phase Locking values between pairs of IMF4 components. A 10-fold cross-validation classification experiment is repeated 100 times for each of the three cases. The average classification accuracies for the three cases are 71.34%, 77.84%, and 70.21%, respectively. As shown by the results, all classification accuracies are not high. The PLV-based dFNC for classification got a comparable perfor-mance as the amplitude-based dFNC when using the regression method by Rashid et al.. As shown in many previous rs-fMRI studies, PCC-based dFNC states resulted from k-means clustering always have a great similarity. For example, in Espinoza et al.'s work ( Espinoza et al., 2019a ), they identified 5 dFNC states and found that three of them had similar patterns except for the connectivity strength within and/or between resting state networks. This kind of similarity between PCC-based dFNC states can be found elsewhere Damaraju et al., 2014 ;Cai et al., 2018 ;de Lacy et al. 2019 ;Espinoza et al., 2019b ;Espinoza et al., 2019c ). Such phenomenon can also be found in our work: PCC-based dFNC states 1, 2 and 4 show similarities in connectivity within the 7 functional networks and between some of them (e.g., state 2 is similar to state 4 with similar connectivity between SC networks and AUD/SM/VIS networks). On the contrary, as shown in Fig. 3 , the similarity between the 4 PLVbased dFNC states is relatively low. When we cluster all PLV-based dFNC matrices into 5, 6, or 7 states, the similarity between states is still low (These results can be seen in Inline Supplementary Fig. S11). Notably, when the PLV-based dFNC states are clustered into up to 6 states, all subjects still have dynamic windows that are assigned to every dFNC state. Even under the condition of 7 clustered states, only a few subjects lack one of the dFNC states. The PLV-based dFNC matrices provide values for each individual and show promise in identifying brain states in rs-fMRI.
In this work, we find decreased phase synchronization within subcortical network and in inter-networks between subcortical network and other networks from the PLV-based dFNC state 4 for subjects with bad sleep quality. The subcortical network is involved in the mediating of sleep ( Miyamoto et al., 2012 ;Abbott et al. 2013 ), which controls sleepwake transitions and rapid eye movement (REM)/non-REM sleep alternation within sleep ( Pace-Schott and Hobson 2002 ). It has been revealed that the brain regions belonging to subcortical network are susceptible to worsen sleep quality ( Li et al., 2020 ). Previous research has demonstrated that functional connectivities in subcortical network are significantly declined in subjects with sleep deprivation, and its connections with other brain regions are also decreased ( Zhou et al., 2017 ). These findings are in line with our discovery. We also find decreased phase synchronization in the intra-networks of somatomotor and visual network from the PLV-based dFNC state 4 for subjects with bad sleep quality. This is in accordance with previous sleep -related FC studies ( Picchioni et al., 2014 ;Killgore et al., 2015 ), which evidences that sleep deprivation deteriorated alertness and processing of sensory or visual information.
The other findings are related to the cognitive control network. Firstly, we find decreased phase synchronization in the intra-network of cognitive control and inter-networks between cognitive control network and other networks from the PLV-based dFNC state 2 for subjects with bad sleep quality. Previous researches have demonstrated that PSQIdefined poor sleepers showed alterations in the right superior frontal gyrus, the hippocampus, and the insula. These three regions are implicated in cognitive and emotional functioning ( Scullin 2017 ). It has been identified that the sleep quality indices (PSQI scores) in chronic primary insomnia patients are negatively correlated with reduced resting state functional connectivity between the left orbitofrontal cortex and the left superior temporal gyrus, between the right basal ganglia and the right superior frontal gyrus ( Zhou et al., 2016 ). The orbitofrontal cortex and, superior temporal gyrus are involved in cognitive process. Basal ganglia is associated with functions like cognition and emotion, while superior frontal gyrus is involved in self-awareness. Another research finds that sleep deprivation leads to reduced anticorrelated activity between the default-mode and cognitive control network ( Farahani et al., 2019 ). These findings are consistent with our discovery.
Secondly, we find increased phase synchronization in the intranetwork of cognitive control and subcortical, and in inter-networks between cognitive control network and other networks, for subjects with bad sleep quality. This phenomenon may reflect a compensatory effort to maintain normal brain function and alertness after poor sleep. This kind of compensatory mechanism has been found in many studies ( Varkevisser et al., 2007 ;Cruz Gomez et al., 2013 ;Gui et al., 2015 ;Zhu et al., 2016 ;Hartwigsen, 2018 ). For example, sleep deprivation triggers compensatory mechanisms to overcome fatigue that are manifested by enhanced functional connectivity in specific brain regions, including the thalamus, paracentral lobule, supplementary motor area, postcentral gyrus and lingual gyrus ( Zhu et al., 2016 ). The mechanisms supporting compensatory flexibility of cognitive networks can be found elsewhere ( Hartwigsen, 2018 ).

Conclusion
In this work we have proposed a novel framework to estimate phase locking value based dFNC for group/individual difference identification, and further for behavior and cognitive performance prediction/classification. Through such a framework we can distinguish PLVbased connectivity patterns recurring in entire fMRI session, which are consistent across individuals and groups. These connectivity patterns have different level of global synchronization which are consistent with previous findings of brain functional integration . We also find that the PLV-based dFNC provides a useful tool for both prediction and classification of individual's sleep quality. As shown by the results, the predictions based on certain PLV-based dFNC states and their combination can achieve significantly better results than that based on PCC or PLV-based static FNC, and PCC-based dFNC derived from the band-passed time series or IMF components. The classification results further illustrate the effectiveness of the proposed method. Classification using dynamic features from PLV-based dFNC outperforms those using PCC/PLV-based static FNC and PCC-based dFNCs. In addition, combining PLV-based static and dynamic features does not significantly improve the classification performance over dFNC features alone, suggesting that static FNC does not add new information when combined with dFNC for classification purposes. A similar conclusion about combining static and dFNC features for classification can be found in Rashid et al.'s work ( Rashid et al., 2016 ). In summary, the proposed framework contribute a valuable way to study dFNCs with added value by incorporating phase information, leading to improved classification and prediction of population groups with different behaviors and cognitions.

Ethics statement
Both the original ethical review as well as any further ethical review and data use agreement pertinent to the submitted manuscript.

Declaration of competing interest
None.