On the Spatial Distribution of Temporal Complexity in Resting State and Task Functional MRI

Measuring the temporal complexity of functional MRI (fMRI) time series is one approach to assess how brain activity changes over time. In fact, hemodynamic response of the brain is known to exhibit critical behaviour at the edge between order and disorder. In this study, we aimed to revisit the spatial distribution of temporal complexity in resting state and task fMRI of 100 unrelated subjects from the Human Connectome Project (HCP). First, we compared two common choices of complexity measures, i.e., Hurst exponent and multiscale entropy, and observed a high spatial similarity between them. Second, we considered four tasks in the HCP dataset (Language, Motor, Social, and Working Memory) and found high task-specific complexity, even when the task design was regressed out. For the significance thresholding of brain complexity maps, we used a statistical framework based on graph signal processing that incorporates the structural connectome to develop the null distributions of fMRI complexity. The results suggest that the frontoparietal, dorsal attention, visual, and default mode networks represent stronger complex behaviour than the rest of the brain, irrespective of the task engagement. In sum, the findings support the hypothesis of fMRI temporal complexity as a marker of cognition.

RSNs exhibit a balanced dynamic between order and disorder in time and space, referred to as spatiotemporal complexity. This phenomenon is linked with brain anatomy and has been studied through whole-brain modelling and chaos theory [25,26], as well as are projected onto structurally informed basis vectors and randomization is performed in this joint domain [43].
In this study, we aim to perform an independent assessment on the most commonly reported aspects of temporal complexity in fMRI during rest and task engagement. We use two measures of temporal complexity, i.e., Hurst exponent and multiscale entropy, and compare them in the context of fMRI analysis. First, we validate the monofractal feature of brain hemodynamics during task engagement and rest using Hurst exponent in a population of unrelated subjects from the Human Connectome Project (HCP) [44]. Second, we assess task specificity of fMRI complexity during the task engagement and resting state. To this end, we perform a pair-wise support vector machine (SVM) analysis on the Hurst exponent and a multiscale entropy-based complexity index across all brain regions. Third, we examine the hypothesis of complexity alteration in brain hemodynamics due to task engagement. Fourth, we investigate the agreement between the spatial distribution of Hurst exponent and multiscale entropy in fMRI. Finally, we perform statistical testing on fMRI complexity through a brain structure-informed surrogate technique based on graph signal processing. We look into the mathematical properties of this technique in detail and its consequences for our analysis. Figure 1 illustrates the procedure of performing temporal complexity analysis and graph surrogate generation on a typical fMRI dataset, as adapted in this study.

