Temporal flow of hubs and connectivity in the human brain

Hubs in brain network connectivity have previously been observed using neuroimaging techniques and are generally believed to be of pivotal importance to establish and maintain a functional platform on which cognitively meaningful and energy-efficient neuronal communication can occur. However, little is known if hubs are static (i.e. a brain region is always a hub) or if these properties change over time (i.e. brain regions fluctuate in their 'hubness'). To address this question, we introduce two new methodological concepts, the flow of brain connectivity and node penalized shortest paths which are then applied to time-varying functional connectivity fMRI BOLD data. We show that the constellations of active hubs change over time in a non-trivial way and that activity of hubs is dependent on the temporal scale of investigation. Slower fluctuations in the number of active hubs that exceeded the degree expected by chance alone were detected primarily in subcortical structures. Moreover, we observed faster fluctuations in hub activity residing predominately in the default mode network that suggests dynamic events in brain connectivity. Our results suggest that the temporal behavior of connectivity hubs is a multilayered and complex issue where method-specific properties of temporal sensitivity to time-varying connectivity must be taken into account. We discuss our results in relation to the on-going discussion of the existence of discrete and stable states in the resting-brain and the role of network hubs in providing a scaffold for neuronal communication across time.


Introduction
The brain's structure and function measured with fMRI and M/EEG have increasingly been analysed with network models. One of the important targets in network modeling of brain activity has been to localize and characterize brain regions (nodes) that play a role as hubs ( Sporns et al., 2007 ). Hubs, in the brain, are believed to play an important role for integrating and distributing information by minimizing the distance between nodes and subsequently lower the metabolic cost of information transfer ( van den Heuvel and Sporns 2012 ). To date, several methods to locate hubs in the brain have been suggested. These include (i) rich clubs, an interconnected ensemble of high degree nodes ( van den Heuvel and Sporns 2011 ); (ii) participation coefficient, where nodes that have many between community connections are considered central ( Guimera et al. 2005 ;Power et al., 2013 ;Gordon et al., 2018 ) (iii) within-module degree z-score, where nodes that have many within community connections are considered central ( Guimera et al. 2005 ;Power et al., 2013 ;Gordon et al., 2018 ); (iv) betweenness centrality (BC), which are nodes that are included in many shortest paths between nodes, where a shortest path is the distance for a node to travel to another node (e.g. Buckner et al., 2009 ;Fransson et al., 2011 ). The latter has fallen out of fashion in functional connectivity studies, partly due to skepticism regarding how correlation values translate to path length between nodes ( Power et al., 2013 ). These different methods to identify hubs will identify different nodes as their candidate hubs, as shown in Fig. 1 A.
There has been a focus on time-varying connectivity (TVC) studies of the brain in recent years ( Laumann et al., 2017 ;Lurie et al. 2020 ). In this context, connectivity is quantified across time. However, when considering connectivity across time, there are (at least) two distinct ways that edges can represent connectivity. First, each estimate can represent the magnitude of connectivity (e.g. each edge could be a Pearson r-value representing the amount that two brain areas co-vary over a short time-window). An alternative view is to consider the flow of connectivity, that is the increase and decrease in connectivity through time relative to its baseline or average connectivity values (see Fig. 1 B). While magnitude-based connectivity is currently standard, there are several reasons to prefer flow TVC in some circumstances. Because there is a relationship between the time-averaged magnitude/variance of time-varying connectivity with functional connectivity Fransson 2015a,2015b ), magnitude-based estimates of TVC include both intrinsic (static) and time-varying properties which cannot easily be disentangled. By considering the flow of connectivity instead, the intrinsic connectivity is removed, which helps identify changes in connec- what different candidate hub identification measures (hubs illustrated with darker colors) will quantify as a hub on a static network (assuming all connected nodes are binary edges and colors denote communities). ( B ) Depiction of different ways to quantify time-varying functional connectivity (TVC). The first "Magnitude TVC " quantifies some statistical value. The second "Flow TVC " quantifies the variation of this statistical value (e.g. relative to the time-averaged connectivity). The choice of TVC edge-type changes subsequent network property measures, such as the shortest path (right). ( C ) Using flow TVC, where thick lines are increases in connectivity, the hubs for three different time-points using betweenness centrality are illustrated. tivity that would have been drowned out by the intrinsic connectivity ( Fig. 1 B). To give an intuitive example of this: imagine two roads connecting two cities (one motorway and one rural road). A hundred percent increase of traffic along the rural road will not entail quantitatively more cars on the rural road compared to the motorway (which always will be at capacity). However, considering the flow of traffic, it would signify a substantial change compared to the norm. Quantifying connectivity as the flow of relative information will allow for changes of both the edges and the following properties that can be investigated ( Fig. 1 B). While not necessarily beneficial in all circumstances, time-varying connectivity estimates representing connectivity flow provide an opportunity to identify new connectivity differences that would otherwise go amiss.
Much of the focus on nodal properties through time has focused on magnitude-based time-varying connectivity, which can lead to a direct mapping of the static functional properties (e.g., Fig. 1 A) such as the participation coefficient ( Shine et al., 2016 ;Betzel et al., 2016 ). Moreover, network measures such as the participation coefficient and module degree z-score have been commonly used in the literature to detect the time-varying properties of network communities. However, caution is warranted when static network measures are used to extract timevarying network properties since it is hard to know the ground truth ( Fortunato and Hric, 2016 ). A pertinent example where a static network measure applied in a time-varying context might cause misinterpretations is the usage of the participation coefficient to extract community partitions across time ( Thompson et al., 2020 ).
It can be hard to disentangle the time-varying properties of nodal measures with magnitude-based TVC networks compared with their intrinsic properties. However, when switching to flow TVC, the interpretations of these properties change. This is exemplified in Fig. 1 C, where increases in connectivity are marked. In the three examples, only betweenness centrality will isolate consistently interpretable properties by identifying nodes that, due to their increase in connectivity, will become highways for the shortest paths calculations. PC and z are problematic as they must relate their changes to a community and will be unable to isolate the importance of the link to integration (red candidate hub in Fig. 1 C, left). Changes in connectivity do not need to be modular, and it is possible that a "chain " of nodes increases in connectivity (e.g., Fig. 1 C, right), which RC would not detect using flow TVC. Nor could it detect multiple keys, but unconnected, nodes in the flow (e.g., Fig. 1 C, middle). In sum, BC is a good candidate metric to identify hubs in flow TVC.
In this article, we apply betweenness centrality on flow TVC calculated from resting-state fMRI data. We start by calculating TVC using two methods that differ in their sensitivity to faster (jackknife correlation) and slower (sliding window) changes in signal co-variance. Next, the flow connectivity estimates can be transformed into distances for calculating the shortest paths. Importantly, we introduce and evaluate the effect of a node visitation penalty when computing the shortest paths. This step alleviates previous criticism against applying betweenness centrality on functional connectivity data. We show that the number of active hubs fluctuates across time, but different constellations of active hubs are present at different time scales. In the discussion section, we review the significance of the current findings in the context of the presence of discrete brain states as well as discussing the merits of the proposed method in relation to studies in comparative anatomy.

