Dynamic changes in brain lateralization correlate with human cognitive performance

Hemispheric lateralization constitutes a core architectural principle of human brain organization underlying cognition, often argued to represent a stable, trait-like feature. However, emerging evidence underlines the inherently dynamic nature of brain networks, in which time-resolved alterations in functional lateralization remain uncharted. Integrating dynamic network approaches with the concept of hemispheric laterality, we map the spatiotemporal architecture of whole-brain lateralization in a large sample of high-quality resting-state fMRI data (N = 991, Human Connectome Project). We reveal distinct laterality dynamics across lower-order sensorimotor systems and higher-order associative networks. Specifically, we expose 2 aspects of the laterality dynamics: laterality fluctuations (LF), defined as the standard deviation of laterality time series, and laterality reversal (LR), referring to the number of zero crossings in laterality time series. These 2 measures are associated with moderate and extreme changes in laterality over time, respectively. While LF depict positive association with language function and cognitive flexibility, LR shows a negative association with the same cognitive abilities. These opposing interactions indicate a dynamic balance between intra and interhemispheric communication, i.e., segregation and integration of information across hemispheres. Furthermore, in their time-resolved laterality index, the default mode and language networks correlate negatively with visual/sensorimotor and attention networks, which are linked to better cognitive abilities. Finally, the laterality dynamics are associated with functional connectivity changes of higher-order brain networks and correlate with regional metabolism and structural connectivity. Our results provide insights into the adaptive nature of the lateralized brain and new perspectives for future studies of human cognition, genetics, and brain disorders.

