Graph analysis of nonlinear fMRI connectivity dynamics reveals distinct brain network configurations for integrative and segregated information processing

The human brain is organized into functional networks, whose spatial layout can be described with functional magnetic resonance imaging (fMRI). Interactions among these networks are highly dynamic and nonlinear, and evidence suggests that distinct functional network configurations interact on different levels of complexity. To gain new insights into topological properties of constellations interacting on different levels of complexity, we analyze a resting state fMRI dataset from the human connectome project. We first measure the complexity of correlational time series among resting state networks, obtained from sliding window analysis, by calculating their sample entropy. We then use graph analysis to create two functional representations of the network: A ‘high complexity network’ (HCN), whose inter-node interactions display irregular fast changes, and a ‘low complexity network’ (LCN), whose interactions are more self-similar and change more slowly in time. Graph analysis shows that the HCNs structure is significantly more globally efficient, compared to the LCNs, indicative of an architecture that allows for more integrative information processing. The LCNs layout displays significantly higher modularity than the HCNs, indicative of an architecture lending itself to segregated information processing. In the HCN, subcortical thalamic and basal ganglia networks display global hub properties, whereas cortical networks act as connector hubs in the LCN. These results can be replicated in a split sample dataset. Our findings show that investigating nonlinear properties of resting state dynamics offers new insights regarding the relative importance of specific brain regions to the two fundamental requirements for healthy brain functioning, that is, integration and segregation.


Introduction
Fluctuations of the blood oxygen level-dependent (BOLD) signal from the human brain display spatiotemporally organized patterns of functional connectivity (FC) during rest, forming so-called resting state networks (RSNs) [1,2]. Recently, there has been increased interest in the temporal dynamics of FC [3,4], which have been quantified with complexity measures (inter alia), like approximate entropy (AE) [5], Hurst exponent [6], and sample entropy (SampEn) [7]. In studies relevant to the present investigation, the nonlinear dynamics of correlational time series (CTS), obtained with a sliding window (SW) [8], were assessed by computing their SampEn [9][10][11]. These studies were either focused on differences in SampEn between schizophrenic patients and healthy controls [10,11], or the spatial overlap between the SampEn gradient across the brain and established RSNs [9]. However, this methodology has yet to be employed to investigate topological features of different network constellations, dissociated based on the degree of nonlinearity governing the interactions among their constituent parts.
We therefore combined this approach with graphtheoretic measures (GMs) used extensively to describe topological features of brain networks [12,13]. In most graph-based approaches, the edges/connections between two given nodes (corresponding to brain regions) to be included in a graph are selected according to correlation-strength, and the edges falling below some cutoff are being discarded [14]. Here, we created two subgraphs of the fully connected graph, with nodes corresponding to RSNs, and the edges to the SampEn of the CTS between them. Our choice to use SampEn as a complexity index was motivated by the fact that it has been used extensively as an analysis tool for physiological time series (TS) [15][16][17][18]. More importantly for the purposes of our study, its parameter choices (see Sect. 2.5) have been validated for (raw) BOLD time series [7] and CTS derived from SW analysis [9][10][11].
We kept edges with the highest and lowest SampEn, respectively, to obtain a high complexity graph (HCG) and a low complexity graph (LCG). Compatible with the formal definition of SampEn, our subgraphs can be interpreted as representations of distinct functional configurations of the same network forming at different temporal scales: Higher SampEn of the CTS among two given nodes indicates faster dynamics, reflected by the irregularity of the CTS. On the other hand, lower SampEn (more self-similarity of the CTS) corresponds to slower dynamics. This procedure was motivated by the fact that the importance of a node in a brain network (its 'hubness') seems to be frequency-dependent [19,20], and depends on the temporal resolution of the methodology employed [21]. In line with these and further results indicating the temporal specificity of intra-and inter-RSN interactions [22][23][24], we expected significant differences in topology between HCGs and LCGs, at the node as well as the global level.
To test our hypothesis, we used resting state data of healthy subjects from the human connectome project (HCP) [25].
Short summary: • Human connectome project resting state time series • Sliding window analysis of the BOLD time series • Sample entropy calculation of the sliding window time series • Derivation of high complexity and low complexity graphs from the sample entropy matrices • Topological analysis with global and local graph measures 2 Methods

