Cortical electrophysiological evidence for individual-specific temporal organization of brain functional networks.

The human brain has been demonstrated to rapidly and continuously form and dissolve networks on a subsecond timescale, offering effective and essential substrates for cognitive processes. Understanding how the dynamic organization of brain functional networks on a subsecond level varies across individuals is, therefore, of great interest for personalized neuroscience. However, it remains unclear whether features of such rapid network organization are reliably unique and stable in single subjects and, therefore, can be used in characterizing individual networks. Here, we used two sets of 5-min magnetoencephalography (MEG) resting data from 39 healthy subjects over two consecutive days and modeled the spontaneous brain activity as recurring networks fast shifting between each other in a coordinated manner. MEG cortical maps were obtained through source reconstruction using the beamformer method and subjects' temporal structure of recurring networks was obtained via the Hidden Markov Model. Individual organization of dynamic brain activity was quantified with the features of the network-switching pattern (i.e., transition probability and mean interval time) and the time-allocation mode (i.e., fractional occupancy and mean lifetime). Using these features, we were able to identify subjects from the group with significant accuracies (~40% on average in 0.5-48 Hz). Notably, the default mode network displayed a distinguishable pattern, being the least frequently visited network with the longest duration for each visit. Together, we provide initial evidence suggesting that the rapid dynamic temporal organization of brain networks achieved in electrophysiology is intrinsic and subject specific.

INTRODUCTION

One fundamental goal of neuroscience is to link brain and behavior. The uniqueness of one's mind could be mainly attributed to its connectome shaped by nature and nurture (Mueller et al., 2013; Seung, 2012). Despite gross similarities revealed by population‐average inferences, the whole‐brain functional connectome at rest displays high interindividual variability and has been repeatedly found to be predictive of individual differences across many cognitive domains (Beaty et al., 2018; Mueller et al., 2013; Tavor et al., 2016; Wang et al., 2015). With the development of personalized neuroscience, assessing the stability and reproducibility of brain profiles becomes a central issue. There is evidence that the individual connectome of the brain is as unique to the individual as their fingerprint (Finn et al., 2015; Horien, Shen, Scheinost, & Constable, 2019; Kong et al., 2019). However, these studies were mostly based on static functional connectivity (sFC), which assumes connections as spatially and temporally stationary over a measurement period.

Recently, an emerging body of literature focusing on the concept of dynamic functional connectivity (dFC) has strongly suggested that functional connections are constantly changing over time (Allen et al., 2014; de Pasquale et al., 2010; de Pasquale, Corbetta, Betti, & Della Penna, 2018; de Pasquale, Della Penna, Sporns, Romani, & Corbetta, 2016; Hutchison et al., 2013; Vidaurre, Smith, & Woolrich, 2017). Moreover, taking advantage of the high temporal resolution of non‐invasive electrophysiological recording such as magnetoencephalography (MEG), accumulating evidence has suggested that the human brain dynamically forms and dissolves networks on a time scale of seconds (Brookes et al., 2014; O'Neill et al., 2015). Recent findings have demonstrated that such network modulation can be observed even on a millisecond level (Baker et al., 2014; Vidaurre et al., 2018), with several specific brain networks repeatedly occurring at different time points and typically lasting ~100 ms. These fast temporal organization properties provide effective substrates and flexibility to allow the human brain to efficiently support different cognitive processes (Bressler & Menon, 2010; Bressler & Tognoli, 2006), and also enable our brains to quickly adapt to rapidly changing environments (Beggs, 2008; Matthews & Hampshire, 2016). Importantly, interindividual variability in this organization of dynamic brain networks has been found to be linked with cognition, demographics, and the presence of psychiatric illnesses (Fatima, Kovacevic, Misic, & McIntosh, 2016; Liu, Liao, Xia, & He, 2018; Lopez, Pusil, Pereda, Maestu, & Barcelo, 2019; O'Neill et al., 2017; Pedersen, Zalesky, Omidvarnia, & Jackson, 2018; Reinen et al., 2018). For instance, a recent MEG study focusing on developmental dynamic networks found that such fast network modulation underlies age‐related changes in functional connectivity (FC) and is potentially implicated with childhood disorders (Brookes et al., 2018); moreover, network dynamic analysis can provide measures quantifying features of neural synchrony disorder in Alzheimer's disease (Sitnikova, Hughes, Ahlfors, Woolrich, & Salat, 2018), or enable the tracking or prediction of postoperative cognitive decline (Carbo et al., 2017). However, despite the importance of understanding the relationships between temporal organization properties of brain function and behavior, it remains unclear whether features of such rapid network organization are as reliably unique and stable in single subjects.