Introduction based laterality index that can effectively capture brain lateralization characteristics underlying higher-order cognition [14], defined as the difference between an ROI's BOLD correlation with the GS of left brain and its correlation with the right brain at each time window (see Methods). GS-based laterality index of an ROI reflects whether the activity of the ROI is more synchronized with the left or the right hemisphere. The correlation between BOLD signal of an ROI and the left/right hemispheric global signal (GS L/ GS R ) represents its synchronization with the left/right hemisphere, respectively. Higher correlation of an ROI with GS L compared to GS R indicates left-hemispheric lateralization and vice versa. DLI of an ROI (a time series of laterality index) is then obtained by calculating the laterality index for each time window; see Fig 1A. Positive DLI of a region indicates stronger interaction with the left hemisphere (i.e., leftward laterality), while a negative one indicates rightward laterality.
We demonstrate that the GS-based laterality index we adopted is highly similar to a conventional laterality measure autonomy index (AI) [12,33] that is widely used. Using resting-state fMRI data from HCP, we found high correlation coefficient between the whole-brain laterality profile obtained by the mean laterality index (MLI, average across all time windows) and that obtained by AI (Pearson r = 0.64 ± 0.12, all participants have p < 0.05). The MLI has good replicability over 4 sessions involving left and right scans (correlation across sessions: 0.66 ± 0.018); see S1 Fig. Furthermore, GS-based laterality index is efficient in that it only takes 2 � n time, compared to the traditional ROI-based method that takes n 2 time (n being the number of regions), therefore is economic for multiple time windows.
We further characterized the dynamic changes in laterality of a region from 2 different perspectives: the magnitude and the sign of laterality, by LF and LR, respectively; see Fig 1D. LF is defined as the standard deviation of laterality time series. Standard deviation of a time-varying measure is commonly used in dynamic brain network analysis (e.g., DFC), which reflects its variability [24,[35][36][37]. LR specifically refers to the number of zero crossings of laterality (switch between left and right laterality) in 2 consecutive windows. It is also inspired by the dynamic brain network analysis, i.e., a region may change the module it belongs to (named "nodal flexibility" proposed by Bassett and colleagues), which correlated closely with cognitive ability [38,39].
We illustrate laterality time series and these 2 measures using a brain region with high level of mean laterality (0.26, left laterality; Fig 1D). The standard deviation of time-varying laterality index corresponds to moderate changes, or deviations, with respect to the mean laterality. By contrast, LR reflects larger changes in laterality, corresponding to sufficiently extreme deviations from the mean, which results in sign changes in laterality. Moreover, a large magnitude of the laterality index corresponds to segregation at the level of hemispheres, and a small magnitude (near 0) indicates integration across 2 hemispheres; see Fig 1D and "Network and structural basis of laterality dynamics" (Results) for explanation.
Based on the above 2 dynamic laterality characteristics, we found that the laterality of most brain regions varied considerably over time. Specifically, primary visual (LF = 0.46 ± 0.09;

Spatial clustering of the laterality dynamics
After characterizing the temporal LF across the whole brain, we investigated how spatially distributed regions that support different functional specializations correlate with each other in laterality ( Fig 1B). Through spatial clustering of the laterality time series across all brain regions, we identified 4 major clusters ( Fig 3A): Cluster 1 consisted mainly of regions from the bilateral visual network; Cluster 2 consisted of regions from the bilateral sensorimotor and cingulo-opercular networks; Cluster 3 mainly covered the frontoparietal network (more from the right hemisphere) and attention network; and Cluster 4 mainly consisted of regions from bilateral default mode network, part of the frontoparietal network (mainly the left hemisphere) and regions from the language network ( Fig 3A). Clusters 1 and 2 showed higher levels of LF (reflected by higher standard deviation of DLI and LR) when compared to Clusters 3 and 4 (repeated-measures ANOVA, p < 0.001, pairwise t test with p < 0.01; see S3 Fig).
We then explored the association between laterality dynamics and out-of-scanner cognitive performance. Specifically, we used multiple tasks adopted by HCP that are generally associated with either lateralized (e.g., language: story comprehension task, and attention: Flanker task) or bilateral (e.g., cognitive flexibility: Card Sort task, and working memory: List Sorting task) brain function. First, we found that the LF and LR of the identified spatial clusters correlated in a positive and negative manner, respectively, with cognitive performance. Second, significant correlation was only found for language function (story comprehension tasks), cognitive flexibility (Card Sort Test and pattern comparison task), and processing speed (Pattern Comparison Test) out of all 13 cognitive measures; see S1 Table. Correction for multiple comparisons were performed using the false discovery rate (FDR) method for all 29 laterality measures employed in our analyses (see Methods section for details of these laterality measures). Furthermore, we also adopted a more conservative multiple-comparisons correction scheme for the 13 behavioral measures. That is, based on the FDR correction of the 29 dynamic laterality measures (FDR q < 0.05), we further employed Bonferroni correction for the 13 behavioral measures, i.e., FDR q < 0.05/13 = 0.0038. All the following results were based on the FDR correction for all 29 laterality measures (FDR q < 0.05). Those correlations that passed the stricter correction method (FDR q < 0.05/13) were also marked.
We further investigated the behavioral correlates of individual variations in the state transition properties. Results showed a link between State 2 and cognitive functions: Individuals with less State 2 showed higher language task difficulty (fraction, r = −0.1, p = 0.007; dwelling time, r = −0.09, p = 0.015), higher CardSort performance (fraction, r = −0.08, p = 0.014) and higher ProcSpeed score (fraction, r = −0.09, p = 0.01). In addition, CardSort was positively correlated with the fraction (r = 0.1, p = 0.006) and dwelling time (r = 0.1, p = 0.014) of State 1. See Fig 4C and 4D and S3 Table. Multiple comparisons correction was performed using the FDR method for all laterality measures used in the analyses; see Methods section.

Network and structural basis of laterality dynamics
Next, we investigated how dynamic laterality is related to functional and structural brain network properties. First, we investigated the functional correlates of dynamic laterality, i.e., we investigated how dynamic laterality of a cluster is related to the time-varying functional connectivity of the brain. We then calculated the correlation between time-varying laterality index and time-varying network features (including degree, participation coefficient, and modularity, reflecting brain integration/segregation) and BOLD activity (amplitude of low-frequency fluctuation (ALFF)).
We found that the functional connectivity of the high-order brain networks, especially those in the left hemisphere (e.g., FPN, DMN, and language network, mainly covered by Cluster 4) play an essential role in dynamic changes in lateralization. For the 3 right-lateralized clusters (Clusters 1/2/3), higher level of right lateralization was associated mainly with weaker functional connectivity of these clusters with the higher-order networks on the left hemisphere ( Fig 5B). For Cluster 4 that is intrinsically left-lateralized, its increased left lateralization was accompanied by reduced functional connectivity with multiple lower-order networks at the right hemisphere, including many of the networks covered by Clusters 1/2/3 (the right subfigure of Fig 5B). Stronger functional connectivity of a cluster with its ipsilateral higher-order networks also plays a role in the increased lateralization of the 4 clusters but has less contribution ( Fig 5B; details can be found in S4 Fig).
For functional network features and ALFF, we found that the whole-brain lateralization (the mean absolute value of DLI across all 360 regions) correlated negatively with the averaged degree (r = −0.43 ± 0.1, 91% p < 0.05), participation coefficient (r = −0.2 ± 0.12, 81% p < 0.05), interhemispheric connection (r = −0.52 ± 0.11, 92% p < 0.05), and ALFF (r = −-0.22 ± 0.12, 85% p < 0.05), while showing positive correlations with modularity of the brain network (r = 0.47 ± 0.06, 91% p <0 .05) across time windows; see S5 Fig. This indicates that when the whole brain is highly segregated (modular), it also demonstrates high laterality, while an integrated state of the whole brain (low modularity) are accompanied by low laterality. At the local level, we found that these correlation patterns are most prominent for the high-order brain regions (most significantly in default mode network, the language network, and the frontoparietal network; S5 Fig).
Second, we explored the structural correlates of laterality dynamics using diffusion MRI data. We used diffusion imaging data from a subset of HCP S1200 (HCP Unrelated 100) to constructed 2 kinds of structural connectivity matrices, fractional anisotropy (FA) matrix and fiber number (FN) matrix for each participant. We found significant positive correlations between the laterality correlation matrix and structural connectivity matrix, including FA matrix (Spearman's rho = 0.14 ± 0.02, significant in 99% participants) and FN matrix (Spearman's rho = 0.14 ± 0.02, significant in 99% participants); see Fig [40]. These results suggested that brain regions with similar laterality dynamics tend to have stronger structural connection. Furthermore, LF of a region correlated positively with its FA/FN-degree, i.e., the sum of FA/FN of fibers between this region and all other regions (LF-FN: rho = 0.32 ± 0.09, significant in 99% participants; LF-FA: rho = 0.28 ± 0.08, significant in 99% participants; see Fig 6B).

Heritability of laterality dynamics
Finally, we investigated heritability of lateralization dynamics by twin analyses (using the kinship information). First, we estimated heritability through the ACE model [additive heritability We then calculated the cosine distance between each pair of monozygotic (MZ) twins, dizygotic (DZ) twins, sibling (SI), and unrelated group participants using various dynamic laterality measures across the whole brain. As expected, the laterality correlation matrix showed greater similarity within MZ twins than those within DZ or non-twins (one-way ANOVA, F = 152.7, p < 0.0001; Fig 7D). Similarly, the similarity of MLI and LF (of all regions) in twins, SIs, and unrelated participants showed a decreasing trend, although the difference between MZ and DZ was not significant (for MLI, F = 23.8, p < 0.0001, MZ versus DZ, p > 0.99; for LF, F = 10.8, p < 0.0001, MZ versus DZ, p = 0.99; Fig 7A and 7B). There was no significant difference in LR between the 4 groups (F = 1.48, p = 0.22; Fig 7C). All one-way ANOVA were done using GraphPad Prism 9.3.1 (http://www.graphpad.com). In summary, dynamic laterality measures were generally heritable.

Validation analysis
To rule out the possible influence of potential nuisance variables on our dynamic laterality results, we performed a series of validation analyses: (1) the effect of sample size: all 991 participants were split into two halves, and the main analyses were repeated separately on both halves of the data (S10 Fig); (2) head movements: we repeated the main analyses with multiple head movement parameters (Friston 24 parameters [41]) being regressed out from the BOLD signal (S12 Fig). We also used more stringent exclusion criteria (mean frame distance (FD) < 0.15 mm and FD < 0.1 mm) to repeat our analyses (S12 Fig); (3) GS: we repeated the main analyses, while regressing out the GS of the whole brain (S11 Fig). This resulted in the sum of the GS from the left and right hemispheres being 0, i.e., when the GS L is positive, the GS R is equal to its negative number (rather than a GS of 0 for one hemisphere). Therefore, the correlation coefficients could still be calculated between hemispheric GS and individual ROIs; (4) time window lengths: another 2 window lengths (60 TRs and 90 TRs) were used to investigate the reproducibility of DLI (S13 and S14 Figs); and (5) the potential contaminating effect of an ROI to its ipsilateral GS: when we calculated the correlation coefficient between an ROI and the GS of its ipsilateral hemisphere, the ROI was included in this calculation, which may have "contaminated" its ipsilateral GS. We therefore removed the ROI when extracting the ipsilateral hemisphere GS and repeated our analyses to examine whether our results are affected (S15 Fig). In all the above cases, our main results, including the distribution of dynamic laterality indices across the whole brain, and their correlation with cognitive performance, were largely replicated, which indicated the reliability of the DLIs. See more details in S1 Text.

Discussion
Laterality is traditionally argued to represent a stable, trait-like feature of the human brain. However, considerable amount of work has recently been devoted to characterizing DFC that has been shown to be related to mental flexibility [19] and to predict ongoing mental states or task performance [21], which inspires the question of whether laterality of the brain is also fluctuating over time. This study proposed a framework for dynamic laterality analysis. The time-averaged DLI replicated well previous findings in terms of the static spatial patterns of functional lateralization [12][13][14], while the time-varying laterality measures provided new insights into the dynamic changes of lateralization in resting state. Specifically, different levels of temporal variation of laterality were found across the brain: Regions within the visual and sensorimotor networks generally showed higher variations, while regions from the default mode, frontoparietal, and language network showed relatively weaker variations. This may suggest that laterality in lower-order brain regions is more flexible, consistent with previous research showing that lower-order regions have shorter intrinsic time scale than higher-order regions [42]. That is, the bilateral sensory areas need to process constantly changing sensory inputs and integrate them with information from the contralateral hemisphere in real time to form accurate and coherent percepts.
To systematically characterize dynamic laterality changes, we used 2 measures, i.e., LF and LR, which showed positive and negative correlations with cognitive performance, respectively. These opposing associations suggest that these 2 measures capture different aspects of dynamic laterality. Generally, brain regions in the left and right hemisphere illustrated positive and negative mean laterality, respectively (indicating that they have more ipsilateral connections), with their laterality fluctuating around the mean value (Fig 1). Previous literature has suggested that comprehending literal meanings was associated with left language areas (high left laterality), while dealing with difficult metaphors were shown to involve the right homotopic regions [43]. Greater fluctuations of laterality therefore suggested larger number of states of intra and interhemisphere interactions, which may involve both left hemisphere language areas and the right homotopic regions that may be beneficial for difficult language tasks. For cognitive flexibility (Card Sort task) that generally involves both hemispheres [44] and larger fluctuations of laterality suggested flexible recruitment of both hemispheres that is potentially beneficial.
Although larger fluctuations in laterality correlated with better cognitive performance, extreme changes in laterality, i.e., frequent LR correlated negatively with task performance. LR suggests that laterality of a region fluctuate too much and deviates remarkably from the mean value (the dominant laterality regime; Fig 1D), i.e., a region's functional connectivity switches from being more ipsilateral to highly contralateral. Frequent reversal thus indicates that a brain region constantly leaves its dominant functional regime (i.e., more ipsilateral connectivity). Collectively, these results suggest that moderate changes in laterality may enhance cognitive performance, while extreme laterality changes may hamper cognition. In the Results part, we have shown that low laterality of the whole brain corresponds to high level of interhemispheric communication (or less intrahemispheric information processing), while high laterality of the whole brain corresponds to the opposite case; these results therefore indicate optimal cognitive performance may be related to a dynamic balance between interhemispheric information exchange and intrahemispheric information processing. As we also show that low laterality of the whole brain corresponds to low modularity (integration) while high laterality of the whole brain corresponds to high modularity (segregation/functional specialization), our results also echo the recent findings that cognitive function depends on a dynamic, contextsensitive balance between functional integration and segregation [45][46][47][48][49].
In addition to LF and LR, we furthermore resolved the temporal structures of the wholebrain laterality dynamics. We revealed 3 recurring states (or "meta-states") with distinct laterality profiles, dwelling time, and transition probabilities. These states generally correspond to the multiple lateralization "axes" identified by Karolis and colleagues through dimensionality reduction of 590 meta-analysis maps: Our State 1 is characterized by strong left laterality of the left default mode and language networks, corresponding to the "symbolic communication" axis that involves Broca's and Wernicke's areas. Our State 3 demonstrates strong right laterality of frontoparietal and default mode networks in the right hemisphere consistent with the "active/perception" axis identified in [40].
Laterality of the human brain is hypothesized to be controlled by multiple factors [12], which is supported by the distinct spatial clusters identified by the clustering of laterality time series of all brain regions. Three of the 4 clusters (Clusters 1, 3, and 4) we identified roughly correspond to the top 3 factors identified in [12], which include the visual, attention, and the default mode network, respectively. It should be noted that our clustering analysis was conducted based on laterality dynamics of each individual, rather than on variations across individuals [12].
The 4 spatial clusters also showed distinct patterns of intercorrelations in their time-varying laterality index (Fig 3B). Of particular interest is that the time-varying laterality index of Cluster 4 (the default mode and language network) correlated negatively with those of the other 3 clusters (Cluster 1, visual network; Cluster 2, sensorimotor network; Cluster 3, attention and FPN network) in most participants (Clusters 4 and 1, 99.7% participants with significant negative correlation; Clusters 4 and 2, 99.1%; Clusters 4 and 3, 95.4%). This negative temporal correlation suggests opposite lateralization patterns between Cluster 4 and Clusters 1 to 3 over different time windows. The association analysis between dynamic lateralization and timevarying functional connectivity also suggested the essential role of Cluster 4: The increased level of right lateralization of Clusters 1/2/3 was associated with their decreased functional connectivity with Cluster 4 on the left hemisphere. Similarly, the increased level of left lateralization of Cluster 4 was associated with its decreased functional connectivity with Clusters 1/2/3 at the right hemisphere, suggesting that Cluster 4 might be critical for lateralization state of the whole brain. Furthermore, the more negative the correlation in laterality between Cluster 4 (default mode and language network) and the other 3 clusters (especially visual and sensorimotor network), the better the cognitive performance. This suggested that opposite lateralization pattern between Cluster 4 and Clusters 1 to 3 optimizes parallel processing in the 2 hemispheres [17] and enhances neural capacity. These results could be explained by the causal hypothesis of hemispheric specialization, which states that lateralization of one function forces the other function to the opposing hemisphere, which optimize parallel processing in complex tasks and increases processing efficiency [5,50,51].
Laterality dynamics of the brain may be related to multiple factors. Functionally, greater lateralization of higher-level brain regions was associated with disconnections with the contralateral hemisphere, lower BOLD activity of the whole brain, and greater network segregation, suggesting less exchange of information across hemispheres, which may be related to less energy (glucose) consumption [16,52]. Structurally, LF of a region correlated positively with its degree of structural connections. A brain region with wider structural connections may be modulated by multiple systems that likely demonstrate rich patterns of interhemisphere interaction, thus larger fluctuations in laterality [53,54]. Furthermore, we found that pairs of regions with stronger structural connections are more positively correlated in their time-varying laterality index, indicating the role of structural connection in synchronization of functional lateralization of spatially distributed brain regions. In comparison, brain regions with negative intercorrelation of laterality showed less structural connections, possibly affected by more global factors, such as the regulatory effect of neurotransmitters on large-scale brain networks [39,47,55]. In view of the influence of GS and head movements on the dynamic brain network statistics [56,57], these 2 factors may also affect DLI. To exclude potential influence from these 2 factors, we repeated the above analyses by regressing GS and head movements parameters from the BOLD signal, respectively, and we found that our main results such as the distribution of DLI across the brain and their correlation with cognitive performance remain largely unchanged (see S11 and S12 Figs for details).
At the genetic level, structural brain lateralization is known to be heritable in large population samples [2]. Our research showed that the dynamic characteristics of lateralization were also heritable, especially in those high-order networks such as default mode and frontal-parietal networks, suggesting that laterality dynamics are stable properties of lateralization. From an evolutionary perspective, lateralization arises as a solution to minimize wiring costs while maximizing information processing efficiency in the rapid expansion of the cortex in evolution [58]. Higher-order systems such as default mode and frontal-parietal networks were among the most expanded regions in evolution [59]. Considering the dynamic laterality of these networks correlates significantly with cognitive performance, the heritability of dynamic lateralization in these networks, therefore, may confer evolutionary advantages.

Limitations
In this study, resting-state data were analyzed, so it was not possible to directly relate the results to specific cognitive processes. Future studies should investigate dynamic changes in brain lateralization under task modulation, particularly with multiple task states or various task loads, to better understand the specific cognitive advantages of dynamic brain lateralization. In addition, while the GS-based laterality index has the advantages of being computationally convenient and highly correlated with the traditional ROI-based laterality index, it has the problem of averaging and being too coarse to detect network interactions. Therefore, we further analyzed which specific functional connectivity or network drives the temporal changes of laterality index of the 4 clusters. Finally, age and handedness are important factors potentially affecting lateralization, and handedness is also partly controlled by genetic factors [60]. Future studies are warranted that explore how cognitive deterioration with aging may affect lateralization dynamics.

Conclusions
To conclude, our study demonstrates the dynamic nature of laterality in the human brain at resting state. We characterized comprehensively the spatiotemporal laterality dynamics by identifying 4 spatial clusters and 3 recurring temporal states and demonstrate that the temporal fluctuations of laterality as well as the negative correlation in laterality among different clusters associate with better language task and intellectual performance. We further explored the functional and structural basis underlying such laterality dynamics, and heritability of laterality dynamics. Our study not only contributes to the understanding of the adaptive nature of human brain laterality in healthy population but may also provide a new perspective for future studies of the genetics of brain laterality and potentially abnormal laterality dynamics in various brain diseases.

Ethics statement
This paper utilized data collected for the HCP. The scanning protocol, participant recruitment procedures, and informed written consent forms, including consent to share deidentified data, were approved by the Washington University institutional review board [31]. Data acquisition for the HCP was approved by the Institutional Review Board of The Washington University in St. Louis (IRB # 201204036), and all open access data were deidentified. The data collected by the HCP adhered to the principles of the Declaration of Helsinki. No experimental activity with any involvement of human participants took place at the author's institutions. Our data analysis was performed in accordance with ethical guidelines of the Fudan University ethics committee.

Dataset
We analyzed multimodal brain imaging data and behavioral measures from the HCP 1200 Subjects Release (S1200) [31]. All brain imaging data were acquired using a multiband sequence on a 3-Tesla Siemens Skyra scanner. For each participant, 4 resting-state fMRI scans were acquired: 2 with right-to-left phase encoding and 2 with left-to-right phase encoding direction (1,200 volumes for each scanning session, TR = 0.72 s, voxel size = 2 × 2 × 2 mm). High-resolution T1-weighted MRI (voxel size = 0.7 × 0.7 × 0.7 mm) and diffusion MRI (voxel size = 1.25 mm, 3 shells: b-values = 1,000, 2,000, and 3,000 s/mm 2 , 90 diffusion directions per shell) were also acquired for each participant. A total of 991 participants (28.7 ± 3.7 years old, 528 females) were included in the final cohort utilized in this study according to the following criteria: (1) completed all 4 resting-state fMRI scans; (2) limited in-scanner head motion (mean FD < 0.2 mm); (3) same number of sampling points (i.e., 1,200 volumes per scanning session); (4) without any missing data in regions included in the employed parcellation scheme.

Resting-state fMRI preprocessing
Resting-state fMRI data were preprocessed using the HCP minimal preprocessing pipeline (fMRIVolume) [61] and were denoised using the ICA-FIX method [62]. Then, a spatial smoothing (FHWM = 4 mm) and a low-pass filtering (0.01 Hz to 0.1 Hz) were performed. We did not employ GS regression, as the mean hemispheric time series were used for calculating the lateralization index [14]. The whole brain was parcellated into 360 regions (180 for each hemisphere), using the HCP's multimodal parcellation (HCP MMP1.0) [32]. These regions were grouped into 12 functional networks based on the CAB-NP [34].

Dynamic laterality index
The pipeline for the DLI analysis is illustrated in Fig 1A. The DLI adopted a sliding window approach and a traditional laterality index [14] to capture the dynamics of laterality. DLI at the tth sliding time window is defined by: where ROIi indicates BOLD time series of ROI i, and GS L and GS R indicate the GS, i.e., averaged time series of the voxels within the left and the right hemispheres, respectively. Pearson correlation coefficient r is Fisher-z-transformed. There are 3 main reasons why we adopted the hemispheric GS-based laterality: (1) GS has been shown to be able to reflect and characterize laterality by previous studies [14], and the computational complexity is very low (n � 2, n being the number of regions); (2) the whole-brain laterality pattern obtained using hemispheric GS is highly correlated to that obtained by ROI-based method like AI defined as follows [12,13,33,63]: where P n 1 rðROI i ; ROI j Þ=n and P m 1 rðROI i ; ROI k Þ=m indicate the averaged functional connectivity between BOLD time series of ROI i and all ROIs in the left and right hemisphere, respectively (S1 Fig), and n and m are number of regions in left and right hemisphere.

Spatial clustering of laterality time series across the whole brain
To explore the spatial organization of dynamic lateralization across the whole brain, we performed spatial clustering on the laterality correlation matrix (obtained by calculating Pearson correlation between laterality time series of all pairs of brain regions). Specifically, we applied the Louvain community detection algorithm on the laterality correlation matrix averaged over all 991 participants and 4 runs. We set the γ parameter to 1, ran the algorithm for 100 times, and reported the cluster partitioning with the maximum modularity parameter (Q).

Temporal clustering of whole-brain laterality state
To identify the recurring laterality states of the whole brain from 991 participants and 1,171 time windows, we adopted a 3-stage temporal clustering approach. Firstly, the 1,171 × 4 time windows of the 4 runs of each participant were clustered into 10 states by k-means clustering. Second, all 9,910 states (991 participants × 10 states) were clustered into 1,000 groups by a second level k-means clustering. Finally, we used Arenas-Fernandez-Gomez (AFG) community detection to cluster the 1,000 groups obtained in the second step. AFG allows for multiple resolution screening of the modular structure and a data-driving selection of clustering number. We determined the optimal number of clusters using the modularity coefficients obtained by changing the resolution parameter from 0.1 to 1.5 in steps of 0.1 [64]. The code of K-means and AFG clustering were from MATLAB function kmeans (with the cosine distance metric) and MATLAB Community Detection Toolbox (CDTB v. 0.9, https://www.mathworks.com/matlabcentral/ fileexchange/45867-community-detection-toolbox) [65], respectively. Based on the final clustering, we took the average pattern of each category as centroids and reclassified all the windows of each participant according to the cosine distance between each window and each centroid.

Cognitive performance in behavioral tasks
To understand the cognitive significance of the dynamic brain laterality, we utilized individual performance of a wide variety of cognitive tasks from the NIH Toolbox for Assessment of Neurological and Behavioral function and Penn computerized neurocognitive battery [66]. These involve story comprehension task that is more left lateralized, Flanker Task (inhibitory control and attention) that maybe more right lateralized, and Card Sort (cognitive flexibility) and List Sorting task (working memory) that tend to more bilateral.
The story comprehension task [67] includes a story condition that presents brief auditory stories adapted from Aesop's fables followed by a binomial forced-choice question to check the participants' understanding of the story topic (For example, after a story about an eagle that saves a man who had done him a favor, participants were asked, "Was that about revenge or reciprocity?"), and a math condition to answer addition or subtraction problems. To ensure similar level of difficulty across participants, math trials automatically adapted to the participants' responses. The story and math trials were matched in terms of auditory, duration, attention demand, and phonological input. We used 3 relevant measures: "LanAcc" ("Langua-ge_Task_Story_Acc", the accuracy in answering questions about the story), "LanRT" ("Lan-guage_Task_Story_Median_RT", median response time to answer the questions), and "LanDiff" ("Language_Task_Story_Avg_Difficulty_Level", the average difficulty of all stories a participant can understand).

Correlation between dynamic laterality index and cognitive performance
In calculating the correlation between dynamic laterality measures and cognitive performance (unadjusted scores), we used Permutation Analysis of Linear Models (PALM; http://fsl.fmrib. ox.ac.uk/fsl/fslwiki/PALM) to correct for the bias of significance estimation due to the kinship among participants. PALM used exchangeability blocks that is widely adopted to control the family structure-related bias in HCP data [68,69]. We calculated the Pearson correlation coefficients between dynamic lateralization attributes [MLI, LF, LR of 4 clusters; laterality correlation within and between 4 clusters; fraction and dwelling time of 3 states, state transitions] and cognitive performances, regressing out sex, age, education years, race, body mass index (BMI), handedness, gray matter volume, white matter volume, and head motion (mean FD). We obtained p-values for all correlation coefficients using 5,000 times permutation tests. FDR (q = 0.05) was used for multiple comparison correction (multiply comparison times = 4 MLI + 4 LF + 4 LR + 10 DLI correlation + 3 dwelling time +3 fraction + transition = 29). We show all the results of FDR q < 0.05. Moreover, we also highlight results that survive a more conservative multiple comparisons correction (FDR q < 0.05/13 = 0.0038, where 13 is the number or behavioral measurements used in this study).

Association between dynamic laterality and dynamic functional connectivity
To identify which functional interactions may play important roles in the dynamic change of laterality, we adopted a regression-based approach to investigate which subnetworks (12 pairs of symmetric networks in CAB-NP) may affect the dynamic laterality of the 4 identified clusters. We averaged the functional connectivity within/between the 24 subnetworks; therefore, there are 24 � 23 / 2 = 276 subnetwork-level functional connectivity. Specifically, we built linear regression models of cluster-level DLI time series on subnetwork-level DFC and performed one-sample t test on the regression coefficient (β), to investigate whether the influence of DFC on DLI was significantly greater than 0 (positive coupling) or less than 0 (negative coupling) (Fig 5A).

Association between dynamic laterality and dynamic network properties and amplitude of low-frequency fluctuation
To understand the neural and network basis of dynamic brain lateralization, we correlated the laterality time series of each brain region with time-varying brain network measures (which reflect the organization of the brain like modular, or integration property) and the ALFF averaged over the whole brain for each participant. We constructed functional brain network using Pearson correlation within each sliding time window and computed the following topological measures: 1. Degree centrality (DC) [70,71], measuring the centrality of each node in the network: where G is a graph with the node i and j connected by edge A ij .
2. Participation coefficient (B T ) measuring the contribution of each brain region to wholebrain integration: where κ isT is the strength of the positive connections of region i to regions belong to the module s at time T; κ iT is the sum of strengths of all positive connections of region i at time T. The participation coefficient of a region is close to 1 if its connections are distributed among all of the modules and 0 if all of its links are within its own module. Community division was obtained by community detection in each time window, and then the participation coefficient was calculated [25].
3. Modularity (Q) measuring the tendency of the network to segregate into several independent modules [72,73]: where m is total number of edges in the network; A is adjacent matrix; k i is degree of node i; C i is community assignment of node i; when and only when 4. Interhemispheric connection [74], measuring the averaged interhemispheric connection strength.
Among these indicators, the averaged DC is measured at the global network level, B T and Q are measured at the module level, and interhemispheric connection is measured at the hemispheric level DC, B T and Q were calculated using the Brain Connectivity Toolbox (BCT; http:// www.brain-connectivity-toolbox.net/) [75]. We removed all negative connections in line with the general assumptions of graph theory [71].
ALFF measures spontaneous energetic activity within each time window [76] calculated as follows: where fft is Fourier transform; abs is operator of taking absolute value; BOLD is BOLD signal within a time window in any region; and N is the length of time window. At global level, we calculated the Pearson r between the ALFF/network time series (averaged across the brain) and the absolute laterality time series (averaged across the brain). At local level, we correlated the ALFF/network time series (averaged across the brain) with the absolute laterality time series of each brain region.

Structural network analysis
Because the probabilistic fiber tracking is very time-consuming, diffusion tensor imaging (DTI) data of a subset of HCP S1200 release, the "100 Unrelated Subjects" (n = 100, 54 females, mean age = 29 years) was used for the structural analysis (1 participant was removed due to head motion). Data underwent motion, susceptibility distortion, and eddy current distortion correction [61,77]. We used FMRIB Software Library v6.0 (FSL; https://fsl.fmrib.ox.ac.uk/fsl/ fslwiki/) and MRtrix3, a toolkit for diffusion-weighted MRI analysis (https://mrtrix. readthedocs.io/en/dev/) [78], to construct the SC matrix. The processing steps were as follows: (1) tissue segmentation based on T1-weighted structural image; (2) calculation of 4D images with gray matter, white matter, and cerebrospinal fluid using multi-shell, multi-tissue constrained spherical deconvolution [79]; (3) generation of white matter-constrained tractography using second-order integration over fiber orientation distributions (iFOD2), a probabilistic tracking algorithm (1,000,000 streamlines for each participant, max length = 250 mm, cutoff = 0.06) [80]; (4) use of FSL's FNIRT to reverse-register the MMP parcellation to individual space; (5) calculation of FA maps using dtifit of FSL; and (6) measurement of the average FA on all the streamlines connecting any 2 parcels. This process finally resulted in 360 × 360 FA matrices and FN matrices for 99 participants. We calculated the FA-degree and FN-degree based on this structural connectivity (SC) matrix, i.e., the sum of FA and FN of each row of the SC matrix. The Spearman correlation between laterality correlation matrix and SC matrix and the correlation between MLI/LF/LR and FN-/FA-degree were calculated for each participant.

Heritability analysis
We used twin information in the HCP data to evaluate the heritability of the dynamic laterality measures. In the 991 participants, 103 pairs of MZ twins, 98 pairs of DZ twins, 195 pairs of SI were included, and their zygosity information was used to fit ACE model, which is able to divide total variance in a phenotype (σ 2 ) into 3 additive parts: (A) additive common genetic factor (heritability), (C) common environment, and (E) unshared environment, or other source of measurement error [81]: The ratio of the variance explained by A, C, and E to the total variance was used as estimates of heritability (h 2 ), environmental factors (c 2 ), and error (e 2 ), all of which sum to 1. All 4 models (for MLI, LF, LR, and laterality correlation matrix) were fitted and estimated using APACE (the Advanced Permutation inference for ACE models; www.warwick.ac.uk/tenichols/apace) [82]. Before model fitting, sex, age, education years, race, BMI, handedness, gray matter volume, white matter volume, and mean FD were regressed out from phenotype as covariables. The significance of heritability for each brain region was obtained using permutation test (1,000 permutations per region).

Code availability
The code used in this study to calculate DLI is available on https://github.com/XinRanWu/ Dynamic_Laterality. Higher absolute value represents higher left/right lateralization. We selected data from 3 participants, sub-100206, sub-100307, and sub-100408, to show the relationship between global lateralization (quantified by the averaged absolute value of laterality) and dynamic ALFF/graph theory indicators. The results show that higher global lateralization is associated with lower ALFF and stronger whole-brain dissociation (higher modularization Q, lower participant coefficient, centrality degree, and interhemispheric connectivity). One asterisk ( � ), p < 0.05; 2 asterisks ( �� ), p < 0.01; 3 asterisks ( ��� ), p < 0.001; 4 asterisks ( ���� ), p < 0.0001; ns, nonsignificant. The underlying data for this figure can be found in S5 Data. ALFF, amplitude of low-frequency fluctuation; DLI, dynamic laterality index. The r on vertical axis represents Pearson correlation coefficient between dynamic laterality indictor patterns using window length of 60 TR/90 TR and indictor patterns of 30 TR (used in the text) among all participants. Each dot in the graph represents one of the 991 participants. As can be seen, MLI has the highest reproducibility (close to 1), followed by laterality correlation, LF, and LR. With the increase of the window length, the correlation between the dynamic laterality indicators and the results of the window length of 30 TR decreased gradually. REST1/REST2 indicates the scan time (in first or second day). LR/RL indicates the scanning direction. LR, left to right. RL, right to left. LF, laterality fluctuations; LR, laterality reversal; MLI, mean laterality index. (PNG) S15 Fig. (A-F) The influence of removing the ROI from its ipsilateral GS in calculating the laterality index. That is, we reconstructed GS by removing the ROI A from its ipsilateral GS.  Table. The association between dynamic lateralization characteristic of 4 clusters and cognitive abilities. Each row of the table corresponds to a cognitive test score in the HCP. Pic-Vocab, Picture Vocabulary Test; ReadEng, Reading Test (reading decoding skill); Cardsort, Dimensional Change Card Sort Test (executive function, specifically tapping cognitive flexibility); Flanker, Flanker task (executive function, specifically tapping inhibitory control and attention); ProcSpeed, Pattern Comparison Processing Test (speed of processing); PicSeq, Picture Sequence Memory Test (episodic memory); VSPLOT_TC, Total Number Correct of Penn Line Orientation. PMAT24_A_CR, Number of Correct Responses of Penn Matrix Test (nonverbal reasoning); ListSort, List Sorting Working Memory Test; IWRD_TOT, Total Number of Correct Responses of Penn Word Memory; LanAcc, the accuracy of answering questions about the story of language task; LanRT, median response time to answer the questions of language task; LanDiff, the average difficulty of all stories of language task for each participant, which represents the overall language comprehension ability. r, Pearson correlation coefficient; p, the permutation test significance obtained by PALM. The bold text represents the association remains significant after FDR correction (q = 0.05, comparison number = 29). The asterisk ( � ) represents the association remains significant after FDR correction (q = 0.05/ 13 = 0.0038). See HCP Data Dictionary for more details. FDR, false discovery rate; HCP, Human Connectome Project; PALM, Permutation Analysis of Linear Models. (XLSX) S2 Table. The association between laterality correlation among 4 clusters and cognitive abilities. Each row of the table corresponds to a cognitive test score in the HCP. r, Pearson correlation coefficient; p, the permutation test significance obtained by PALM. The bold text represents the association remains significant after FDR correction (q = 0.05). The asterisk ( � ) represents the association remains significant after FDR correction (q = 0.05/13 = 0.0038). FDR, false discovery rate; HCP, Human Connectome Project; PALM, Permutation Analysis of Linear Models. (XLSX) S3 Table. The association between characteristics of 3 laterality states and cognitive abilities. Each row of the table corresponds to a cognitive test score in the HCP. r, Pearson correlation coefficient; p, the permutation test significance obtained by PALM. The bold text represents the association remains significant after FDR correction (q = 0.05). The asterisk ( � ) represents the association remains significant after FDR correction (q = 0.05/13 = 0.0038). FDR, false discovery rate; HCP, Human Connectome Project; PALM, Permutation Analysis of Linear Models. (XLSX) S1 Data. "S1_Data.xlsx" includes all individual observations used in Fig 2A. (XLSX)