Participants
The whole sample was made up of 812 (N = 812) (410 women; n = 410) healthy subjects from the HCP between the ages of 22 to 37 years. The sample was then split into two sets of 406 participants each to obtain a replication dataset.

Resting state fMRI data
The subjects imaging data consisted of a subset of the high-level resting state connectivity analyses outputs, based on the 2017 HCP data release. BOLD TS contained 2400 time-points, corresponding to two concatenated 15 min sessions, to control for any influences of the phase-encoding direction, which was left to right, and right to left, respectively, for the two sessions. The data were collected with the following parameters: gradient-echo EPI sequence, TR = 720 ms, TE = 33.1 ms, flip angle = 52°, field of view = 208 9 180 (RO x PE), matrix 104 9 90 (RO x PE), slice thickness = 2 mm, 72 slices, 2.0 mm isotropic voxels, multiband factor = 8, echo spacing 0.58 ms, and bandwidth = 2290 Hz/Px. Then data were (pre)processed according to the minimal preprocessing pipeline of the HCP [26,27] as well as subjected to artifact removal [28,29]. Subsequently, a group-based principal component analysis (PCA) was carried out [30]. The group PCA output was used as input for group ICAs to obtain parcellations at different dimensionalities based on data from the whole sample [27]. This was done with the help of FSL's MELODIC toolbox [31,32]. Finally, individual TS for the ICA components were reconstructed with a dual regression technique [33], resulting in one representational BOLD TS per ICA component for every subject.