Given that functional networks are suggested to be dominated by common organizational principles and stable individual features with minor day‐to‐day variability (Gratton et al., 2018; Satterthwaite, Xia, & Bassett, 2018) and that temporal dynamics is correlated with individual behaviors, we hypothesized that features of individuals' temporal organization of brain networks are stable across different sessions and can be used to identify subjects from a large group. Provided that under resting state both microscale neuronal populations (Berkes, Orbán, Lengyel, & Fiser, 2011) and macroscale networks (Vidaurre et al., 2017) are organized in a temporally coordinated way, we modeled the functional connectomic dynamics of spontaneous brain activity as recurring networks fast shifting between each other in a coordinated manner. Specifically, we analyzed two 5‐min MEG resting‐state datasets from 39 healthy subjects over two consecutive days. MEG measures extracranial fields induced by currents in the brain, and mathematical modeling of these fields allows reconstruction of neuronal electrical activities on a millisecond time scale. Our study used the Hidden Markov Model (HMM) to model temporal organization of large‐scale brain networks, namely, the temporal characteristics and interactive coordination of discrete and repeating large‐scale networks (or states) over time (Baker et al., 2014; Vidaurre et al., 2017). In detail, individual temporal organization of functional networks was quantified with features of the network‐switching pattern (i.e., transition probability, mean interval time) and the time‐allocation mode (i.e., fractional occupancy, mean lifetime). In addition, to describe the global time‐varying profiles of brain network organization, each subject's static and dynamic FC profiles on each day were also calculated respectively, following the data analyses in previous functional magnetic resonance imaging (fMRI) studies (Finn et al., 2015; Liu et al., 2018). Time‐resolved statistics of network temporal organization and static and dynamic FC profiles were used to identify subjects from the group. We further probed the differences between temporal characteristics of recurring brain networks.

MATERIALS AND METHODS

Data acquisition

Fifty‐eight healthy subjects participated in this study, all of whom were recruited from local universities. Fifty‐one subjects who did not undergo electrooculography (EOG) recording during MEG experiments were used in the main identification analyses, while the remaining seven subjects (19–25 years old, mean age 20.1 years; three males) underwent EOG recordings and were used to investigate whether subject‐specific eye movement patterns drive identification accuracy. The procedures of this study were approved by the Peking University Institution Review Board, and all subjects provided written informed consent. MEG data were recorded using a whole‐head Triux system (Elekta Neuromag, Helsinki, Finland) at a sample frequency of 1,000 Hz in a magnetically shielded room. Each subject was scanned with MEG over two consecutive days, and each session consisted of a 2‐min eyes‐open resting‐state recording, followed by a 5‐min eyes‐closed resting‐state recording. Only the eyes‐closed data were used for the present analysis. A wooden head supporter was used to help minimize head movements during scanning. The head position relative to the sensor array was recorded continuously by four head position indicator (HPI) coils. To allow coregistration of the MEG data with brain anatomy, three anatomical landmarks (nasion, left and right preauricular points), HPI coils and the scalp surfaces (~200 points) were digitized using a 3D digitizer (Polhemus, VT) at the start of each scan session. After MEG acquisition, T1‐weighted anatomical magnetic resonance images (MRI) were acquired using a FSPGR sequence (repetition time [TR] = 6,676 ms; echo time [TE] = 2.92 ms; flip angle=12; voxel size 1 × 1 × 1 mm) on a 3T GE Discovery MR750 MRI scanner (GE Healthcare, Milwaukee, WI). Coregistration of MEG data to anatomical MRI was then implemented by matching the digitized head shape to the head shape extracted from the MRI.

MEG preprocessing and source reconstruction

MEG data were initially inspected visually using Brainstorm (Tadel, Baillet, Mosher, Pantazis, & Leahy, 2011). MEG channels that exhibited excessive noise (mainly caused by the operating condition of hardware) were first identified and removed (removed number: mean 4.0 ± SD 0.7). Considering that segmented data may not completely reflect the actual dynamic nature, only full 5‐min episodes were included for analyses. We manually filtered data of all subjects so that, if either of the two data episodes collected across two days from the same subject displayed obvious artifacts, this subject would be excluded. This filtering process eventually excluded 12 subjects. In the remaining 39 subjects (18–32 years old, mean age 22.6 years; 20 males), MEG data were preprocessed to remove the external interference using the temporal extension of Signal Space Separation (tSSS) method (Taulu & Hari, 2009) with the MaxFilter 2.2 software (Elekta‐Neuromag, Helsinki, Finland), using a subspace correlation limit of 0.9 and a sliding window of 10 s.