Data and Preprocessing
We obtained the rest and task fMRI and diffusion-weighted scans of 100 unrelated subjects (ages [22][23][24][25][26][27][28][29][30][31][32][33][34][35] from the HCP1200 release [44]. Each subject underwent a number of fMRI recording sessions including four rsfMRI runs, seven task fMRI runs, and a diffusion MRI run. Each fMRI dataset had a voxel size of 2 × 2 × 2 millimetres and the repetition time (T R ) of 720 milliseconds in a 3-T Siemens Skyra scanner. We utilised two rest runs with left-right phase encoding as well as four task runs with a minimum length of 3 min or 250 T R 's, i.e., Language, Motor, Social, and Working Memory. The other three task recordings, i.e., Gambling, Emotion, and Relational, were shorter than 3 min and therefore excluded from the analysis. Note that each task-based fMRI recording had a specific task design with a different number of conditions and trials (Table 1). See [45] for the description of the fMRI tasks. We included the rsfMRI datasets with left-right phase encoding only due to the known issue of asymmetric drop-out between left-right and right-left phase encoding of rest runs in the HCP database [46] and its potential impact on task specificity analysis in Section 2.3. Since the length of rsfMRI in all subjects was considerably longer than all task fMRI recordings (14.4 min versus 3 to 4 min), we used the first 399 T R 's of rest runs in order to make the extracted complexity measures comparable. The fMRI datasets were preprocessed using SPM8 through a procedure described in [47]. First, fMRI datasets underwent a spatial smoothing by a 5 mm isotropic Gaussian kernel. Six motion parameters as well as average cerebrospinal fluid signal and white matter signal were then regressed out. We did not apply any further bandpass filtering on the fMRI time series at this stage. A parcellation mask [48] was used to parcellate the grey matter into 360 cortical regions of interest (N ROI = 360). The corresponding diffusion-weighted datasets were preprocessed through the steps outlined in [47] and used to extract the structural connectivity matrix of each subject for graph signal processing.

Temporal Complexity Analysis of fMRI
One of the most common ways to measure scale invariance in time series is by using the Hurst exponent H [49]. It determines whether there is a predominant time-scale or frequency component in the underlying dynamical process. The value of H varies between 0 and 1 where H ≤ 0.5 implies short-memory or fast return to the mean (such as white noise), and H ≥ 0.5 represents long-memory or a trending behaviour with random turning points. The value of H = 0.5 represents a random walk whose time points have no correlation with their past values. Scale-free signals have a long memory, because all of their time scales and spectral components contribute equally to their dynamics. This leads to a power law relationship in the spectral power of scale-free signals in the form of P( f ) ∝ f −β where P( f ) is the power spectral density at frequency f , and β is a non-zero positive real number referred to as spectral exponent. For some scale-free processes such the fractional Brownian motion, there is a theoretical relationship between the spectral exponent β and the Hurst exponent H via the equation β = 2H + 1 [50]. In this study, we used detrended fluctuation analysis (DFA) [51] to estimate the Hurst exponent and log-linear line fitting to the power spectral density function to estimate the spectral exponent of fMRI time series. The DFA algorithm has been widely used in the past for the analysis of fMRI fractality [30,52]. Previous studies have reported an approximate range of 0.5-0.9 for the Hurst exponent of fMRI time series which highlights them as a class of signals with complex dynamic and long memory [30,52]. Given the direct link between the Hurst exponent and signal entropy [53], we also looked into the complex behaviour of fMRI using multiscale entropy [54]. This measure is based on the sample entropy [55] at several time scales of a signal x (here, mean fMRI time series at a particular ROI). Multiscale entropy characterises white noise by a very large entropy value (usually above 4) at the first time scale which rapidly decreases across coarser time scales. This leads to small areas under the multiscale entropy curve for highly irregular signals such as white noise. On the other hand, complex signals such as red/pink noise represent a relatively flat or monotonically increasing pattern across most of the time scales with a larger area under the curve compared to white noise [14]. Therefore, the area under the curve of multiscale entropy can be considered as a complexity index for temporal complexity analysis of biosignals [56]. See Appendix A for a detailed description of multiscale entropy.

Task Specificity of fMRI Complexity
In order to test the dependency of the complex properties of brain hemodynamics on task engagement, we used a set of SVMs with the radial basis function (RBF) kernels in order to evaluate the separability of the Hurst exponent and multiscale entropy over two rest runs and four task runs of fMRI. The SVM analysis is a supervised learning method which is widely used for performing classification and regression studies using fMRI datasets. We extracted ROI-wise complexity measures of N subj = 100 subjects and six fMRI runs. For each pair-wise comparison between the two fMRI runs, we considered each subject as an observation and the brain maps as feature vectors of size N ROI × 1 where N ROI = 360. This yielded 100 feature vectors of size N ROI × 1 for each fMRI run. Here, we investigated whether mental tasks can be discriminated in a population of healthy subjects using the ROIwise temporal complexity of fMRI. We quantified the performance of SVM classifiers using the percentage of their classification loss or the proportion of observations misclassified by the model. This pair-wise SVM analysis led to two symmetric accuracy matrices of size 6 × 6 for the two temporal complexity measures. The SVM classifiers were evaluated via 5-fold cross-validation and their hyper-parameters were optimized through grid search. We minimized the issue of overfitting in the SVM classification analysis by optimizing the box constraint parameter in MATLAB's fitcsvm function. This parameter controls the maximum number of support vectors, thereby preventing overfitting. The larger the box constraint parameter, the less support vectors are assigned to the SVM classifier.

Spatial Distribution of fMRI Complexity across Grey Matter
To be able to interpret the complexity measures of fMRI across brain regions, one must perform statistical testing. In this study, we generated N Surr = 100 surrogate fMRI datasets for each subject and each fMRI run through a graph signal processing framework which combines the brain structure and function in order to generate the surrogates of fMRI whose null distribution preserves the spatial smoothness of the fMRI signal on the structural connectome at each time point [42,43]. Our motivation for adapting this technique was to incorporate the underlying anatomical aspects of fMRI in the significance testing step, an important piece of information which is usually neglected in the fMRI complexity analysis studies. As discussed in Appendix B, the graph surrogate method: (i) preserves second order statistical moments across fMRI time points, i.e., temporal correlation; (ii) randomizes functional connectivity; (iii) randomizes the spatial variation of scale-free dynamics of fMRI at each single ROI; and (iv) randomizes the spatial variation of the scale-free dynamics between ROI pairs.
In order to obtain the group-level maps of brain complexity, we applied binomial testing on the subject-specific brain maps. First, each individual map of complexity was thresholded at a significance level of α subj = 0.01 in order to obtain a binary map at the subject level. Then, the binomial distribution P(n) of having n detections was used at each ROI to examine the significant number of suprathreshold regions across subjects at the significance level of α group = 0.001. It was equivalent with 10 detections (i.e., the number of suprathreshold regions) for a population of 100 subjects. The results were corrected for multiple comparisons at the number of regions and fMRI runs tested, i.e., 360 × 6 comparisons.

FMRI Represents Complex Behaviour during Rest and Task
As seen in Figure 2, the normalized group mean power spectral density of all fMRI recording sessions (four fMRI runs and four task runs) show log-linearity in the frequency domain. This feature was more pronounced in rsfMRI datasets and the spectral exponents were RSN-specific with default mode, frontoparietal, and dorsal attention networks presenting the highest exponents in most cases ( Figure 2B). The β exponents of different RSNs, however, were shown to be task-dependent. For example, the Social task led to the highest spectral exponent across subjects at the dorsal attention network or Rest runs led to the highest exponents at the default mode and frontoparietal networks. A striking observation was related to the existence of dominant peaks in the log-linear power spectral density functions of task fMRI in contrast to rsfMRI (Figure 2A). In order to rule out the influence of the task designs in this spectral feature of task fMRI recordings ( Figure 2C), we regressed out the task timings trial-wise from the data and checked the spectral log-linearity of the residuals only. As Figure 2D illustrates, the peaks still remain in the log-linear power spectral density functions of task fMRI residuals, although they are slightly suppressed. This suggests that these spectral peaks are independent from the task design and likely explain why the Hurst exponent is generally smaller during task engagement and larger at rest [21,30,39].

Task Engagement Lowers Complexity of BOLD Activity
As summarized in the classification loss matrices of Figure 3, the fMRI of resting state and task engagement can be classified with high accuracy using both temporal complexity measures. However, the dynamics of fMRI sessions at rest are not distinguishable from each other (note the high classification loss between rsfMRI recordings in Figure 3C,D). Here, each dataset has been characterized as a set of N subj feature vectors of size N ROI × 1 where N subj = 100 and N ROI = 360. The distribution of temporal complexity across brain areas in different rest and task fMRI runs (brain maps of Figures 3A and 4) suggests that the spatial profile of complex dynamics in fMRI is affected by task engagement and resting conditions. As evident in the histograms of individualized mean Hurst exponents across brain areas in Figure 3B, each task results in a regionally specific reduction in dynamic complexity, an observation in line with previous findings in the literature [30]. This establishes a framework for comparing the task burden in subjects. For example, the Hurst exponent histograms suggest that the working memory task is likely the most demanding brain state in the population of this study, because its associated histogram covers the lowest interval of self-similarity.

Complex Dynamics Exist in the Brain Structural-Functional Coupling
In Figure 5, we show the spectral power and multiscale entropy patterns of the projected fMRI datasets onto brain structure at rest and during task engagement. As Figure 5A illustrates, the distribution of spatial energy across the brain connectome follows a log-linear relationship. However, despite what we observed for the spectral power of fMRI in Figure 2A, a considerable difference between the power spectral distributions of graph signals extracted from different tasks and rest sessions is not evident here. This can also be observed in the multiscale entropy patterns and complexity indices of brain graph signals in Figure 5B,C which are quite comparable over four tasks and the average rest. In particular, there is an inflection point in the scatter plot of complexity indices ( Figure 5C) associated with the colour transition in the multiscale entropy patterns of Figure 5B from dark blue (low randomness) to green (high randomness). These elbow or knee points are also indicative of log-linearity in the complexity domain of brain structure-function.
Complex dynamics of brain graph signals Figure 5. (A) Logarithmic plots of the power spectral density functions of brain graph signals (i.e., the projection of the fMRI data at rest and task onto brain structure) versus brain spatial harmonics. The plots of two rest runs are averaged. Each grey curve belongs to a single subject and the red curves represent group mean. All curves are normalized to 1. (B) Multiscale entropy patterns of the graph signals, colour coded by their associated brain spatial harmonyic (C) The complexity indices associated with the multiscale entropy curves of (B). Figure 6A illustrates the group-mean spatial patterns of multiscale entropy-based complexity index across brain areas for the average rest runs and four tasks. Furthermore, Figure 6B presents the pie chars of the percentage of suprathreshold ROIs in 7 RSNs associated with Figure 6A. All maps were thresholded at the subject-level p-value of 0.01 and family-wise error-corrected at the group-level p-value of 0.01 using the graph surrogate data generation method introduced in [43]. Furthermore, Table 2 summarizes the contribution of 7 RSNs in the suprathreshold brain regions of Figure 6. In all fMRI runs, the visual areas were amongst the regions with highest complex dynamics. The next mostly engaged RSNs in all runs were the dorsal attention network with a maximum cover of 7.5% followed by the frontoparietal and default mode networks with a maximum cover of 8.9%, all during the Working Memory task. Several regions across dorsal lateral prefrontal cortex (DLPFC) showed up as the most frequently observed areas with the highest complexity across all tasks and rest runs. These regions included 8Av, 9m, 8BL, 10d, 9a and p10p. Other most repeated suprathreshold regions included: left RSC (anterior cingulate cortex), left TE1a (middle temporal gyrus), PGi (temporo-parieto occipital junction), PGs (inferior parietal cortex) and right 45 (inferior frontal cortex). The limbic network represented the least complex regions with most random behaviour in all task and rest runs. See [48] for the anatomical description of these brain labels.

Spatial Patterns of Complex Dynamics in fMRI
Thresholded fMRI complexity maps using the graph surrogate data method Figure 6. (A) Spatial distribution of group-mean fMRI temporal complexity across brain areas for 4 task runs and the average rest run. All maps are thresholded using the graph surrogate data generation [43] at the subject level p-value of 0.01 and family-wise error corrected at the p-value of 0.01. (B) Pie charts are the percentage of suprathreshold ROIs in 7 RSNs after graph surrogate testing, normalized by the number of ROIs. See Table 2

Discussion
This study reinforces the existence of complex dynamics in brain function [22] and provides further evidence for the hypothesis of distinct complexity features in human behaviour and cognition [30]. Our results suggest that: (i) task-based and rsfMRI signals exhibit temporal complexity, inferred by Hurst exponent and multiscale entropy; (ii) rest and task periods of brain function can be distinguished from each other with high accuracy based on their temporal complexity profiles; (iii) cognitive load can suppress the complex dynamics of fMRI in contrast to the resting state; (iv) spatial distribution of Hurst exponent and entropy-based complexity index in fMRI are highly correlated; and (v) the visual, frontoparietal, default mode, and dorsal attention networks represent maximal complex behaviour compared to the rest of the brain in most mental states.
Temporal changes of neural activity in the brain are observed on the scale of milliseconds in single-cell spiking to the order of seconds in BOLD fluctuations [57]. In contrast to the traditional neuroscientific view which would mostly consider this variability as a random disturbance and measurement noise, many systematic patterns of information have been detected in neural variability over multiple temporal and spatial scales [54,58,59]. Neural variability at the BOLD level not only reflects inter-subject differences such as behavioural traits, but it also captures within-subject changes such as the dynamical states of the brain and cognitive load [57]. An important feature of neural variability in the brain is its temporally complex behaviour. Although temporal variability and temporal complexity are closely related, these two concepts are not necessarily the same. In other words, temporal complexity always exhibits variability, but a variable process is not necessarily complex. Also, high complexity is not necessarily equivant with high entropy. As a matter of fact, Gaussian white noise represents the highest signal entropy, highest variability, and highest unpredictability amongst all signal types, but it is not deemed as temporally complex according to our definition. A complex pattern of brain activity is rich in information over time and represents a balanced dynamic between order and disorder within brain networks or between different brain regions [13,14]. Recent theories in neuroscience have proposed that the temporal complexity of brain function is likely associated with information processing in the brain [13,59]. The first theory [60] indicates that a healthy brain retains an optimal degree of instability which allows it to enter to different dynamical states and sample the external world. In this way, the brain learns how to optimize responses to environmental stimuli. The second set of theories [61] argue that a moderate level of randomness is necessary in the neural system because it increases the probability of neuronal firing in subthreshold neurons. In contrast to the second proposal, a third view [62] is that the temporal complexity of brain function can enhance or suppress the likelihood of neural synchrony between cortical areas. All of these theories speak to the direct relationship between the temporal complexity of brain function and information transfer across neural populations, cortical regions and functional networks.
Considering the key role of complex processes in brain mechanisms, one would expect to find realizations of temporal complexity in brain hemodynamics as well. To investigate this possibility, it is crucial to quantify the temporal complexity of neural variability. Signal entropy measures were found to be relevant for investigating the complexity of brain function. This is mainly because these measures do not rely on distributional assumptions (unlike variance-based measures), and are not restricted by sinusoidal signal waveforms (unlike frequency-based measures). In this study, we utilized two measures for the temporal complexity analysis of fMRI in order to cross-check the fMRI complexity analysis results: multiscale entropy as an entropy-based measure which evaluates similar patterns of information throughout signals and the Hurst exponent as a variance-based measure which quantifies the memory of signals and their tendency to return to their mean values. Multiscale entropy has been shown to be sensitive to RSN-specific hemodynamics, and reproducible across healthy subjects [13,14]. The direct relationship between multiscale entropy and self-similarity [53] makes it a good candidate for investigating the properties of the Hurst exponent in fMRI.
It is important to note that the existence of self-similarity and fractality in BOLD signals has been subject to discussion in recent decades, mainly due to the limited temporal resolution of fMRI and sluggishness of the hemodynamic responses in the brain [63]. A fundamental question here is, even assuming the presence of temporal complexity in brain hemodynamics, that of whether the recorded BOLD signal inside an MRI scanner is able to adequately capture it. To address this question, the following issues should be considered: (i) low temporal resolution of the measured BOLD changes with a typical T R value in the scale of hundred milliseconds to seconds; (ii) the use of short-length fMRI time series (less than 200 T R 's) in some studies; (iii) potential impact of preprocessing steps and residual scanning artifacts on the nonlinear dynamics of fMRI; (iv) upon the agreement on the complex dynamics of fMRI, whether it is monofractal or multifractal [21,64]. It has been shown that non-fractal time series with inadequate memory can still exhibit log-linear spectral power and be falsely identified as complex processes [63]. Therefore, one may argue that the observed complex behaviour of fMRI is not biological, but it simply originates from the signal aspects of fMRI such as its inadequate length in the previous studies. In fact, it has been shown that the accuracy of fMRI fractal analysis can be affected by the number of brain volumes [64]. As Figure 4 shows, a linear association across brain regions was evident between the Hurst exponent and the area under multiscale entropy curves of rsfMRI at the group level. Although the two measures operate at different ranges and utilize different methodologies to quantify temporal complexity, their spatial agreement across brain regions is relatively high (Pearson correlation above 0.8). However, the distinction between the dynamics of rsfMRI and task fMRI is more apparent in the Hurst exponent brain maps (Figure 3) in contrast to the entropy-based complexity index brain maps (Figure 4). It speaks to numerous monofractal analysis techniques which one can utilise to estimate the temporal complexity of fMRI, though each method may vary in its strengths, biases, and sensitivity to the biological changes of brain function. In spite of these technical differences, several simulated and experimental studies supported the hypothesis of a power law distribution for the fMRI power spectrum over the frequency band of 0.01 Hz to 0.1 Hz [21,52,65,66]. See [67] for a recent review on this topic.
The magnitude of H across brain areas in our study ( Figure 3) is comparable with the previously reported findings showing a typical range of ≈0.5-1 for H [30,68,69]. The spatial extent of fMRI temporal complexity in our results is maximal across the frontoparietal, dorsal attention, visual and default mode networks and minimal across deep brain areas such as the limbic network (see Figures 3 and 4). The considerable overlap between the analysis results of this study with the existing literature suggests that the fMRI signal length (minimum of 263 T R 's-see Table 1) has been enough to replicate the previous hypotheses about the temporal complexity of brain function. In fact, a previous study has shown that valid Hurst exponents can be obtained from fMRI time series as brief as 40 s (≈56 T R 's) through the DFA algorithm [30]. Furthermore, this implies that the fMRI preprocessing steps have not significantly manipulated the true dynamics of fMRI datasets in this study. According to our results, the spatial patterns of the Hurst exponent in fMRI are highly correlated with the associated patterns of entropy-based complexity index. This indicates that the two measures of fMRI temporal complexity converge to a similar outcome even though their computation is completely different.
The results of this study suggest that the temporal complexity of some regions such as motor areas is task-dependant. Higher values of complexity index and Hurst exponent of a given brain region or RSN (Figure 7) imply that the associated fMRI signals are temporally redundant and more predictable. Regions/networks with slower dynamics are likely responsible for the process of internal stimuli with low surprise and high adaptability. However, external stimuli such as sensory inputs and auditory inputs may reduce the adaptability of brain dynamics, increase surprise and shift the temporal complexity of brain function towards faster dynamics. An exception would be the visual network which shows variable dynamics from slow to fast across different tasks (see Figure 7), likely due to the larger capacity of visual areas in contrast to the other networks and brain regions. It is important to note that despite the presence of complex dynamics in fMRI, the relationship between fMRI temporal complexity and intrinsic functional connectivity is scale-dependant. It has been shown that the weighted sum of functional links to a given brain node, referred to as functional connectivity strength or FCS, is associated with the functional significance of that node in support of the information transfer across the brain [70]. The link between the FCS and temporal complexity of rsfMRI has been shown to be scale-dependant [13,14]. In particular, an inverse relationship has been reported between the FCS and temporal complexity of RSNs at fine time scales of multiscale entropy (τ ≤ 5 at a T R of 0.72 s equivalent to time periods shorter than 3.5-4 s), while it turns to a proportional relationship at coarse time scales (τ > 6 or time periods greater than ≈4 s). This scale-dependant relationship varies for different RSNs. For example, frontoparietal and default mode networks represent the highest correlation between resting state FCS and temporal complexity at fine scales and lowest correlation at coarse scales, while it is the opposite for somatomotor, sub-cortical and visual networks [13,14]. On the other hand, the fine time scales of fMRI were associated with the dynamics of local neural populations, whilst the coarse time scales are likely related to long-range functional connections [13]. Altogether, these observations would suggest that, in order to obtain a comprehensive picture about complex dynamics of fMRI, one must consider the anatomical locations (i.e., spatial distribution) of brain regions. This speaks to the necessity of appropriate spatiotemporal methods for the significance testing of fMRI temporal complexity which can account for brain structure and function at the same time.
Surrogate data testing is a powerful method for characterizing the statistical properties of time series. In this approach, we want to compare the measure of interest extracted from the original data, i.e., the alternative hypothesis H 1 , to the distribution of the same measure obtained from a large number of surrogate data, i.e., the null hypothesis H 0 . In the context of fMRI temporal complexity, the null hypothesis could be that fMRI time series at different brain regions are generated by some non-complex processes and their functional relationships are also non-complex. If the complexity indices of fMRI signals fall within the null distribution H 0 , this means that these indices only rely on the statistical properties preserved under the null H 0 . Otherwise, we can reject the null hypothesis and interpret fMRI temporal complexity as revealing statistical properties beyond H 0 . A critical question here is how to specify the null hypothesis of scale-free dynamics in fMRI and how to remove this signal feature of interest from the original data in order to generate surrogates [71,72]. In this study, we chose the brain graph randomization technique [43] which shuffles the power spectral density of fMRI in the graph domain for each time point, while keeping its underlying anatomical properties. As analytically discussed in the Appendix B, the surrogate data in this study preserve the pair-wise correlation between fMRI time series between brain regions. However, the functional connectivity, spatial variation of complex dynamics and spatial variation of the complex dynamics of functional connectivity are randomized in the surrogate time series. Therefore, the complex dynamics of brain ROIs with suprathreshold complexity indices is significantly different from the chance level and likely play a key role in moderating the information transfer across the whole brain. According to Table 2 as well as Section 3.4, the supporting role of brain regions is maximal within the frontoparietal, dorsal attention, visual, and default mode networks such as DLPFC and PGi, regardless of the resting state or task engagement. This finding is in line with the previous studies showing the dominant role of these brain areas in internal processing and the facilitation of information exchange in functional brain networks [73,74]. A number of caveats need to be considered in interpreting the results of this study. From the technical perspective, one should be aware of the limitation of Hurst exponent and multiscale entropy in capturing multivariate relations between fMRI time series at multiple ROIs. Both complexity measures used in this study treat ROI-wise fMRI signals as a set of individual and independent time series. However, mean fMRI time series at different ROIs are often statistically dependent and correlated due to the smearing effect of hemodynamic changes in the brain and a possible effect of fMRI preprocessing steps such as spatial smoothing. Therefore, it is plausible to use the multivariate versions of temporal complexity measures for fMRI data analysis [75,76]. A systematic comparison between the univariate and multivariate measures of complex dynamics and temporal complexity remains for our future work. Another consideration should be given to the nature of BOLD signals as an indirect measure of neural activity and the reduced amount of information it carries, in contrast to other direct measurements with higher temporal resolution such as local field potentials. Although fMRI has a greater capacity to cover larger brain areas than the localized measurements of neural activity and can provide large-scale information about the complex properties of brain dynamics, its temporal complexity must be treated as an indirect property of brain function. In this study, we assumed monofractality in fMRI and compared the entropy-based complexity index of fMRI with the classical Hurst exponent at different ROIs. It would be informative to check the possible links between the multiscale entropy and multifractality of fMRI using longer fMRI datasets and higher temporal resolutions during rest and task engagement. We also focused on the complexity of fMRI and its spatial distributions in the time domain only. However, this is only one dimension along which the complexity of brain dynamics can be measured. A more holistic view would be to consider the complexity of brain function in both time and space and adapt multivariate measures for the complexity analysis of fMRI time series [77].

Conclusions
Temporal complexity is a reproducible aspect of fMRI during rest and task engagement. This feature of brain function is task specific and can be suppressed by cognitive load. FMRI complexity is a discriminative feature between rest and task in the brain functional domain, but not in the brain structural domain. A brain structure-informed statistical testing of fMRI complexity reveals several areas with suprathreshold temporal complexity within the frontoparietal, visual, dorsal attention, and default mode networks.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. Multiscale Entropy Analysis
Multiscale entropy analysis [54] is based on the calculation of sample entropy [55] at multiple time scales. Sample entropy treats each short piece of x as a template to quantify a conditional probability that two templates of length m, which are similar to within a tolerance level r, will remain similar when m becomes m + 1. Note that self-matches are not considered in calculating this conditional probability. A template X m i is defined as (in all equations, scalar variables are in normal font, while vector variables are in bold.): where N is the number of time points in x and m is the embedding dimension. Two templates X m i and X m j are considered as neighbours if their Chebyshev distance d(X m i , X m j ) is less than a tolerance parameter r. This leads to an r-neighbourhood conditional probability function C m i (r) for any vector X m i in the m-dimensional reconstructed phase space: where B m i (r) is given by: where Ψ(.) is the sign function with 0 for non-positive inputs and 1 for positive inputs. Sample entropy is then given by: where B r m is the average of B m i (r)'s over all templates: Sample entropy is always non-negative, but it can also become undefined. It is important to multiply the tolerance parameter r by the standard deviation of x to account for amplitude variations across different signals [55]. Multiscale entropy extracts sample entropy after the coarse graining of the input signal x at a range of time scales τ. A coarse-grained vector x(τ) = {x i (τ)} is defined as: In this study, we used these parameter values: m = 2, r = 0.5 and τ max = 10. We reduced the dimensionality of multiscale entropy patterns to a single value by calculating the area under each multi-scale entropy curve over all scales, divided by the maximum number of scales (i.e., τ max ). This leads to a complexity index M i for each parcellated rsfMRI dataset as follows [56]:

Appendix B. The Graph Surrogate Method for fMRI Complexity Analysis
In this section, we explain a surrogate data method, proposed in [43], which we hypothesize is a relevant technique for generating null distributions in the hypothesis testing of fMRI complexity. This approach uses the graph signal processing framework [42] in order to combine the brain structure and function. Let Y ∈ R N ROI ×N T be a parcellated rsfMRI dataset from N ROI brain regions of interest (ROIs) which were uniformly sampled at N T time points.

Appendix B.1. Combining Brain Structure and Function
In order to model the spatiotemporal correlates of fMRI, we consider a basis set U G ∈ R N ROI ×N ROI which contains the spatial harmonics of brain function and another basis set U F ∈ R N T ×N T which models the dynamics of fMRI using a Fourier transform with N T bins. Joint projection of Y into the structural and spectral domains of the brain can be modelled as: whereŶ is a N ROI × N T matrix and g is a nonlinear mapping from the time domain into the joint complexity-structure domain. In this study, we assume that g is a linear function and also, the structural and spectral domains of brain function are linearly separable. In that case, Equation (A8) is written as:Ŷ where U G is the transpose of U G and * denotes the complex conjugate operator. This general formulation considers the spatial and spectral properties of brain function through U G and U F , respectively. If U F is an identity matrix, Equation (A9) leads to the representation of graph signal processing for fMRI datasets.

Appendix B.2. Spatial Harmonics of Brain Structure
In order to obtain the structural basis set of brain function (i.e., U G in Equation (A9)), one can incorporate a brain structural graph from the corresponding diffusion MRI [42]. This graph can be characterized as G = (V, E) where V is a set of N ROI vertices and E is a set of weighted edges associated with a symmetric and real-valued adjacency matrix A G ∈ R N ROI ×N ROI . Each element in A G represents the number of white-matter pathways between two brain regions. The symmetric normalized Laplacian matrix of A G is then used to form a structural space onto which the brain function is projected: where I ∈ R N ROI ×N ROI is the identity matrix and D G ∈ R N ROI ×N ROI is the degree matrix of A G , a diagonal matrix whose non-zero elements are defined as: where a ik is the ik'th element of A G . According to the spectral graph theory [42,78], spatial harmonics of the fMRI temporal correlation matrix G can be obtained through the eigendecomposition of L sym G as follows: where the diagonal matrix Λ G includes the eigenvalues (or spatial harmonics) of L sym G and U G contains the corresponding spatial eigenmodes. Given that A G is real valued and symmetric, both Λ G and U G will be real valued. The matrix U G constitutes an N ROIdimensional brain structural space or structural space, in short. The assumption here is that brain structure is fixed in time. Therefore, all elements of U G are time invariant. The left-side multiplication of U T G Y in Equation (A9) is usually referred to as graph signals in the literature [42].

Appendix B.3. Graph Surrogate Generation
We used the graph surrogate data generation method in [43] to develop the null distributions of fMRI complexity. This technique is based on the idea of sign randomization in brain structural eigenmodes. Let A ∈ R N ROI ×N ROI be the fMRI temporal correlation matrix associated with the vectorized fMRI data Y ∈ R N ROI ×N T . The spatial eigenmodes of A could be obtained through graph signal processing as the columns of the square and orthogonal matrix U G ∈ R N ROI ×N ROI .
One can randomize the eigen matrix U G and generateŨ G by keeping the absolute value of its elements, while shuffling the signs of eigenvectors (i.e., the signs of all elements across a column are either flipped or not). In this case, a surrogate form of Y with a similar size, called hereafterỸ, is obtained as [43]: Each element ofB can be written as: where u ik is the ik'th element of U G andũ kj is the kj'th element ofŨ G . In Equation (A13), U T G Y randomizes the fMRI matrix Y in the structural brain graph domain and left-side multiplication to U G takes the randomized graph signals back to the functional domain. Note that, in the case of U G =Ũ G (i.e., if there is no shuffling in the graph domain), while in the case of graph shuffling, we have: where a ∈ R is a real number.

Appendix B.4. Regarding Linearity
The surrogate data matrixX in Equation (A13) can be also written as: whereb i is the i'th column inB. By assuming the index j as a temporal argument in Equation (A17), the i'th element ofX at time t is given by: This implies that each element of the surrogate dataX at time t is a linear random combination of its corresponding fMRI time point over all ROIs. Therefore, the graph surrogate method [43] preserves the linear relationships between original fMRI time points.

Appendix B.6. Regarding the fMRI Temporal Correlation Matrix
To unravel the pair-wise relationship between the time points of graph surrogates, a similar cross-correlation analysis can be performed by simply transposingX in Equation (A19). This leads to a super adjacency matrixC T ∈ R N T ×N T as follows: This is based on the fact that both U G andŨ G are orthogonal, so the termsŨ G U T G and U GŨ T G are identity matrices. SinceC T is independent fromX, it does not change throughout the surrogate generation procedure. This suggests that the fMRI temporal correlation matrix embedded in the surrogate data matrixX is preserved through the method described in [43].
Appendix B.7. Regarding Temporal Complexity The complex properties of a multivariate time series can be studied by looking into the log-linearity of its auto-and cross-spectral power functions in the frequency domain. A log-linear pattern of spectral power means that the frequency components have been distributed according to an exponential law S( f ) = 1/ f β where β is referred to as the spectral exponent. This parameter is related to the memory of the underlying signal measured by Hurst exponent H, varying between 0 and 1. A value near 0 reflects the short memory (i.e., quickly coming back to the signal mean or noise-like behaviour), while a value close to 1 represents long memory (i.e., slow return to the signal mean or a smooth dynamic). In the spatial case of fractional Brownian motion, an exponent of H = 0.5 leads to a random walk.
In order to check the complex properties of graph surrogates inX, one can use the linear expansion in Equation (A18) and obtain the Fourier transform ofx i (t) as follows: where x k (t) is the mean fMRI time series at the k'th ROI. The power spectral density function ofx i (t) is then obtained as: Each single multiplication term F m ( f )F * n ( f ) in Equation (A22) is either equal to the power spectral density of individual surrogate time series (i.e., S mm ( f )) or coherency between pairs of surrogate signals (i.e., S mn ( f )) as follows: Therefore, the spectral power of each surrogate time series inX is a random linear combination of all spectral powers and all possible pair-wise coherencies between original fMRI time series across ROIs. Now, let us assume that there is some level of complex dynamics within the original fMRI time series, that is, the spectral power and coherency across a sub-group or all temporal dimensions of X follow the power law: S mn ( f ) ≈ 1/ f β mn , m = 1, ..., N ROI , n = 1, ..., N ROI .
With this assumption, it is clear from Equation (A22) that the power spectral density functions of surrogate dataX (i.e.,S ii ( f ), i = 1, ..., N ROI ) do not necessarily preserve the complex properties of the original fMRI data at individual ROIs or interaction between regions. If there is no variation in complex properties across brain regions (i.e., if the spatial distribution of spectral exponents is the same across all regions), the Hurst exponent of surrogates would be the same as the original data. This is an example in which the complex properties of the original dataset are completely preserved by the graph surrogate data method [43]. In practice, however, the fMRI temporal correlation matrix A imposes a non-uniform distribution of the complex properties across regions which can be picked up by the graph surrogate method. Note that the coherency between different graph surrogates (i.e.,S ij ( f ), i, j = 1, ..., N ROI , i = j) gives a shuffled version of functional connectivity in the time domain. To show this, we take advantage of the cross-correlation and Wiener-Khinchin theorems outlining the reciprocal relationship between correlation and coherence in the time and frequency domains. The Wiener-Khinchin theorem states that forx i (t), the power spectral density functionS ii ( f ) is equal to the Fourier transform of its auto-correlation functionR ii (τ): where τ denotes the delay and F is the Fourier transform operator. Furthermore, the cross-correlation theorem indicates that the correlation between surrogatesx i (t) andx j (t) in the time domain is equivalent to their coherence in the frequency domain: whereR ij (τ) is the cross-correlation ofx i (t) withx j (t) at the delay τ and * denotes the complex conjugate operator. This suggests that the graph surrogate method [43] tends to remove the spatial variation of possibly complex properties of the original data at each single fMRI time series, which preserves the temporal relationships but randomizes their functional connectivity.