RSN selection
The RSN selection from the 50 components was carried out in an automatized way by using the component labeler of the GIFT toolbox (https:// trendscenter.org/software/gift/), which uses the templates from the Functional Imaging in Neuropsychiatric Disorders Lab at Stanford University, USA (https://findlab.stanford.edu/functional_ROIs.html) [34]. Each components' spatial map was correlated with these templates and subsequently assigned to the best fitting template (i.e., the corresponding RSN), based on the maximum correlation value (only fits above r = 0.15 were included). This resulted in 23 (N = 23) cortical/subcortical components, corresponding to the following networks (NW): Auditory NW (AUD; n = 1), basal ganglia NW (BG; n = 1), sensorimotor NW (SM; n = 2), visual NW (VIS; n = 5), visuospatial NW (VSP; n = 2), DMN (n = 2), executive control NW (EXC; n = 3), language NW (LNG; n = 2), precuneus NW (PREC; n = 1), and salience NW (SAL; n = 3). VIS was further divided into one pVIS (n = 1) and four higher visual NWs (hVIS; n = 4). The DMN was split up into one dorsal (dDMN; n = 1) and one ventral (vDMN; n = 1) subcomponent. The EXC was partitioned into one rightlateralized (rEXC; n = 1) and two left-lateralized NWs (lEXC; n = 2). Finally, SAL was divided into two anterior NWs (aSAL; n = 2) as well as one posterior NW (pSAL; n = 1). Components with extensive cerebellar activations were excluded from the analysis. The component labeled vDMN by the GIFT toolbox (r = 0.39) also displayed considerable activations in the precuneus and was therefore denoted vDMN/PC, upon visual inspection. The component labeled PRE by the GIFT toolbox (r = 0.33) displayed considerable activations in the cuneus and was therefore denoted PREC/CUN. One component was identified as a thalamic NW (TH; n = 1) upon visual inspection and included in the analysis as such. The abbreviations for the components associated with their corresponding RSNs will be used as proxies for the components from now on. For example, vDMN/PC stands for the IC associated with the vDMN/PC. For a depiction of the RSNs spatial maps, see Fig. 1.

Sliding window
We used a 60-s tapered window (rectangle convolved with a Gaussian: standard deviation [SD] = 1TR), after filtering the BOLD TS from 0.017 to 0.1 Hz to remove frequencies lower than 1/w, where w is the window size. This was done to limit the influence of temporal spurious fluctuations in FC [35]. The window was slid forward in steps of 1TR, resulting in a series of 2317 correlation matrices of size 23 9 23, i.e., in 253 (23 9 22/2) CTS. See Fig. 2 for a graphical description of the workflow described in Sects. 2.4-2.6.

Sample entropy
To compute the SampEn of a given CTS x = [ x 1 ; x 2;ÁÁÁ ; x n ] with length N, one creates an embedding vector with m running data points from x: v i ¼ x i ; x iþ1 ; . . .; x iþmÀ1 ½ , with m being the embedding dimension. Then, for each i (1 B i B Nm) define where r = erx is a tolerance value, e a scaling parameter, and rx the SD of x. H(Á) is the Heaviside function.
and k Á k 1 is the Chebyshev distance, i.e., Averaging over all m embedding vectors, we have that.
Then, the SampEn is defined as.
resulting in a nonnegative number, with higher values indicative of more complexity, and lower values indicative of more regularity of the underlying CTS. We did this for every CTS of every subject, resulting in a 406 9 253 SampEn matrix. We choose the standard parameter values that are used in the literature, i.e., m = 2 and e = 0.20 [9,10,36], to ensure comparability of our results with these past investigations.
To validate our choice of complexity measure (SampEn), we also calculated an alternative, amplitude based, complexity index, namely AE (for a detailed description/derivation of AE, see [36,37]). Subsequently, the 253 SampEn indices were correlated with their corresponding AE counterparts for every subject, with Pearson's r. The correlations were generally high (range: 0.77-0.93, mean correlation: r = 0.87 [first fisher-z and then back transformed]), and statistically significant for every subject (Bonferroni corrected), as well as statistically different from zero (one-sample T test).
Since basic BOLD signal properties can influence complexity measures [38], we calculated the temporal signal-to-noise ratio (tSNR) (mean relative to SD) from the (pre-processed) unfiltered BOLD TS, as well as for the filtered version (bandpass from 0.017 to

Graph measures
To create the two graphs of the network (high complexity vs. low complexity), we binarized the SampEn matrices after selecting the highest (HCG) and lowest (LCG) SampEn connections, respectively. To avoid that the values of our graph measures (GMs) would be biased by a specific threshold, we computed all GMs for a range of density values. The reported results were averaged across thresholds, following the approach for correlation-based graph analyses suggested by [39]. For a graphical description of our workflow, see Fig. 2.

Global efficiency
Global efficiency (GE) is an index measuring the efficiency of information exchange in a network, i.e., its functional integration [40], and is defined as the inverse average shortest path length: where d ij represents the shortest path from node i to node j [41]. We computed the average GE of every subject, for the HCG and LCG, respectively. We then compared their means with a T test for dependent samples.

Betweenness centrality
To assess the importance of a given node, we computed its betweenness centrality (BC) [14], defined by the following expression: where r ij is the total number of shortest paths from node i to node j, r ij v ð Þ is the number of those paths passing through v, and n is the number of nodes in the graph. We did this for all 23 nodes of the HCG and LCG, respectively, for every subject. As null models we computed 1000 random graphs with the same node The sample entropy was calculated for every FC TS. c, d: From the resulting sample entropy matrix two binary adjacency matrices were created, by keeping the edges (connections) with the highest and lowest sample entropy values, respectively. e: Two undirected graphs (high complexity, and low complexity) were created from the adjacency matrices and subjected to further analysis. (Color figure online) degree distribution [42], for the HCG and LCG of every subject. Following a procedure described by [12], we then averaged the BC values of the 23 nodes across the 1000 graphs across all subjects and subtracted this cutoff value from the average BC of every node from the non-random graphs. A node was determined to be a hub if this value was at least one SD above the average difference across all nodes. This way we obtained the hub structure for the HCG and LCG, respectively.

Modularity
As a measure of functional segregation, we computed the Newman-Girwan modularity (MOD) of the graphs, which, given a partition of the network nodes, results in an index Q (if Q = 0 the community structure of the graph is not stronger than one would expect by chance; Q [ 0 if stronger, and Q \ 0 if weaker) [43]. To obtain the optimal partition and deal with the known issue of degeneracy [44], we used the Louvain algorithm [45], as part of a procedure described in [46]: After running the algorithm 100 times, we computed the (23 9 23) association matrix, where entry A ij represents the number of times node i and node j are assigned to the same community across the 100 iterations. From this matrix, we subtracted a null model association matrix obtained from random permutations of the original partitions and kept the values above zero. Finally, we obtained the final partition by running the algorithm on the resulting association matrix. We computed Q for the HCG and LCG for every subject and subsequently compared their means with a T test for dependent samples. To gain more insight into the properties of the hubs, we obtained with the procedure described in 2.6.2 we calculated the within-module z-score (WMZ). The WMZ value of node i is defined as: where m i is the module of node i, k i m i ð Þ the withinmodule degree of i (the number of edges between i and all other nodes in m i ), and k m i ð Þ and r k m i ð Þ are the mean and SD of the within-module m i degree distribution [47]. We also computed the participation coefficient (PC), which, for node i, is equivalent to: where M is the set of modules (given by the partition of the nodes), and k i m ð Þ is the number of edges between i and all nodes in module m [47]. Nodes with high WMZ score and low PC contribute mainly to the modular segregation of a network (provincial hub [PH]), whereas a node with high PC contributes to (inter-modular) integration (connector hub [CH]) [41]. Subsequently, we performed multiple comparisons of these measures between our hubs (HCG and LCG).     As a cutoff, we defined that for a hub to be considered a PH its WMZ score had to be at least a SD above the mean WMZ score across all hubs, which was the case for BG and TH only.

Betweenness centrality/hubs
To compare the hubs in terms of their PC we used a nonparametric Friedman test for repeated measures, which resulted in a significant column effect (Chisquare value = 132.27, p \ 0.001, df = 5). Multiple comparisons (Bonferroni corrected) revealed that BG and TH had significantly higher mean ranks (BG = 4.11; TH = 4.09), compared to the other hubs, which did not differ significantly in their mean rank, with SM2 = 3.13, hVIS1 = 3.04, LNG1 = 3.27, and VSP2 = 3.36. To compare the hubs in terms of their WMZ we used a nonparametric Friedman test for repeated measures, which resulted in a significant column effect (Chi-square value = 131.76, p \ 0.001, df = 5). Multiple comparisons (Bonferroni corrected) revealed that BG and TH had significantly higher mean ranks (BG = 4.07; TH = 4.13), compared to the other hubs, which did not differ significantly in their mean rank, with SM2 = 3.32, hVIS1 = 3.31, LNG1 = 3.07, and VSP2 = 3.11.

Replication data set and additional analyses
In the (split sample) replication data set, all results from the original set could be replicated. We also wanted to assess if our main results were independent of our method to derive RSNs (dual regression, automatized assignment of ICA-derived components to specific functional NWs), and our decision to binarize the SampEn matrices in deriving HCGs and LCGs. We therefore repeated our analysis on a subset Fig. 4 Graphical description of the topology of the high complexity network (left; red) and the low complexity network (right; blue) at a density of 18 percent, i.e., with edges representing the top (high complexity network) and bottom (low complexity network) 18 percent functional connectivity time series in terms of sample entropy. Hubs are surrounded by a black box. Note: The graphs were created by averaging the node degrees across all subjects and then rounding to the nearest integer. Subsequently a representational graph was created with the Havel-Hakimi algorithm [85,86]. This procedure was not part of any calculations to obtain the results of this study. AUD = auditory network; BG = basal ganglia network; TH = thalamic network; SM1, SM2 = sensorimotor networks; VSP1, VSP2 = visuospatial networks; LNG1, LNG2 = language networks; PREC/CUN = precuneus/cuneus network; pVIS = primary visual network; hVIS1, hVIS2, hVIS3, hVIS4 = higher visual networks; dDMN = dorsal default mode network; vDMN = ventral default mode network; rEXC1 = right executive control network; lEXC1, lEXC2 = left executive control networks; aSAL1, aSAL2 = anterior salience networks; pSAL = posterior salience network. (Color figure online) (n = 100) of subjects after choosing a different (a priori) parcellation, based on [49], and subsequently computed weighted versions of our main GMs (GE, MOD, and BC). The resulting effect sizes were equivalent in direction, and either larger (GE; Cohen's d [ 0.80, large size effect), or analogous (MOD; Cohen's d = 0.18, small size effect), compared to the ones from the main analyses. In terms of BC, thalamic and basal ganglia regions retained their hub status in the HCGs, as did higher visual, sensorimotor, and attentional (akin to VSP2) regions in the LCGs. We therefore concluded that our original results were unbiased by subject selection, parcellation, and binarization of the SampEn matrices.

Global efficiency and modularity
The main finding of our investigation is that when RSNs interact with high levels of complexity and irregularity they display a topology that permits more efficient integration of information, compared to when their interplay is governed by lower levels of complexity and changes at a slower pace. On the other hand, the low complexity network (LCN) shows a higher degree of functional segregation, compared to the high complexity network (HCN). Importantly, this pattern of results is replicable across datasets, parcellations, and different ways to derive the HCGs and LCGs from the SampEn matrices. These outcomes are in line with past investigations reporting topological differences of functional RSN configurations at distinct frequencies [19,20,50]. Our outcomes are also compatible with evidence from fMRI that the brain during rest switches between periods where its network structure exhibits different degrees of MOD [51], and with evidence based on magnetoencephalography that the same holds true for GE [12]. Since our results were obtained with a (to our knowledge) new combination of methodological approaches (namely complexity analysis of CTS and graph-theoretic measures), they underscore the validity of the studies mentioned above, in the sense of multiple operationalization. Additionally, they could open the door to gaining new insights into the mechanisms underlying mental illness (see Sect. 4.3 for a discussion).

Hub structure
BG and TH emerged consistently as hubs of the HCN across datasets, and both qualified as a CH as well as a PH, according to our criteria. It has been suggested that TH acts as global/kinless hub capable of interacting with multiple cortical networks [47,52,53], allowing for efficient and integrative information processing [54,55]. Of note, age-related decreases in GE are related to decreased local efficiency of BG and TH structures [56], and, in accordance with its global hub properties in our study, lesions to TH also result in decreased MOD on a network level [39]. The structural basis for the CH properties of TH and BG is well established, of course as constituent elements of the cortico-basal ganglia-thalamo-cortical loop [57,58], but also through extensive reciprocal interconnections of TH (mainly with prefrontal cortex) through which it can exert influence on cortico-cortical activity [59,60].
When compared to the HCN, the overall structure of the LCN seemed to be less determined by properties of a few hubs, since PPC and BC values of the LCN hubs were significantly lower than for BG and TH in the HCN. The average number of clusters was significantly higher in the LCN, as was the average MOD, indicative of a topology that favors functional segregation, juxtaposed with the HCN. The spatial layout of most LCN hubs encompassed regions that have consistently been associated with networks exhibiting functional hub properties [12,[61][62][63]. These were precentral and postcentral gyri (SM2), cuneus (hVIS1), and inferior parietal lobule (VSP2). Interestingly, a similar hub structure emerges when the graph-analyses are based on bandpass-filtered TCs (between 0.01 and 0.075 Hz) [19], although global measures of network topology were not evaluated in this study. This temporal dimension of the (functional) network architecture of the brain is also reflected in our results.

Clinical applications
Our results and methodology can potentially offer new insights concerning neuronal correlates of neuropsychiatric illnesses, such as major depressive disorder (MDD) and schizophrenia (SCZ). With regard to MDD, numerous studies have found diverging FC between patients and healthy controls [64][65][66], but often the temporal specificity of these patterns is not taken into account. However, it has been shown that differences in hub structure between MDD and healthy populations are a function of the dynamics with which these brain regions interact [19]. In light of our results, it would be interesting to investigate changes that are specific to the HCN and LCN, respectively, at the global (GE and MOD), as well as the node level (e.g., the hub structure). Since MDD populations are not homogenous in terms of their symptoms and impairments [67], one can speculate that a HCN/LCN specific analysis could shed light on the physiological underpinnings of these different manifestations of the disease; however, testing these hypotheses has to be the focus of future investigations.
The same holds true for the study of SCZ patients, where affective and psychotic symptoms can affect subpopulations to a different degree [67,68]. Additionally, given that BG dysfunction in SCZ has been established by a plethora of studies and meta-analyses [69][70][71], the global hub properties of BG found in the present investigation would predict decreased GE, as well as MOD, in SCZ patients, compared to controls. There is indeed evidence that this seems to be the case [72][73][74][75], and interestingly, TH dysfunction has also been found in SCZ [76,77], which is also in line with our results. Of note, after calculating the SampEn of CTS obtained from a SW (analogous to our study), [11] reported increased complexity within the visual system of SCZ patients, as well as within and between subcortical and cortical structures [10]. In light of our results, and the findings mentioned above, one can speculate that this might be reflective of a compensatory mechanism [78], to account for decreased effectiveness of BG and TH in ensuring proper information processing. Apart from MDD and SCZ, we think that our approach in general can offer new insights to the workings of the healthy, aging, and diseased brain, and opens the door for future investigations of nonlinear brain network interactions across different populations.

Limitations
There are some caveats inherent to our approach. First, the use of a SW to assess FC dynamics requires determining the values of two parameters (window size and SD) a priori; however, a substantial body of literature exists that suggests that our choices are close to optimal for BOLD time series [35,[79][80][81]. It should be noted that a number of other ways for estimating dynamic FC exist, such as spatial distance, innovationdriven co-activation patterns, and jackknife correlation [82][83][84], that offer a frame by frame temporal resolution, whereas (by its nature), SW analysis rather captures slower FC dynamics. Applying our approach to these alternative methods is an exciting prospect for future investigations. Secondly, in terms of our graph analyses, we concede that the choice of the cutoff values (determining the density of the graphs) is somewhat arbitrary, but by averaging the GMs across a range of cutoff values we are confident that our results are not biased by a specific threshold. Finally, another issue pertains to the influence of noise. By definition, SampEn is sensitive to noise in two ways: Unstructured noise would result in higher SampEn values, whereas structured noise would deflate Sam-pEn values. Nevertheless, the fact that SampEn values were generally not significantly correlated with the tSNR of their BOLD TS indicates that they were not simply a function of the noise content of the underlying TS.

Summary
We have shown that subcortical networks are the central nodes in a large-scale brain network whose inter-node interactions are marked by fast dynamics and high levels of complexity. The topology of this HCN favors the integration of information, juxtaposed to a LCN governed by slower dynamics. In contrast, the LCNs structure displays higher levels of MOD than the HCN. These results are important, since they show that investigating nonlinear properties of resting state dynamics can reveal the relative importance of specific brain regions to different fundamental requirements for healthy brain functioning.
Funding Open Access funding enabled and organized by Projekt DEAL. No funding was received to assist with the preparation of this manuscript.
Data availability The HCP data are freely available at http:// www.humanconnectomeproject.org/data/.

Conflict of interest
The authors have no conflicts of interest to declare that are relevant to the content of this article.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.