Following tSSS, the data were bandpass filtered into 0.5–48 Hz band, and then an atlas‐based beamforming approach (Hillebrand, Barnes, Bosboom, Berendse, & Stam, 2012) was used to perform source reconstruction implemented in FieldTrip (Oostenveld, Fries, Maris, & Schoffelen, 2011). The cortex of each subject was parcellated into 78 regions of interest (ROIs) using the Automated Anatomical Labelling (AAL) atlas (Tzourio‐Mazoyer et al., 2002). This was done by spatially normalizing each subject's anatomical MRI to the MNI template using a nonlinear transformation implemented with SPM12 (https://www.fil.ion.ucl.ac.uk/spm/software/spm12) and labeling all voxels according to the AAL atlas. Then the inverse registration to individual anatomical space was performed. For the forward solution, a single shell model (Nolte, 2003) was created for the head model and each ROI was divided into numbers of voxels with side lengths of 5 mm. For beamforming, data covariance was defined within a 0.5–48 Hz frequency range and a time window spanning the whole 5‐min eyes‐closed experiment. Dipole orientation was determined through a nonlinear search for an optimum signal‐to‐noise ratio. The beamformer weights were normalized by its vector norm and then multiplied with the sensor data matrices to obtain the source timeseries.

The source timeseries were then filtered into five classical MEG/EEG subfrequency bands (delta [0.5–4 Hz], theta [4–8 Hz], alpha [8–13 Hz], beta [13–30 Hz], and lower gamma [30–48 Hz]) using a basic finite impulse‐response (FIR) filter in EEGLAB (https://sccn.ucsd.edu/eeglab/index.php). Each ROI contained a specific number of voxels, all with their own source timeseries. The amplitude envelopes of these timeseries were calculated by Hilbert transform for all frequency bands (i.e., the broadband and five subfrequency bands), and these envelopes were then averaged across all voxels in a particular ROI to obtain a representative timeseries. Therefore, for each subject in each scan session, a set of 78 timeseries for each frequency band was finally obtained, and then this set of timeseries was downsampled to 100 Hz for computational efficiency. Notably, a potential confound in beamformer‐derived source timeseries is that two signals (e.g., from two regions) may exhibit spurious correlations purely as a result of source leakage (spatially blurred representation of the source distribution) (Brookes, Woolrich, & Barnes, 2012; Colclough et al., 2016; Palva & Palva, 2012). Therefore, a multivariate symmetric orthogonalization (Colclough, Brookes, Smith, & Woolrich, 2015) was employed on the source timeseries of all bands. In addition, for subsequent HMM analysis, the timeseries for each subject in each session were demeaned and normalized to unit variance and then concatenated temporally across all subjects and two sessions into a single group matrix.

Hidden Markov Modeling

Fast transient dynamic networks were identified using an HMM (Baker et al., 2014; Vidaurre et al., 2017; Woolrich et al., 2013) on the group concatenated MEG source timeseries. The HMM assumes that the timeseries data can be described as a hidden sequence of a finite number of discrete states, and at each time point, the brain is in one of these states. The model assumes that the states cannot be imaged directly but can be inferred from measured observables (here, the concatenated timeseries), and that the state being active at time t relies only on which state was active at time t−1 (i.e., Markovian). Each state here is represented by a multivariate Gaussian distribution, defined by the mean and covariance, which represents brain networks of distinct activity and FC, respectively. An HMM with 12 states was inferred, and variational Bayes inference (Rezek I, 2005) was used on the HMM to derive the full posterior distribution on the model parameters (implemented in the HMM‐MAR toolbox; https://github.com/OHBA-analysis/HMM-MAR). Accounting for variations in the inference due to different initializations, five realizations were implemented, and the model with the lowest free energy was chosen. Viterbi decoding was performed to determine which of the derived states the brain was most likely in for each time point (i.e., state time course). Whereas the states were estimated at the group level, the state time course that indicated the dynamic characteristics of brain networks changing with time was obtained for each individual in each session.

The individual state time course indicated the temporal structure of brain network states during a specific experimental period (e.g., one scan session), and four statistics for each subject per session were estimated to quantitatively describe the individual temporal organization characteristics of the inferred states (Baker et al., 2014). The individual transition probability (TP) matrix is the probability of transition from one state to another state. Different from the group‐level TP, which could be calculated from concatenated single‐group data, we calculated the TP matrix from individual actual nonprobabilistic state time courses (i.e., Viterbi decoding path) to indicate the switching pattern of dynamic networks for each subject in each session: where indicates which of the states is active at time t, and TP describes the probability of transitioning from state i to state j between time t−1 and t.

The mean interval time (IT) describes the average time between two adjacent occurrences of each state: where is one if and is zero otherwise, and T is the length of the state time course.

Fractional occupancy (FO) is defined as the proportion of time spent visiting each state:

The mean lifetime (LT) is defined as the average time spent in one state before transitioning to another state:

For these four statistics, we classified TP and IT as the network‐switching pattern, which indicates how states are sequentially ordered and how long it takes for a state to recur; FO and LT belonged to the time‐allocation mode, which delineates the time duration of each visit and the total duration proportion of all visits for different states.

Analysis of static and dynamic FC profiles

To explore the uniqueness of global time‐varying profiles of individual network organization, we also calculated static and dynamic whole‐brain FC profiles similar to previous fMRI studies (Finn et al., 2015; Liu et al., 2018). For the static network construction, Pearson correlation coefficients between the source timeseries of each pair of 78 AAL nodes were calculated. For dynamic network construction, connectivity was estimated as a function of time, with a widely used sliding time‐window approach (Hutchison et al., 2013; O'Neill et al., 2017). The whole 5‐min source timeseries was segmented into overlapping time windows, and each window width was 10 s with a 1 s sliding step. Then, 290 time windows were obtained for each subject in each scan session. For each time window, Pearson correlation was used to calculate the FC between each pair of 78 regions, which means 3,003 pairs of connections (i.e., 78 × 77/2). The FC matrix was vectorized to a 1 × 3,003 row vector, and then the vector was concatenated in time to form a 290 ×3,003 dynamic FC matrix. Therefore, each column of the dynamic FC matrix represented the time course of an individual connection between two AAL regions. The dynamic FC variability and stability profiles were calculated to describe the time‐varying characteristics of each connection. In detail, the dynamic FC variability indicated the overall fluctuating level for each connection over time, and the stability indicated the tendency to remain stable during a short period (see details in Supporting Information). After processing as above, two 1 × 3,003 row vectors were obtained for each subject per scan session.

Identification analysis

To quantify the uniqueness of individual time‐resolved features extracted from the HMM and static and dynamic FC matrices, an individual identification analysis presented by Finn (Finn et al., 2015) was performed. Specifically, for each subject, we compared each statistic (i.e., TP, IT, FO, and LT computed using the HMM, static FC profiles and dynamic FC variability and stability) of this subject from Day 1 to the same statistic of all subjects on Day 2. We computed the similarity using the Pearson correlation coefficient across the vector form of each statistic for each comparison. The predicted identity of this subject was assigned the same label as the subject on Day 2, who was optimally similar to this subject. If the predicted identity matched the true identity, the identification results of this comparison were assigned a score of 1; otherwise, a score of 0 was assigned. Then, Day 1 and Day 2 were reversed, and all the processes described above were performed again. Finally, the identification accuracy rate was defined as the sum of the identification score in the two days divided by twice the total number of subjects. Having obtained the identification accuracy for each statistic, a nonparametric permutation test was performed to compute its statistical significance. In detail, the vector forms of each statistic of all subjects from Day 1 were compared to a randomly permuted list of subjects from Day 2, such that each subject on Day 1 was assigned a “correct” identity corresponding to a different subject on Day 2. The identification processes described above were reperformed and the identification accuracy was recorded. Then, Day 1 and Day 2 were reversed. The permutation test was conducted 10,000 times, and the p value was defined as the proportion of the 10,000 corresponding random identification accuracy values that were equal to or greater than the actual value.

Quantifying the contributions of state temporal characteristics to identification

An analysis quantifying the contributions of specific states to individual identification was performed in a manner consistent with the prior fMRI study (Finn et al., 2015). We first computed a group consistency measure (Φ). Given two row vectors and that represented one of the temporal statistic (i.e., IT, FO, LT) values of the HMM inferred states (based on the broadband source signals) in two sessions after z‐score normalization with each element corresponding to one specific state, the definitions of Φ were as follows: where i indexes subject, k indexes state, N is the total number of subjects (N = 39) and ∑ indicates the sum of over all subjects. A high Φ value for a given state indicated that the state's temporal characteristics were consistent both within a subject and across the group.

To further identify the state that truly contributed to individual identification, we calculated the differential power (DP):where , denotes the consistency of a specific state between different subjects, and | indicates the probability that the consistency value between two different subjects is higher than it is within the same subjects for a specific state. States with larger DP values indicate greater contributions to individual identification.

Factors affecting identification accuracy

We evaluated the reliability of our identification results using different processing procedures or varying some parameters in our analysis. First, we estimated new HMMs from scratch using 6–16 states for the broadband source signals and repeated the identification analysis to examine the effect of different numbers of states on identification accuracy. The average value of the identification accuracies (i.e., the comprehensive index, CI) of the four HMM statistics (TP, IT, FO, LT) was used to uniformly describe the results. Second, for the dFC profile analysis, an important issue is that the optimal choice of time‐window width remains debatable (Hutchison et al., 2013; Leonardi & Van De Ville, 2015; O'Neill et al., 2018). Two different time‐window widths (5 s and 15 s) were used to construct the dynamic FC matrix in the broadband and a similar identification analysis was reperformed.

Moreover, to exclude the possibility that significant identification was caused primarily by subject‐specific head movement patterns, we performed identification based on head motion feature only. Specifically, we calculated framewise displacement (FD) similar to fMRI studies (Ide & Li, 2018; Power, Barnes, Snyder, Schlaggar, & Petersen, 2012) for each subject in each session. For each time point t, the FD was calculated as: , where ( and ( are the translation and rotation transformation parameters; r = 50 mm, a constant that approximates the mean distance from the cerebral cortex to the center of the head and used to covert rotations into displacement; and (other parameters, i.e., , could be expressed in the same way). Furthermore, we computed the mean and SD of the FD across all subjects and two sessions, specified 50 bins spanning the mean ± 3 SD, and obtained head motion distribution vectors accordingly. Then, we used the head motion distribution vectors alone to perform identification analysis. Similarly, to investigate whether successful identification was simply driven by subject‐specific eye movement patterns, we also calculated EOG distribution vectors using the seven subjects for whom EOGs were recorded. For each subject in each session, the EOG value at each time point was divided by the mean EOG values of that session, and then the mean and SD of EOG values across all subjects and sessions were calculated. EOG distribution vectors were estimated and used to perform identification in the same way as head motion.

RESULTS

Nonrandom dynamic switching between brain networks is subject specific

An HMM with 12 states was inferred from the resting‐state MEG data of 39 subjects over 2 days. Each subject had a characteristic state time course that manifested as different states being active at particular time points. Figure 1a shows an illustrative 2,000‐ms section of the state time course for one example subject based on the broadband source signals, and mean activation maps of all 12 HMM states in the broadband are shown in Figure 1b (thresholded mean activation maps are shown in Figure S1). The spatial maps of the HMM inferred states were interpreted through combining fMRI research literature with the spatial correlation between the activation pattern of each state and a widely‐used seven‐network template (Yeo et al., 2011), and further these states were divided into two groups (see details in Supporting Information). One group is composed of higher‐level cognitive regions, such as the DMN (State 1; Raichle et al., 2001), language network (State 2; Friederici, Chomsky, Berwick, Moro, & Bolhuis, 2017), decision‐making areas (State 3; Wikenheiser & Schoenbaum, 2016), and salience network (State 4; Seeley et al., 2007; Sridharan, Levitin, & Menon, 2008). State 5 appeared to be whole‐brain negative, indicating that the whole‐brain activities were lower than the average activity level (deactivation state). Meanwhile, the other group consists of sensory and motor areas, such as bilateral auditory network (States 6 and 7), visual network (States 8–10), and somatomotor network (States 11 and 12).

Figure 1

Twelve HMM network states based on the broadband (0.5–48 Hz) source signals. (a) Example 2,000‐ms section of the state time course for a subject. (b) Mean activation maps of all 12 HMM network states in the broadband. The first five states correspond to the higher‐level cognitive group, and States 6–12 correspond to the sensory‐motor group

Using each subject's temporal structure of the HMM network states, we illustrated that the switching pattern between brain networks is not a stochastic feature but rather an innately representative feature of one's brain connectivity. In particular, we extracted the TP (Figure 2a) and IT (Figure 2b) to quantify individual network‐switching patterns. For the TP, we constructed two TP matrices for each subject corresponding to the two scan days; as shown in Figure 2a, the matrices of the same subject on the 2 days had considerable similarity (R = .95), while the similarity of the matrices on the diagonal was relatively lower (R = .62 or .72). All subjects' TP matrices and IT values for 12 networks were used to identify subjects from the group, with identification accuracies in the broadband of 37.2 and 35.9%, respectively (Figure 3a). In specific frequency bands, the identification accuracies of TP and IT were 15.4 and 14.1% in the delta band, 20.5 and 17.9% in the theta band, 37.2 and 30.8% in the alpha band, 42.3 and 34.6% in the beta band, and 11.5 and 11.5% in the gamma band, respectively (Figure 3a). All identification accuracies were significantly higher than in the random identification case (all p < .0001 except p < .05 for TP and IT in the delta and gamma bands).

Figure 2

Schematic diagrams of the network temporal organization features inferred from the HMM in the broadband. (a) Schematic diagram of transition probability matrices from two example subjects (rows) for two consecutive days (columns). The matrix shows the probabilities from one state transitioning to any other different one. R: Pearson correlation coefficient between two specific matrices. (b–d) Schematic diagrams illustrating the mean interval time (IT), fractional occupancy (FO) and mean life time (LT) of each state from six example subjects for two consecutive days. Different colors denote the 12 states, and the length of each color block illustrates the relative time length (of IT and LT) or proportion (of FO) for a certain state

Figure 3

Identification accuracies of statistics for each frequency band. (a) Identification accuracies of the network temporal organization features inferred from the HMM. Transition probability (TP) and mean interval time (IT) indicate the network‐switching pattern between the 12 network states, and fractional occupancy (FO) and mean life time (LT) indicate the time‐allocation mode between the 12 inferred states. (b) Identification accuracies of static and dynamic functional connectivity (constructed with the sliding time‐window method) profiles. The asterisks denote that the accuracy is significantly higher than in the random identification case (permutation test; **p < .0001, *p < .05). Static: static functional connectivity; var: dynamic functional connectivity variability; sta: dynamic functional connectivity stability. Broadband: 0.5–48 Hz; δ: 0.5–4 Hz; θ: 4–8 Hz; α: 8–13 Hz; β: 13–30 Hz; γ: 30–48 Hz

Nonrandom time allocation for distinct brain networks is subject specific

Two statistics from the temporal structure of the HMM network states were used to quantitatively evaluate the time‐allocation mode: FO (Figure 2c) and LT (Figure 2d). Our results demonstrated that properties of the time‐allocation mode of recurring brain networks can successfully recognize individuals, meaning the time‐allocation mode is subject specific. Specifically, the FOs of 12 HMM networks for the same subject were very similar, whereas different subjects' FOs were obviously distinctive. The identification accuracies of FO and LT were 52.6% and 33.3%, respectively, in the broadband (Figure 3a). The identification accuracies in specific frequency bands for FO and LT were 20.5 and 10.2% in the delta band, 23.1 and 14.1% in the theta band, 43.6 and 32.0% in the alpha band, 30.8 and 37.2% in the beta band, and 17.9 and 12.8% in the gamma band, respectively (Figure 3a). All the identification accuracies were significantly higher than in the random identification case (all p < .0001 except p < .05 for LT in the delta, theta, and gamma bands).

MEG static and dynamic FC profiles identified subjects with significant accuracy

The static and dynamic FC profiles captured by MEG identified subjects from the group with significant accuracy. The identification accuracies obtained from static and dynamic FC profiles are summarized in Figure 3b. For sFC, the identification accuracy in the broadband was 74.4%, and the identification accuracies in subfrequency bands were 30.8% in the delta band, 34.6% in the theta band, 46.2% in the alpha band, 70.5% in the beta band, and 67.9% in the gamma band. For dFC variability and stability, the identification accuracies in the broadband were 50.0 and 47.4%, respectively, while the identification accuracies in the subfrequency bands were as follows: delta band, 19.2 and 12.8%; theta band, 21.8 and 17.9%; alpha band, 26.9 and 21.8%; beta band, 32.0 and 14.1%; and gamma band, 12.8 and 10.2%, respectively. All of the identification accuracies above were significantly higher than in the random identification case (all p < .0001 except p < .05 for variability in the gamma band and stability in the delta, beta, and gamma bands).

DMN displayed patterns that were obviously distinguishable from other networks

Our results also substantiated the hypothetical role of the DMN as a special network during resting state among these transiently switching networks. Figure 4a illustrates all subjects' IT, FO and LT of the 12 inferred HMM states across two scan sessions for the broadband data. For the network‐switching pattern, the mean amount of time spent between two consecutive visits to a particular state (IT) was between 700–5,000 ms on average; as for the time‐allocation mode, these states occupied 2.9–10.5% of the total time (FO) and were stable for 70–130 ms on average (LT). In terms of all states, both the IT and LT of the DMN (State 1) were the longest (IT: 4910 ms, LT: 130 ms) and occupied the lowest percentage (2.9%) of the total duration (FO); namely, while the DMN was much less frequently visited under resting‐state conditions than the other states, each occurrence of the DMN was on average longer than that of the other network states. In addition to the DMN, within the higher‐level cognitive group of networks (States 1–4 and deactivation State 5 included in this group), the deactivation state had the largest FO (10.5%) and a relatively long LT (100 ms). In contrast, the sensory‐motor group of networks (States 6–12) displayed subtler variations than the higher‐level cognitive group.

Figure 4

Temporal characteristics and contributions to individual identification of the HMM states based on the broadband source signals. (a) Boxplots of temporal characteristics of the 12 network states. The mean interval time (IT), fractional occupancy (FO) and mean life time (LT) for each state are shown across all 39 subjects and two scan sessions. The dashed line divides 12 states into the higher‐level cognitive group (left side) and the sensory‐motor group (right side). The asterisks in each diagram denote that the corresponding value of a state differs significantly from all the other states across the two groups (paired t‐test; **p < .0001, *p < .05). (b) Differential power and group consistency (Φ) of network states in terms of the IT, FO, and LT. The group consistency was transformed to natural logarithmic values (ln Φ) for demonstration

Temporal characteristics of the DMN were more subject general than those of other networks

To explore the extent to which different states contribute to individual identification, as shown in Figure 4b, we calculated the DP and the group consistency (Φ) for each state's temporal statistics (IT, FO, and LT) in the broadband. For DP, the DMN had the lowest DP values for FO and LT and a relatively low value for IT, whereas the deactivation state had the highest DP values for FO and relatively high values for IT and LT. Regarding the group consistency (Φ), across the three temporal statistics, the DMN had the highest Φ value, and the deactivation state also showed a relatively high value. In brief, the temporal characteristics of the DMN were more subject general, while the characteristics of the deactivation state were more subject specific when compared to those of the other resting‐state networks.

We explored several factors that may affect identification accuracy. First, for the HMM analysis, we varied the number of estimated HMM states from 6 to 16 for broadband source timeseries. Apart from the number of states being six for LT, the identification accuracies of all number of states were significantly higher than in the random identification case (Figure S2; all p < .0001 except p < .05 for number of states being seven for LT). In addition, the identification accuracy of the CI tended to increase with the number of states, probably because when the division of states became more fine‐grained, individual differences in temporal organization were better reflected. Second, different time‐window widths (5 s and 15 s) were used to construct the dFC matrix in the broadband, and the identification results are shown in Figure S3, which were all significantly higher than those in the random identification case (all p < .0001). Moreover, to investigate whether significant identification was simply driven by subject‐specific head motion or eye movement patterns, we tested whether subjects could be identified based on the distribution of their head framewise motion or eye movement. The identification rates on the basis of head motion and eye movement were 5.1 (p > .25) and 14.3% (p > .16), respectively. These low accuracies indicated that our significant identification was unlikely due to subject‐specific motion patterns in the scanner. Together, these results demonstrated the reliability of our identification results.

DISCUSSION

Taking advantage of the high temporal resolution of MEG, we show for the first time that the temporal organization of large‐scale networks is reliably subject‐specific on a subsecond time scale. We demonstrated that network temporal organization, which was quantitatively analyzed using systematic time‐resolved features extracted from the temporal structure of the HMM states and global time‐varying profiles of network dynamics, displayed stable interindividual variations and can be used to identify individual from the group with significant accuracy. Our results suggest that personal uniqueness is imprinted on the large‐scale network's rapid temporal organization, which lays the prerequisite for investigations into how individual‐specific temporal dynamics of brain networks give rise to different cognitive processes.

Among all the statistics used in our study, the four HMM inferred statistics (TP, IT, FO, and LT) and the dynamic variability and stability were the measurements that describe the temporal dynamic features of brain networks (i.e., dynamic measurements), while the static FC profile captured the average connectome level during a dynamic process (i.e., static connectivity). We found that static connectivity identified subjects with a much higher accuracy than dynamic measurements. On the one hand, the lower identification associated with the dynamic measurements could be a result of ongoing mind wandering activities (Kucyi, 2018; Kucyi & Davis, 2014), in that different ongoing thoughts may exist in two scan sessions even for the same individual and that the dynamic measurements reflect more such ongoing states than the static one. On the other hand, it is also possible that the limitations of current methodologies failed to generate time‐evolving connectivity that could fully capture the core traits of dFC. For instance, one major consideration about the HMM is the temporally mutual exclusivity (i.e., only one state is active at each time point), which is not necessarily a physiologically meaningful feature of the brain. However, the HMM states may be thought of as representing the most dominant state at each given time point, and multiple states rendering simultaneously could be represented through the FO of different states at different time scales. Based on the dynamic structure obtained by the HMM, our results indicate that the temporal organization features of networks could be used to identify subjects from the group with significant accuracy.

In addition, we found that different subfrequency bands showed distinct identification accuracies. Specifically, alpha and beta bands displayed the highest accuracies for the four HMM statistics (TP, IT, FO, LT) and tended to have relatively high accuracies for dFC variability and stability. Previous studies have revealed that brain rhythms are heritable phenotypes, which are under strong genetic control, and different genetic origins of brain rhythms result in different heritabilities (Buzsaki, Logothetis, & Singer, 2013). Meanwhile, several EEG studies have held that alpha and beta bands of brain FC at rest exhibit higher heritability than other frequency bands (Smit et al., 2010; Smit, Stam, Posthuma, Boomsma, & de Geus, 2008). Based on these findings, our result that alpha and beta bands displayed higher recognition accuracy might highlight the interesting relationship between alpha and beta band FC and genetics, yet no direct genetic evidence is available here. Future studies may be conducted to probe into this issue and to further investigate whether alpha and beta bands are superior candidate frequency bands in developing dFC features as possible brain measurements. Together, these findings pave the way for thoroughly understanding the neuronal mechanisms behind the emergence and representations of individual brain connectome differences.

Furthermore, we linked the differences between the temporal characteristics of 12 recurring HMM networks and their known functions established at the group level. Particularly, a vast body of literature has reached a consensus that the DMN is a hallmark of spontaneous fluctuations and plays a special role in brain functions, instantiating processes that support emotional processing, self‐referential mental activity, and the recollection of prior experiences (Fox & Raichle, 2007; Raichle, 2015). Current results show that the DMN was the least frequently visited state and the one with the longest duration for each visit; also, in accordance with previous studies, the DMN occupied only a small part of the total scan duration (Baker et al., 2014; Vidaurre et al., 2018). Such a standout pattern significantly distinguished it from other transient networks, and we speculate that this pattern is in line with its special role during resting state (Raichle, 2015). In addition, the deactivation state, wherein the whole‐brain activities were lower than the average level, also demonstrates an intriguing pattern because it has the highest FO among all states and was more subject specific compared to other states. We hold that the whole‐brain deactivation phenomenon may indicate that the brain is truly “at rest” in this state. Taken together, our results potentially contribute to a better understanding of the temporal characteristics of network organization under resting state.

Notably, in our study, the brain parcellations of individual subjects were derived from group‐level brain atlases. Although this parcellation process is a common approach, recent literature has shown that the interindividual variability in brain anatomy is not negligible; for example, the position, size, and shape of brain regions can vary between individuals (Amunts et al., 2005; Li et al., 2019). Future studies can build on individual brain anatomy or functional segmentation, which may lead to improved identification accuracy due to its ability to capture the idiosyncrasies of subjects' functional network construction.

Based on the general idea that the spurious correlation caused by source leakage manifests as a zero‐lag linear coupling, solutions for leakage reduction generally remove any zero‐lag connection between sources. These solutions might be oversimplified for individual‐level studies, considering that zero‐lag connections may contain FC information vital to identification, such as the zero‐lag synchronization of neural oscillations (Gollo, Mirasso, Sporns, & Breakspear, 2014; Roelfsema, Engel, Konig, & Singer, 1997). At the core, source‐leakage effects are rooted in MEG source reconstruction issues. In our study, we used beamformer for inverse modeling of our data, as beamformer fits resting‐state source localization for its interference rejection properties (Brookes et al., 2011; Hillebrand, Singh, Holliday, Furlong, & Barnes, 2005; Van Veen, van Drongelen, Yuchtman, & Suzuki, 1997). However, this approach typically yields a spatially blurred representation of the source distribution referred to as source leakage. Future studies should be conducted to specify the role of zero‐lag synchrony in whole‐brain interaction and to advance the methodological solutions for disambiguating source leakage effects from true zero‐lag interaction. As using standard beamformer source reconstruction can hardly solve these problems (Colclough et al., 2016), future solutions may be improved in terms of using methods for signal reconstructions that can suppress contributions from other cortical sources (Hauk & Stenroos, 2014; Wens et al., 2015) and approaches based on regularization priors that could obtain spatially sparse source distribution (Huang et al., 2014; Uutela, Hamalainen, & Somersalo, 1999).

In conclusion, our study provides cortical electrophysiological evidence that the individual rapid temporal organization of brain networks is reliably unique and stable on a subsecond level, and also that the differences between temporal characteristics of recurring brain networks are related to their known functions established on a group level. Future studies should be conducted to establish precise associations between the temporal organization features of networks and subjective experience and objective behavior. Also, similar studies can be improved by recording EOG in all subjects and then removing eye‐related artifacts. Meanwhile, because recent research has proposed that resting state may not be the best experimental condition for investigating interindividual differences (Greene, Gao, Scheinost, & Constable, 2018), the more appropriate conditions under which brain‐behavior associations across different domains can be better explored should be identified. We expect that our findings will contribute to the developing field of personalized neuroscience and facilitate investigations into the possible neuromarkers of cognition, behavior and diseases.

CONFLICT OF INTERESTS

The authors have no financial or competing interests to declare.

Appendix S1. Supporting Information.

Click here for additional data file.