Data used
The 100 independent subject dataset from the Human Connectome 500 subject release was used ( van Essen et al., 2012 ;Smith et al., 2013 ) to study the flow of brain connectivity in the form of timeresolved measurements of nodal betweenness centrality during rest. Subject recruitment procedures and informed consent forms were approved by the Washington University institutional review board. The dataset is publicly available on the ConnectomeDB database ( https://db.humanconnectome.org . Whole-brain echo-planar imaging was performed at 3T using a 32 channel head-coil with TR = 720 ms, TE = 33.1 ms, flip angle = 52°, in-plane FOV 208 ×180 mm, 72 slices, 2.0 mm isotropic voxels, multi-band acceleration factor = 8. Further information regarding the MR acquisition parameters can be found in ( U ğurbil et al., 2013 ). As a training dataset for the methodological development and parameter fitting, we used the 'LR' resting-state data set. We used an independent dataset (the 'RL' dataset) as a discovery dataset.
The fMRI data used had undergone image preprocessing and removal of artifacts from non-neuronal origins employing of the FIX ICA (FM-RIB's Independent Component Analysis-based X-noisifier) data artifact rejection process which removed ICA components from the data that were considered to constitute signal contributions from white matter, cerebrospinal fluid, head movement, cardiac and respiratory sources ( Glasser et al., 2013 ;Griffanti et al., 2014 ;Salimi-Khorshidi et al., 2014 ).

fMRI data preprocessing
We used the Power ( Power et al., 2011 ) parcellation of the cortex and subcortical regions to extract BOLD signal intensity time-courses from 264 spherical Region-of-Interests (ROI, radius = 5 mm), where each ROI is assigned to a subnetwork. The ten subnetworks studied were DMNdefault mode (58 nodes), SM -sensorimotor (35 nodes), VIS -visual (31 nodes), FP -frontoparietal attention (24 nodes), SA -saliency (18 nodes), CO -cingulo-opercular (14 nodes), AU -auditory (13 nodes), Sub -subcortical (13 nodes), DA -dorsal attention (11 nodes), VA -ventral attention subnetwork (9 nodes). Thirty-eight nodes were unassigned (U). To minimize the effect of movement ( Power et al., 2012 ;Dijk et al., 2012 ), image "scrubbing " was performed. The framewise displacement (FD) for each image volume was computed (rejection at FD > 0.5). Rejected volumes were deleted and missing data-point were estimated using a cubic spline interpolation ( Thompson and Fransson, 2015b ). The method of estimating deleted data with a cubic spline interpolation after image scrubbing has been done in previous time-varying connectivity studies, see for example Allen et al., 2014 . On average, 1.45% of the image volumes per subject were interpolated (across both 'LR' and 'RL' datasets). Subsequently, fMRI signal intensity time series were band-pass filtered (0.01-0.2 Hz). fMRI resting-state data from one subject contained missing data and was not used in the analysis. Additionally, we performed a linear regression on each BOLD signal time series in which framewise displacement (FD) values and the mean ROI signal (taken across all ROIs at each point in time) were regressed out (see Fransson et al., 2018 for motivation).

Calculation of time-varying connectivity time series
There exist a variety of methods to compute time-varying estimates of connectivity for BOLD fMRI time series. Importantly, several studies have shown that methods differ in their sensitivity to detect slowly occurring versus relatively faster, event-like changes in functional connectivity Ochab et al., 2019 ;Xie et al., 2019 ). Since we in this work were interested in capturing both slow and fast evolving changes, we opted to use two different TVC methods. To this end, we used the Jackknife correlation method which has previously been shown to have superior performance to estimate single time-point resolved changes in functional connectivity  and has been used in multiple TVC studies since (e.g. Fransson et al. 2018 ;Shine et al., 2019 ;Ghahari et al., 2020 ).
The jackknife correlation between two time-series x and y at t can be expressed as which is the Pearson correlation between x and y overall time-points except for x and y at time-point t . The ̄ and ̄ are the expected values, excluding the data at time point t : The minus sign in Eq. (1) corrects for the sign inversion caused by the jackknife correlation. Our motivation for the jackknife correlation method was based on the study that compared multiple methods to derive point-based estimates of time-varying connectivity . Due to variance compression of the jackknife correlation, individual connectivity estimates yield little meaning by themselves, but, importantly, they become meaningful in relation to the rest of the time series of connectivity values. The jackknife correlation method derives changes in connectivity relative to the other time-points. Hence, it is an apt method to derive the flow of connectivity. To capture slower changes in time-varying connectivity, we used the used sliding window (SW) method ( Kiviniemi et al., 2011 ;Leonardi and Van de Ville 2015 ) using a window length of 90 TR time frames (65 s).
In order to create flow connectivity representations, time series for both TVC methods were Fischer transformed and then converted into Zvalues by subtracting the mean and dividing by the standard deviation. This entails that each time series has a mean of 0 (average connectivity value) and a standard deviation of 1. All time-varying functional connectivity methods were derived using teneto ( Thompson et al., 2017 ; https://github.com/wiheto/teneto ).

Mapping correlation-strength into distance weights for time-varying fMRI data
Since previous criticisms of betweenness centrality have raised concerns regarding the biological plausibility of using a connectivity matrix of correlation coefficients to calculate shortest paths in a network ( Power et al., 2013 ), we applied a non-linear transform to these values to create data-driven distributions optimized to isolate hubs. We tested multiple different functions and scaling parameters on the training dataset. The scaling aimed to maximize the variance of the betweenness centrality estimates across nodes. The rationale here was that, if there are nodes that have efficient functional pathways, this optimization function will identify them. Using the results from the training dataset, we opted to use an exponential correlation-to-distance mapping function: ( ) = − + with the distance penalty constant k set to 1.25 and 1.5 for the jackknife correlation and sliding window TVC methods, respectively. A full account of how correlation strength was mapped into distance, including the derivation of the node penalty, is given in Supplementary methods and Supplementary Figures S1-S5.

Surrogate data
In order to demonstrate time-varying properties of interest, we computed estimates of time-varying BC based on surrogate data to infer whether the temporal profiles of changes in nodal BC display timevarying properties or not. To this end, we used the Fourier phase randomization method (same phase randomization applied to all nodes) to produce surrogate BOLD signal intensity time series because it preserves the stationary statistical properties of the signal such as the degree of auto-correlation, co-variance structure and the shape of the power spectrum properties whereas time-varying changes in the co-variance structure is destroyed ( Lancaster et al., 2018 ). Surrogate data from phase randomization of BOLD data from the 100 subjects were used as input to both the sliding-window and jackknife method and surrogate TVC fMRI data was obtained. Next, the same correlation strength to distance mapping function was used to the surrogate TVC data from which the surrogate-based BC values were computed.

Identifying hubs
We use the term "candidate time-varying hubs " (sometimes just "candidate hubs " for brevity) and "active hubs " and these have two distinct meanings. Candidate time-varying hubs are those nodes which appear to have a hub-like role in their time-averaged properties compared to the other nodes. To identify the candidate time-varying hubs we apply a threshold of one standard deviation above the average BC value. This selection procedure helps us select which nodes to focus on. Since these threshold boundaries are somewhat arbitrary, we also provide supplementary figures for increasing/decreasing the threshold. Further, when contrasting the candidate hubs to the surrogate data, note that it is not the BC-magnitude that is evaluated in the contrast but rather whether a node maintains its topographic role as candidate hub role in the surrogate data.
The candidate hubs consider the entire time series. Active hubs however are when the candidate hubs have high BC value during the timeseries. Again we apply a threshold to determine these time-points. An active hub is when a candidate hub's time-varying BC value is larger than its standard deviation of its own BC value.

Identifying candidate time-varying hubs
As a starting point, we first calculated the temporal mean of BC (averaged across 1200 (jackknife correlation) and 1110 (sliding window) time-points, respectively) and subsequently averaged across subjects.
The results shown in Fig. 2 suggests that the set of candidate hubs that are associated with large BC values (see Supplementary Figures S7 and  S8 for results employing lower and higher thresholds to mark the candidate hubs) across time is clearly sensitive to the TVC method used. So, by judging from the temporal mean of BC, there is a marked difference in anatomical localization of candidate hubs, a difference that is dependent on the method's sensitivity to the speed of change in brain connectivity. The jackknife TVC method, which is sensitive to quicker fluctuations in connectivity, singles out candidate hubs that mainly reside in the medial prefrontal and lateral parietal Default Mode, left Fronto-Parietal and Salience networks ( Fig. 2 A-C). On the other hand, when using the sliding window method that are sensitive to slower changes in brain connectivity, we primarily found high BC values in subcortical regions ( Fig. 2 E-G). A total of 11 (3 DMN, 2 FP, 2 SA, 1 VIS) and 23 (3 DMN, 3 SA, 2 CO, 2 AU, 13 Sub (all nodes)) nodes were flagged as "fast " and "slow " flow hubs, respectively ( Fig. 2 A and E), depending on which method was used.
In the surrogate jackknife data (marked in red in Fig. 2 ), we observe no candidate hubs. Thus, we conclude that it is very likely that the null hypothesis does not hold for the actual jackknife data entailing that there are indeed candidate time-varying hubs. This finding implies that on a shorter time-scale, our results suggest that there are "event "like signal changes in TVC fMRI that can be attributed to time-varying changes in co-variance. In contrast, the surrogate sliding window data yielded a result that is more similar to the empirical data but with marked lower overall values of the temporal mean of BC. Eight of 23 candidate hubs were retained in the surrogate data to be a candidate hub, primarily subcortical nodes, and these nodes should not be considered significant.
Additionally, we were here also interested if the average (across all subjects) of the number of active hubs varied during the resting-state fMRI experiment. The active hubs are shown in Fig. 2 D, (jackknife correlation results, fast hubs) and Fig. 2 H (sliding window results, slow hubs). Approximately 25% of all of the fast hubs (3/11) were, on average, active at any given time during the resting-state fMRI experiment ( Fig. 2 D). For slow hubs, slightly less than 50%, on average, (10/23) were active at any given point in time ( Fig. 2 H).
At first glance, it might seem that a fewer number of hubs are located in the left compared to the right hemisphere. This observation is particularly so for the hubs based on the jackknife correlation method shown in Fig. 2 A and C. However, when employing a more lenient threshold, as shown in Supplementary Figure S7, a more symmetrical pattern of hubs across hemispheres is obtained. This result suggests that the apparent asymmetrical pattern of hubs shown in Fig. 2 is more related to variability in sensitivity than caused by differences in communication routes between hemispheres. MNI coordinates and anatomical labels for fast and slow candidate hubs are given in Supplementary Tables 1  and 2.

Identifying active hubs
Although Fig. 2 D and H might suggest that the number of active hubs for both fast and slow hubs is constant across time, we next focused on the temporal and spatial behavior of the fast and slow candidate hubs a function of time. The temporal pattern of the number of active hubs using the jackknife and the sliding window approach is exemplified for a single subject in Fig. 3 A and D, respectively. Interestingly, the results for a single subject suggest that the number of active hubs oscillates in time, particularly for the case of slow hubs ( Fig. 3 D), but also to a lesser extent for fast hubs ( Fig. 3 A). Similar oscillations in the number of active hubs were found in all subjects, with typical fluctuations ranging between 0 and 8 (fast hubs, jackknife method) and between 3 and 20 (slow hubs, sliding window method). When computing the frequency of the time series of active hubs, we note that both TVC methods yield lowfrequency temporal oscillations for the number of active hubs during the . Additionally, the temporal mean of BC was calculated for surrogate data (phase randomized) and the results are shown in red in panels A and D. Nodes for which the temporal mean of BC were larger than the mean + 1 std taken across all subjects (candidate hubs) are marked with asterisks (blue = empirical data, red = surrogate data) in panels A and E (11 and 23 nodes). The anatomical positions of the candidate hubs are displayed in panels C and G . Panels D and H show the average number of active candidate hubs as a function of time (all subjects, shaded area shows 95% confidence intervals). The results for active candidate hubs across time for phase randomized data in panel H (marked in red) were based on the nodes marked with asterisks in panel E. Note that the inclusion of phase randomized data with the sliding window method was to show the temporal properties of the possible active hubs in the surrogate data. Note that the temporal mean BC for surrogate data using the jackknife method did not produce any significant hubs (panels A and D). See Supplementary Tables 1 and 2 for  information on anatomical localization. resting-state fMRI experiment (single subject Fig. 3 BD, the average of all subjects in Fig. 3 CF).
Of note, the results in Fig. 3 E and F suggest that slow oscillations in the number of active hubs are, to some extent, also present for surrogate data when using the sliding window approach. Importantly, as the difference spectrum shown in Fig. 3 H suggests, the slow oscillations in the number of active hubs exceeds that what can be expected from chance alone (see also Supplementary Figures S9 and S10 for corresponding results using a lower (mean + 0.5 std) and a higher (mean + 1.5 std) threshold, respectively).
The results in Fig. 3 support the notion that the number of active hubs oscillates in the low-frequency range during rest. Further, our findings show that with flow TVC, there are regular intervals when the brain switches between a larger group of hubs to a smaller group of hubs. One can hypothesize that switching between these different number of hubs provides a beating connectivity pulse of integration and segregation between brain regions.

Sub-network-dependent properties of active hubs
For the results shown in Fig. 3 , we addressed the question of recurrent patterns of active hubs exclusively from a quantitative perspective. These results do not tell us if the spatial configuration of active hubs during troughs and peaks as exemplified in Fig. 3 A and D differed in any systematic way. First, we ask the question if there exists a difference in subnetwork membership for active hubs during epochs of low versus high hub activity. If so, it would imply that epochs of low and high hub activity are not only quantitatively different but also different in terms of their spatial affinity to different subnetworks in the brain. To address this question, we picked out all points in time for which the number of  Fig. 2 A and E. Results marked in red signify surrogate data (phase randomized). Panel A shows the number of active hubs as a function of time when using the jackknife method to compute TVC in a single subject. Similarly, panel D shows the results when using the sliding window approach to compute TVC for the same subject. Panel B shows the corresponding amplitude spectra for fluctuations in the number of active hubs (jackknife method). Similarly, panel E shows the amplitude spectra for the same individual (panel D) when using the sliding window approach. Panel C shows the amplitude spectra for fluctuations in the number of active hubs when averaged across all subjects and using the jackknife method. Panel E shows the corresponding results for the sliding window method. For the sliding window method, the difference in amplitude spectra between empirical and surrogate data is shown in panel H (shaded area shows 95% confidence intervals). Note that the results for surrogate data using the jackknife approach in panels A-C are lacking, since no nodes were significantly active across time (see also Fig. 2 A). This is also the reason why the difference in amplitude spectra between surrogate data and empirical data for the jackknife method is not shown.

Fig. 4. Presence of slow temporal fluctuations in the total number of active hubs.
Time-points where the number of active hubs falls within the lower 10percentile or upper 10-percentile are exemplified for one subject in panels A and B (left, blue dots and red dots). The right column shows in the form of pie charts the distribution of subnetwork membership for active hubs during time-points that fall within the lower 10 and upper 10 percentiles averaged across all subjects. The number within parenthesis for each subnetwork is the mean percentage based on all time-points and all subjects. Note that all percentages given in the pie charts are adjusted for the number of nodes in each network. The results shown in panel A pertains to fast hubs based on jackknife TVC, whereas panel B display results for slow hubs that were derived using sliding window TVC. See also Supplementary Figure S15 for a presentation of phase randomized data using the sliding window TVC method.
active hubs was within the lower 10-percentile and upper 10-percentile in all subjects. For both percentiles, we calculated the relative number of active hubs in each subnetwork. The results are shown in Fig. 4 A (jackknife correlation method, fast hubs) and Fig. 4 B (sliding window method, slow hubs). Additionally, the results for the sliding window approach applied to surrogate data are shown in Supplementary Figure  S15.
The results in Fig. 4 B suggest that the relative distribution of slow hubs across subnetworks are closely matched during peaks and troughs. Moreover, the spatial location of active hubs during both troughs and peaks do not disperse from the localization pattern obtained for the overall mean number of active hubs (these numbers are given within parenthesis in the pie-charts shown in Fig. 4 ). A notable exception to this observation is the Cingulo-opercular subnetwork that during both peaks and troughs is slightly overrepresented compared to the time-averaged share of active hubs.
In the case of fast hubs ( Fig. 4 A), the lower 10-percentile included only time-points where no active hubs were detected, whereas the spatial distribution of active hubs for the upper 10-percentile were again closely matched with the time-averaged distribution. So, although we have established that the number of active number of hubs oscillates across time ( Figs. 3 and 4 ), neither fast nor slow hubs showed any signs of deviations from the mean spatial patterns of active hubs.

No co-occurence of active hubs argues against discrete connectivity states
With the sets of active hubs for each TVC method in mind, we investigated whether any spatial pattern of active hubs reoccurs over time. If true, this would suggest that the brain during rest revisits specific spatial constellations of hubs. Importantly, such a finding would support the brain cycles through multiple connectivity patterns of discrete states during rest. This notion of multiple discrete states has been suggested repeatedly in the literature ( Allen et al., 2014 ;Gonzalez-Castillo et al., 2015 ;Hansen et al., 2015 ;Karahanoglu and Van De Ville, 2015 ;Vidaurre et al., 2017 ) and not merely forced into discrete clusters by the algorithms employed.
However, while our results for the temporal behavior of hubs have not provided any unequivocal support for existence for recurring, durable, and spatial well-defined global patterns of shortest paths in brain connectivity, the possibility that such recurring patterns of hub activity exist cannot be ruled out. It may not be distinguishable from each other at the level of subnetworks shown in Fig. 4 . It is conceivable that spatial patterns of time-lasting and recurring hub activity can be characterized only at the level of individual nodes. To investigate this possibility, we computed spatial overlap matrices based on BC values for all candidate hubs at all points in time in individual subjects. The degree of spatial overlap of active hubs during different time-points was computed using the Jaccard index and displayed in the form of a time x time symmetrical matrix. A matrix entry of 0 means that no active hub was present at both points in time, whereas a value of 1 means that all active hubs (100 percent) present at t1 were also present at t2.
The degree of spatial overlap for active hubs at different time points are shown in Supplementary Figure S11-13. The results using the jackknife correlation method shows that for neighboring points in time, the degree of spatial overlap for active hubs is slightly above 40 percent and then rapidly decreases to chance levels for larger differences  Figure S11B). Conversely, the results from using the sliding window method show a substantial spatial overlap for active hubs in time, starting at more than 80 percent for adjacent timepoints and then displaying a slow return to chance levels after 80 -100 TRs (Supplementary Figure S11D). This behavior of hub configurations across time is not surprising given that TVC in the sliding window case is calculated for 90 consecutive time-points. Thus, the large degree of spatial overlap at different time-points does not reflect long-lasting constellations of active hubs that slowly vanish and are reconfigured into new configurations. Rather, it is a consequence of the usage of the sliding window method to calculate TVC. Importantly, neither the results based on the sliding window nor the jackknife correlation method shows any evidence of hubs' recurring spatial configurations.
Thus, our results have not provided support for any kind of presence of global configurations of hubs that reoccur through time. Nevertheless, an interesting question can be asked regarding pairwise co-occurrences of hubs in time. In other words, we are interested in whether there is a pattern of how often a pair of candidate hubs are active hubs at the same point in time. If some pairs of hubs are relatively more commonly engaged at the same time than other pairs of hubs, this would indicate a preference for specific pairwise engagement of hubs. The results for co-occurrence, when averaged across all time-points and subjects, is shown for jackknife TVC based results in Fig. 5 AB and for sliding window TVC methods in Fig. 5 CD, respectively. Interestingly, for the case of slow hubs, a hub located in the right anterior insula (located in the SA sub-network) shows the highest degree of co-occurrence with several other hubs in the DMN as for the thalamus and putamen within the Sub-cortical subnetwork.
In sum, our results suggest, while there are no reoccurring patterns of time-varying whole-brain hubs, there are reoccurring pairs of hubs. This finding would indicate regular "activation chains", as described in Fig. 1 C, occurring at regular intervals (supported by Fig. 1 C). However, these occurrences are local (i.e., specific hubs increase their activity) instead of global nature (all hubs increasing their activity).

Discussion
Several conclusions can be drawn from this work. First, our results suggest that the flow of connectivity in hubs does not fit previously reported observations of TVC as a succession of smooth (e.g., Allen et al., 2014 ) or instantaneous (e.g., Yaesoubi et al., 2015 ) transitions between discrete states. Rather, our results suggest that there are no recurring spatial constellations of nodes that act as hubs for changes in brain connectivity. In the current context, this implies that we did not find evidence for any specific global set of shortest paths that re-occur together repetitively during a resting-state fMRI scan. It could be argued that state-like features in flow TVC may not necessarily be distinguishable in terms of shortest paths. However, as outlined in the introduction and Fig. 1 , we think that this scenario is unlikely since betweenness centrality has a unique capability to detect consistently changed properties of flow TVC that is not dependent on differences in static functional connectivity, but rather targets on-going changes in integration and segregation between networks.
It is important to point out that our results pertain to time-varying activity in hub activity only. Therefore, we do not claim that our findings are incompatible with previous findings that have suggested that the resting brain traverses smoothly ( Allen et al., 2014 ) or instantaneously ( Yaesoubi et al., 2015 ) between discrete brain states that can be defined spatially as well as temporally. However, from the perspective of timevarying hub activity, we find that our observations fit well with the view offered by recent studies of metastable configurations of brain connectivity during rest ( Cavanna et al., 2018 ;Tognoli and Kelso, 2014 ).
While whole-brain and recurrent spatial configurations hubs could not be detected, we have presented results that suggest that the number of active hubs changes periodically over time. Thus, time-periods with a high number of active hubs are interspersed with intervals that are marked by fewer active hubs. These quantitative fluctuations of the number of active hubs were observed at the whole-brain level without any significant changes in the involvement of specific subnetworks. Although no direct links can be made at this point, it is interesting to note that bi-modal properties in TVC have previously been reported for modularity measures in resting-state fMRI ( Betzel et al., 2016 ;Shine et al., 2016 ). However, at a local level, we found recurrent patterns of pairs of nodes that acted has hubs together in time. This observation, particularly noticeable in the Saliency and Subcortical subnetworks, shows that pair-wise co-occurrences happened on average 5-10 ( Fig. 5 A and B, fast hubs) or 15-20 percent ( Fig. 5 C and D,slow hubs) of the total scanning time. This regular pattern of short-lived and pair-wise hub cooccurrences in the functional brain connectome can be viewed as the establishment of highways that serves to provide temporary fast (i.e., short) routes for the flow of connectivity related to increases in neuronal traffic whenever needed ( Fig. 1 C).
Another important observation is that shortest path hubs can be identified at both quicker and faster time-scales. Using the jackknife correlation method, we identified fast hubs located primarily in the medial frontal, and bilateral parietal default mode left frontoparietal and saliency network. In contrast, using the sliding window method that is sensitive to the flow of brain connectivity at slower time-scales, we detected hubs primarily in subcortical structures with no spatial overlap to the hubs found at the faster time-scale.
Our finding that nodes in subcortical regions show a strong preference for slow time-varying changes in hubness is noteworthy. Although subcortical resting-state networks were not targeted in early studies (e.g., Damoiseux et al., 2006 ), it was early on shown that the thalamus is an integrated part of the default mode network ( Fransson, 2005 ). Recent research has confirmed that the thalamus and other subcortical regions such as the nucleus accumbens and regions in the basal forebrain are strongly interconnected with cortical structures in the default mode network ( Alves et al., 2019 ). Interestingly, the role of subcortical structures for cognitive information processing has recently been investigated. Specifically, it has been suggested that subcortical structures play a major role by providing "short-cuts " between lower sensory areas and higher-order cortical areas for certain decision-making processes ( McFayden et al., 2020 ). It is therefore not surprising that subcortical regions frequently act as hubs of brain communication during rest.
Regarding the preference for slower changes in hubness across time for subcortical structures, we note that this finding is in agreement with our previous study on the frequency dependency of resting-state functional connectivity ( Thompson et al., 2015b ). In that study, we created "frequency-graphlets " that essentially use power-spectrum correlations to measure the frequency-dependency on node centrality and global efficiency. In brief, we showed that the subcortical network had a peak of connectivity at around 0.016 Hz when measured across the frequency interval 0.01 to 0.1 Hz. Further investigations have warranted that target the relationship between the temporal properties of subcortical hubs and behavior.
We have in this study showed the feasibility and utility of calculating hubs with flow TVC. In doing so, we have utilized the concept of computing shortest paths on networks that pertain to changes (i.e. flow) of brain connectivity across time. In this context, it is important to point out that shortest paths intrinsically are linked to routing in communication networks since they provide the most efficient use of resources linked to the flow of information in the network. This link makes the implicit assumption that the brain can actively route information selectively through the shortest paths. Whether this assumption is a realistic hypothesis for communication in the brain has recently been debated. For example, it can be argued that the shortest-path route assumption presupposes that neuronal processes have access to knowledge regarding the global brain network topology and that the routes of shortest-paths are dependent on how networks are defined ( Avena-Koenigsberger et al., 2018 ). On the other hand, evidence from studies that combined data from both structural and functional connectivity measures have emphasized the role of shortest paths to ensure reliable and efficient communication ( Goni et al., 2013 ;Hermundstad et al., 2013 ;Bettardini et al., 2017 ). Moreover, the assumption that conservation of energy expenditure is to some degree present in the human brain resonates with previous investigations on the relationship between brain metabolism and cognitive work ( Raichle and Mintun 2006 ).
Interestingly, our approach to measuring shortest paths and BC that penalizes paths that traverse many nodes at the expense of paths that includes few nodes have gained support by previous comparative studies in neuroanatomy and brain connectivity. An example for which empirical evidence has been found for neuronal paths of communication that involve few nodes is favored in nature rather than paths that contain many steps is the well-studied nematode C. elegans. A recently published study used state-of-the-art techniques to investigate the structural connectome in C. elegans at an extremely high spatial resolution ( Cook et al., 2019 ). Importantly, although the authors found that the structural connectome of C. elegans was sparse (3.2% of all possible connections were present), they also found that the connectivity patterns in C. elegans enable information from single sensory neurons to potentially reach from 70% to 98% of all other cells in the network within two synaptic steps ( Cook et al., 2019 ). Their finding that the trafficking of neuronal information between any two neurons to a large extent can be carried out using only two or fewer nodes provides support for our approach to penalizing multi-node paths when calculating the shortest paths.
If we direct our attention to the macaque brain, similar findings have been reported. In the seminal paper on the hierarchical processing of the visual cortex by Felleman and van Essen ( Felleman and van Essen, 1991 ), they showed that neuronal communication in the macaque visual cortex might transcend a hierarchy that spans up to 10 levels. However, they also conclude that the majority of hierarchical levels traverses only 0, 1, or 2 levels. Importantly, the authors argue that neuronal signals can span the entire visual cortex, including all its hierarchies, and using only one relay. Their results show that it is possible to communicate not only within the visual cortex and also to regions in the parietal and frontal areas with very few nodes in between. Admittedly, the Cook study was carried out in a nematode with few neurons compared to the human brain, and the macaque brain is smaller and less complex than the human brain. But we do not see any compelling reasons why this fundamental aspect of brain organization should not have been preserved throughout the evolutionary paths that have formed the structure and function of the central nervous systems in mammals. We find that previous studies performed in different contexts provide experimental evidence and a reasonable biological basis for our idea to approach the problem of finding candidate hubs by calculating distances between brain regions that include a node penalty when calculating the shortest paths.
Networks are models of data. The edges represent some relationships between the nodes. In functional connectivity, Pearson correlations are often used to quantify this relationship. However, when con-sidering multiple time points and analyzing how networks change over time, looking at the difference in edge connectivity instead of the overall statistical relationships helps in understanding time-varying behavior. Flow TVC, as described here, thus differs from the magnitude connectivity as it removes the intrinsic connections between some nodes that always persist. This intrinsic relationship should not be forgotten, but ignoring this long term relationship from the time-varying estimate, helps find which nodes or edges change together regardless of their long term relationship (recall the car lane example from the introduction). It deserves to be mentioned that the concept of time-varying connectivity contains both stable (static) intrinsic connectivity, time-varying connectivity, and noise. Flow TVC aims to remove the static connectivity component. If one can perfectly remove the static component and the noise level is kept constant, then flow TVC will decrease the signal-tonoise ratio of the signal time series. This limitation is an unavoidable consequence of the methodology, which entails that larger sample sizes and replications of identified effects to be highly important in future studies.
It could be argued that the failure to identify discrete patterns could be due to noise driving the BC fluctuations. Moreover, it could be argued that the jackknife correlation method is more sensitive to outliers than the sliding window approach. We, therefore, checked that all jackknife based estimates for betweenness centrality did not correlate with the frame-wise displacement values (see results in Supplementary Figure S14). Moreover, it could be argued that the jackknife method is more sensitive to image scrubbing and interpolation than the sliding window approach. We, therefore, re-analyzed the dataset by leaving out the scrubbing and interpolation step. The results presented in Supplementary Figure S16 shows a pattern of hub activity that is very similar to that obtained with the scrubbing procedure included ( Fig. 2 A-C).
In sum, we could, in our study, not detect discrete brain states during rest based on shortest path lengths and betweenness centrality. However, discrete brain states have previously been shown using TVC fMRI combined with methodological approaches that target other aspects of time-varying brain connectivity. As pointed out earlier, the discrepancy between our results and previous studies of time-varying brain connectivity does not provide contradictory evidence on discrete brain states during rest. It is prudent to keep in mind that the current work is based on a novel method to compute hubs. Research on time-varying properties of hubs in the brain is still in its infancy. The time-varying properties of hubs and their relationship brain states will likely require a simultaneous assessment of various network measures that together will provide additional information.
Nevertheless, from the perspective of the time-varying activity of hubs, our observations suggest that the flow of brain connectivity can be viewed as occurring in a state-space, where each state-vector is related to the shortest paths for neuronal communication. Changes in flow TVC and accordingly, the number and location of active hubs, are carried out in a meta-stable landscape that encompasses basins of attraction that define preferred routes of communication for different spatial networks and time-scales. Our metastability view on the flow of connectivity and the time-varying properties of active hubs also resonates well with a recent EEG study that shows that EEG micro-states are not discrete activities of pools of neurons but rather spatially as well as temporally continuous ( Mishra et al., 2020 ). Moreover, our results might apply to physiologically well-defined transitions, such as different stages of sleep or during anesthesia. A recent study in rats shows that the computational role of hubs changes during REM compared to non-REM sleep ( Clawson et al., 2019 ).

Data availability
Data used in this study is publicly available at https:// db.humanconnectome.org .

Code availability
The code for analyzing the data of study is available at Teneto: https://github.com/wiheto/teneto .

Credit authorship contribution statement
Peter Fransson: Conceptualization, Methodology, Investigation, Software, Analysis, Writing -original draft, Writing -review and editing, funding acquisition William H Thompson: Conceptualization, Methodology, Investigation, Software, Writing -review and editing