# Multifractal Dynamic Functional Connectivity in the Resting-State Brain

- Department of Physiology, Semmelweis University, Budapest, Hungary

Assessing the functional connectivity (FC) of the brain has proven valuable in enhancing our understanding of brain function. Recent developments in the field demonstrated that FC fluctuates even in the resting state, which has not been taken into account by the widely applied static approaches introduced earlier. In a recent study using functional near-infrared spectroscopy (fNIRS) global dynamic functional connectivity (DFC) has also been found to fluctuate according to scale-free i.e., fractal dynamics evidencing the true multifractal (MF) nature of DFC in the human prefrontal cortex. Expanding on these findings, we performed electroencephalography (EEG) measurements in 14 regions over the whole cortex of 24 healthy, young adult subjects in eyes open (EO) and eyes closed (EC) states. We applied dynamic graph theoretical analysis to capture DFC by computing the pairwise time-dependent synchronization between brain regions and subsequently calculating the following dynamic graph topological measures: Density, Clustering Coefficient, and Efficiency. We characterized the dynamic nature of these global network metrics as well as local individual connections in the networks using focus-based multifractal time series analysis in all traditional EEG frequency bands. Global network topological measures were found fluctuating–albeit at different extent–according to true multifractal nature in all frequency bands. Moreover, the monofractal Hurst exponent was found higher during EC than EO in the alpha and beta bands. Individual connections showed a characteristic topology in their fractal properties, with higher autocorrelation owing to short-distance connections–especially those in the frontal and pre-frontal cortex–while long-distance connections linking the occipital to the frontal and pre-frontal areas expressed lower values. The same topology was found with connection-wise multifractality in all but delta band connections, where the very opposite pattern appeared. This resulted in a positive correlation between global autocorrelation and connection-wise multifractality in the higher frequency bands, while a strong anticorrelation in the delta band. The proposed analytical tools allow for capturing the fine details of functional connectivity dynamics that are evidently present in DFC, with the presented results implying that multifractality is indeed an inherent property of both global and local DFC.

## Introduction

Functions of the brain emerge from dynamic interactions between the elements of its complex neuronal networks (Chialvo, 2010; Werner, 2010). This phenomenon is present across a broad range of spatial scales from the microanatomical level of individual neurons through neuronal cell assemblies to macroanatomical brain regions (Sporns, 2011). Pioneering works by Friston et al. (1993) and Biswal et al. (1995) paved the way to the emergence of a new field of neuroscience aiming at describing brain function through its anatomical and functional connectivity (FC) (Sporns et al., 2005; van den Heuvel and Hulshoff Pol, 2010; Friston, 2011). The key concept underlying the latter is assessing the statistical interdependence of neural activity registered at disparate regions of the brain, as it is assumed to be proportional to the degree of their functional cooperation (Friston et al., 1993). FC studies did not only reveal the existence of several resting-state brain networks such as the default mode network (Raichle et al., 2001; Greicius et al., 2003) or the task-positive network (Fox et al., 2005), but also showed that FC properties responded to changes in physiological conditions e.g., sleep (Horovitz et al., 2009; Liu et al., 2015) or cognitive stimulation (Esposito et al., 2006; Racz et al., 2017). In addition, altered FC was also found in pathological conditions like those in degenerative dementias (Pievani et al., 2011), schizophrenia (Liu et al., 2008), or multiple sclerosis (Cader et al., 2006).

Until recently, most studies typically considered connections in functional networks and thus FC itself, too, as being stationary despite the fact that a dynamic approach might provide a more detailed and more realistic description of brain connectivity (Hutchison et al., 2013a). Indeed it has been shown that FC dynamically fluctuates even in the resting state (Chang and Glover, 2010) and also during task modulation (Sakoglu et al., 2010). Since then, investigating dynamic functional connectivity (DFC) has become one of the most rapidly evolving fields of neuroscience with a steadily expanding body of literature (Hutchison et al., 2013a; Calhoun et al., 2014; Preti et al., 2017).

In FC studies most often functional magnetic resonance imaging (fMRI) is used to monitor brain activity with high spatial resolution and precise anatomical localization (Hutchison et al., 2013a; Preti et al., 2017). Fluctuations in FC are usually captured with a sliding-window (SW) approach, however other approaches such as point process analysis (Tagliazucchi et al., 2012a) or paradigm free mapping (Gaudes et al., 2013) have also been presented. During SW analysis, FC is calculated from a small data segment (i.e., within the actual window), then the window is advanced by a predefined time step and the process is repeated until the whole signal is covered. To assess FC within the actual window, statistical interdependence is usually estimated by bivariate statistical methods as Pearson-correlation (Hutchison et al., 2013b), but multivariate methods such as spatial independent component analysis (Allen et al., 2014) or time-frequency methods (Chang and Glover, 2010) are also often used. DFC is then described through–including but not limited to–properties such as the number of stable global states, their variability, and transition probabilities (Allen et al., 2014; Calhoun et al., 2014; Damaraju et al., 2014). As graph theory provides a useful tool in characterizing complex networks of the brain (Bullmore and Sporns, 2009) along several topological aspects (Rubinov and Sporns, 2010), dynamic graph theoretical analysis is also frequently applied (Tagliazucchi et al., 2012b; Yu et al., 2015; Du et al., 2016; Racz et al., 2018). Finally, some DFC studies focus–instead of on global network topology–only on one or a few individual connections between specific regions (Rack-Gomer and Liu, 2012) or intrinsic connectivity networks (Chang and Glover, 2010; Allen et al., 2014). The fluctuating nature of DFC is then usually captured in measures such as standard deviation (Kucyi and Davis, 2014; Falahpour et al., 2016) or coefficient of variation (Gonzalez-Castillo et al., 2014), however these descriptive measures may be insensitive to finer temporal structuring, which may well be present in DFC.

Although large-scale DFC attracted increasing attention only recently, the dynamic nature of the functional coupling between neuronal cell assemblies had been addressed earlier (Friston, 2000). In fact, functional connectivity for modalities like electroencephalography (EEG) and magnetoencephalography (MEG) was reported having non-linear characteristics (Stam and van Dijk, 2002; Stam et al., 2003). Also, it was shown that several properties of DFC did not have a characteristic time-scale, instead they showed *scale-free* (fractal) dynamics; Gong et al. (2003) presented that fluctuations in phase synchronization between brain regions were scale-free with the characterizing exponent being stable across multiple subjects. Stam and de Bruin (2004) investigated DFC in terms of global synchronization and found that in the alpha and beta bands it scaled with a higher exponent with eyes closed than open. EEG microstates—periods where EEG topography remains constant for 80–120 ms (Lehmann et al., 1987)—also exhibited fractal dynamics as reported by Van de Ville et al. (2010). While these studies evidenced fine, complex temporal structuring present in functional connectivity dynamics both on global (state) and local (individual connection) levels, to the best of our knowledge still only a few studies investigated the scale-free nature of DFC.

Global scale-free (i.e., monofractal) behavior is most commonly characterized by the Hurst exponent (*H*) in the time-, and by the negative power spectral slope (i.e., scaling exponent, β) in the frequency domain (Eke et al., 2000, 2002). *H* and β are inherently interrelated (Eke et al., 2002) as they both characterize the global long-range correlation (LRC) in a signal. This is established by the Wiener-Khinchin theorem stating that the power spectrum is equivalent to the Fourier-transform of the linear autocorrelation function (Kantz and Schreiber, 2004). Describing dynamics through only *H* or β implicitly assumes that frequency components of the power spectrum are independent/random, and information encoded in the phase angles is not considered. According to the definition by Schreiber and Schmitz (2000) this property holds only for linear dynamics. As mentioned above, functional connections in the brain has been shown to be non-linear, which calls for more in-depth analysis techniques capable of providing a detailed-enough description of their dynamic characteristics. At this end, *multifractal* analysis considers scaling as a local instead of a global property of the signal (Mandelbrot, 1986; Tel, 1988; Theiler, 1990), yielding a set of exponents characterizing the scaling in the signal (Kantelhardt et al., 2002; Barunik and Kristoufek, 2010). Moreover, using the method of decomposing a signal into the sign and magnitude time series of its increments (Ashkenazy et al., 2001), it has been shown that multifractal properties of a signal correlated well with its degree of non-linearity (Ashkenazy et al., 2003; Gomez-Extremera et al., 2016; Bernaola-Galvan et al., 2017). Since the seminal work of Ivanov et al. (1999), a diverse set of physiological processes were shown to exhibit multifractal dynamics such as human heart rate variability (Ivanov et al., 2001; Ashkenazy et al., 2003), motor coordination (Ihlen and Vereijken, 2013) or gait dynamics (Ashkenazy et al., 2002). Multifractal analysis of human heartbeat dynamics was also able to capture the separate effects of sympathetic and parasympathetic blockade (Nunes Amaral et al., 2001) as well as reveal the impact of congestive heart failure as the loss of multifractality, substantiating future clinical and diagnostic applications (Ivanov et al., 1999, 2001, 2004). The multifractal nature of neural dynamics was also reported using several different modalities (Shimizu et al., 2004; Wink et al., 2008; Ihlen and Vereijken, 2010), therefore a multifractal approach appears a proper choice when investigating the supposedly rich dynamic properties of functional connectivity.

A recent study using dynamic graph theoretical analysis of multichannel functional near-infrared spectroscopy (fNIRS) data demonstrated that DFC in the human pre-frontal cortex (PFC) expressed multifractal properties (Racz et al., 2018). In this paper we make an attempt on expanding some of the limitations of this previous study: firstly by investigating DFC based on whole-head measurements instead of relying on those only in the PFC, and secondly by considering not only the temporal evolution of global network properties but describing the dynamic fluctuations of individual connections in the network as well. We estimate dynamic functional connectivity based on whole-head EEG measurements using the synchronization likelihood method (Stam and van Dijk, 2002) and apply dynamic graph theoretical analysis. By doing so, we calculate the temporal evolution of three network measures—Density, Clustering Coefficient, and Efficiency—, in order to characterize separate topological aspects of the dynamic networks. Then, both global DFC (as captured in the fluctuations of these network measures) and individual dynamic connections (captured as the fluctuating synchronization levels between regions) are made subject to multifractal time series analysis to reveal their dynamic properties. We performed EEG measurements in eyes open (EO) and eyes closed (EC) states in male and female subjects, which allowed for exploring differences related to state, gender, and network measure. Furthermore, analysis of individual connections between different brain regions allowed us to show if they express not only mono- but indeed multifractal character and also if they show any particular topological pattern regarding their dynamic properties. Our findings imply that multifractal analysis of the dynamics of global functional connectivity as well as that of individual functional connections may provide a valuable tool when extracting information on the temporal structuring of DFC and carry potentials for experimental and/or clinical applications as well.

## Materials and Methods

### Participants, Experimental Protocol and Data Acquisition

A total of 24 young, healthy volunteers (age: 24.25 ± 2.4 ranging from 20 to 29 years, 12 female) participated in this study. This number of subjects was determined by statistical power analysis of preliminary measurements. The study was approved by the Semmelweis University Regional and Institutional Committee of Science and Research Ethics (ethical approval number: 2017/94) and all subjects provided written informed consent. No participants had reported history of any neurological or psychiatric disorders. Test subjects were instructed not to consume any stimulant (e.g., caffeine) at least 4 h prior to participating and to have at least 6 h of sleep the night before. During the measurement, subjects were seated comfortably in an armchair in a light- and sound-attenuated, electrically sealed room, instructed to remain still, retrain from structured thinking while remaining awake. Resting-state EEG measurements were performed with eyes open (EO) while visually fixating on a dot on a computer screen, followed by another recording session with eyes closed (EC), resulting in four analysis groups: female eyes open, female eyes closed, male eyes open and male eyes closed (F_{EO}, F_{EC}, M_{EO} and M_{EC}, respectively). The recorded signals were visually inspected online, and a session was completed once an artifact-free 305 s record in both EO and EC state have been obtained, which was achieved within 20 min with all subjects.

Measurements were carried out using an Emotiv Epoc+ wireless EEG system (Emotiv Systems Inc., San Francisco, CA, USA), acquiring signals from 14 brain regions according to the 10–20 international system, including AF3, F3, F7 FC5, T7, P7, O1, O2, P8, T8, FC6, F8, F4, and AF4 with additional CMS/DRL reference electrodes at P3 and P4 positions. Data was sampled at 2,048 Hz and internally band-pass filtered between 0.2 and 45 Hz with additional notch filters at 50 and 60 Hz, then down-sampled to an effective temporal resolution of 128 Hz. Electrode impedances were kept under 20 kΩ during signal acquisition. All measurements were performed with maximal contact quality confirmed by the provided Emotiv Xavier TestBench^{TM} software (version 3.1.18).

### EEG Data Pre-processing

EEG data was pre-processed off-line in Matlab 2012a (The Mathworks, Natick, MA, USA) using the freely available EEGLAB toolbox (Delorme and Makeig, 2004) along with custom scripts and functions. Independent Component Analysis (ICA) was performed (Hyvarinen and Oja, 2000) to remove signal components related to eye movement, blinking, muscle contraction and general noise. These components were identified by their power spectra, visual appearance, and spatial power distribution over the cortex and rejected from the data before performing inverse ICA. Subsequently, the ICA-pruned datasets were band-pass filtered in the traditional frequency bands used in EEG-analysis: δ, 0.5–4 Hz; θ, 4–8 Hz; α, 8–13 Hz; α, 13–30 Hz, and γ, 30–45 Hz according to Stam and de Bruin (2004). All further analyses were also performed on broadband (unfiltered), ICA-pruned signals as well. The first and last 2 s of each measurement segment were rejected, resulting in datasets with length of 38528 data points for both states of every subject.

### Synchronization Likelihood

The synchronization likelihood (SL) method (Stam and van Dijk, 2002) was used for pairwise estimation of dynamic functional connectivity. SL identifies non-linear statistical interdependencies between a pair (or in case of global SL, a larger set) of signals. It is by its nature dynamic (i.e., estimates synchronization as a function of time), normalized and seems unaffected by non-stationarities (Stam and van Dijk, 2002). These properties make SL a suitable tool in FC studies using EEG measurements, as EEG signals (and brain activity in general) are often considered non-stationary (Kaplan et al., 2005; Freeman and Quian Quiroga, 2013) and the functional coupling between different neuronal ensembles was confirmed to be non-linear by several studies (Friston, 2000; Stam et al., 2003).

SL measures the general synchronization between discretely sampled processes *x*(*t*) and *y*(*t*), *t* = 1, 2, … *T*. First, the temporal evolution of *x*(*t*) and *y*(*t*) is reconstructed in the state space by temporal embedding (Takens, 1981), where *x*(*t*) and *y*(*t*) is converted into a set of state space vectors *X*(*t*) and *Y*(*t*) as

where *d* is the embedding dimension and *m* is the time lag. Further, let's define the probability for every state space vector *X* (*t*) [and for *Y*(*t*), similarly] that the distance of a randomly selected vector *X* (*t*+*u*) is closer than distance *r*_{x}(*t*) as

where *u* is the temporal distance, |·| is the Euclidean norm, θ is the Heaviside step function, *w*_{1} is the Theiler correction for autocorrelation (Theiler, 1986) and *w*_{2} is a window parameter such as *w*_{1}≪*w*_{2}≪*T*. It should be noted, that *w*_{2} serves as the time window in a SW analysis, and as *u* can be negative as well, the ‘window length’ appears as 2*w*_{2}, with the middle 2*w*_{1} segment discarded to avoid effects of autocorrelation. The distance parameters *r*_{x}(*t*) and *r*_{y}(*t*) are set for every time point *t* that *C*(*r*_{x}(*t*), *X*) = *C*(*r*_{y}(*t*), *Y*) = *p*_{ref}, with *p*_{ref} usually fixed at a value close to 0. Thus, *p*_{ref} basically serves as an internal thresholding variable. Finally, the synchronization likelihood at time point *t* is defined as the conditional probability that *Y*(*t*) and *Y*(*t*+*u*) are closer than *r*_{y}(*t*) given that *X*(*t*) and *X*(*t*+*u*) are closer than *r*_{x}(*t*) and calculates as

*SL*(*t*) is then computed for every time point *t*. Note, that the concept of synchronization likelihood is strongly related to the correlation integral (Grassberger and Procaccia, 1983) and can be considered as a so called “fixed-mass” or “*k*-nearest-neighbor” approach (Theiler, 1990).

In this study, initial parameters of *SL* (*d, m, w*_{1}, *w*_{2}) were set to fit the filter parameters for each frequency band, according to Montez et al. (2006). Parameter settings for each frequency band and those for broadband EEG data are shown in Table 1.

### Dynamic Graph Theoretical Analysis

Synchronization likelihood was computed on the pruned EEG datasets for all pair-wise combination of the channels, yielding a 14 × 14 weighted synchronization matrix for every time point, in which the connection strength is assumed to be proportional to the level of synchronization between brain regions. Each of these matrices capture the actual topology of the underlying network, and calculating different network measures over them yield Network Metrics Time Series (NMTS) that describe the temporal evolution of network topology. Complex networks have several aspects to their topologies such as modularity or small-worldness (Rubinov and Sporns, 2010). The network is required to contain a sufficiently large number of nodes for network descriptors to make sense, i.e., there is no point in calculating for example the node degree distribution on a network with 14 nodes. It has been demonstrated however, that global network measures Density (D), Clustering Coefficient (C) and Efficiency (E) can be used effectively to describe and capture significant topological differences in smaller networks (Racz et al., 2017). We used the weighted formulas to calculate the aforementioned network measures. Weighted Density (often termed also Connectivity Strength) is the fraction of overall connectivity strength present to the maximal possible connection strength in a network (Rubinov and Sporns, 2010) and calculates as

where *n* is the total number of nodes, and *c*_{ij} is the connection strength–in this case, *SL*(*t*) for every *t*–between nodes *i* and *j*. Density is a general measure of “wiring cost” of a network, and is also equal to the average normalized node degree (Rubinov and Sporns, 2010). The Clustering Coefficient of an individual node denotes the fraction of the existing triangles to the maximal possible number of triangles around a node (Rubinov and Sporns, 2010) and was generalized to weighted networks by Onnela et al. (2005). Clustering Coefficient of node *i* calculates as

where *k* is the number of edges connected to node *i* and ${\stackrel{~}{w}}_{ih}$ is the weight of the edge between nodes *i* and *j*. Edge weights are normalized by the largest weight present in the network, therefore on binary networks the formula returns with the original definition, that is also equivalent to the fraction of a node's neighbors that are also neighbors of each other (Watts and Strogatz, 1998). The global Clustering Coefficient is the average taken over individual nodes. C is the most general measure of network segregation and related to “local” information flow in the network (Rubinov and Sporns, 2010). Finally, Efficiency is a global parameter capturing network integration, and is often associated with the speed of information processing in a complex system (Rubinov and Sporns, 2010). It is the average inverse shortest path length between all nodes of a network and computes according to

where ${d}_{ij}^{W}$ is the length of the shortest weighted path between nodes *i* and *j* (Latora and Marchiori, 2001). It is strongly related to the average shortest path length, however often considered as being a superior measure to the latter in describing network integration (Achard and Bullmore, 2007). Efficiency is related to the “global” information flow in the network (Rubinov and Sporns, 2010). For the sake of simplicity, in the following we will drop the superscripts “*w”* and will refer by D, C, and E to their weighted forms, respectively.

Since SL has an internal step of thresholding, to avoid acquiring an overwhelming amount of results we decided not to use any additional threshold values, as in a previous study (Racz et al., 2018) the dynamic properties of DFC did not show any specific relation to the value of threshold. We calculated the time evolution of Density, Clustering Coefficient, and Efficiency [D(*t*), C(*t*), and E(*t*), respectively] for every subject both in EO and EC states. For the calculations of D, C, and E we used functions of the Brain Connectivity Toolbox by Rubinov and Sporns (2010).

### Focus-Based Multifractal Signal Summation Conversion (FMF-SSC)

Multifractal analysis, instead of a single scaling exponent yields a set of scaling exponents, each describing scaling in fluctuations of different sizes in the signal (Kantelhardt et al., 2002). This can be achieved by characterizing scaling at several statistical moments *q*, where small fluctuations are amplified by the negative, while large fluctuations by the positive moments (Barunik and Kristoufek, 2010; Ihlen, 2013). We used the multifractal generalization of the Signal Summation Conversion (SSC) method (Eke et al., 2000; Mukli et al., 2015) to extensively characterize the power-law scaling of the NMTSs and the SL(*t*)s of individual connections. In SSC, the signal is cumulatively summed and standard deviation σ is calculated at different window sizes ranging from a minimal to a maximal scale (*s*_{min} and *s*_{max}, respectively). Within each window the local linear trend is removed before calculation to diminish effects of non-stationarity. The power-law dependence of σ on *s* is captured in the Hurst exponent *H* according to σ(*s*)∝*s*^{H} (Eke et al., 2000). The multifractal generalization of SSC (MF-SSC) consists of repeating the steps of the analysis at different statistical moments *q*, yielding the unified scaling function *S*(*q*,*s*) (Figure 1) (Mukli et al., 2015)

**Figure 1**. The scaling function. *H*(2) characterizes the global i.e., monofractal scaling, while the difference between *H*(−15) and *H*(15) termed Δ*H*15 captures the degree of multifractality. The focus is used as a reference point for linear regression. Scaling function of D*(t)* of subject male 4 during eyes closed is used for illustration.

where *N*_{s} is the number of non-overlapping time windows at scale *s*, and υ is the index of the actual window of calculation. As σ is now dependent on *s* and *q* as well, their relationship is established via the generalized Hurst exponent *H*(*q*) according to σ(*q, s*)∝*s*^{H(q)} and can be acquired by linear regression on the values of *S*(*q*,*s*) (Barunik and Kristoufek, 2010). On empirical signals with finite length the regressed functions converge in one specific point termed *focus* (Mukli et al., 2015). This can be simply shown by replacing *s* in *S*(*q*,*s*) with the total length of the signal as it results in the disappearance of the sum and therefore *q* from equation (7). The focus therefore is used as a reference point during regression, which renders the multifractal analysis of empirical time series very robust and prevents it from ending up with “inverted” or corrupted multifractal spectra (Mukli et al., 2015). Note that this way of handling of empirical signals is not unique to SSC, but can be applied also to other MF methods such as Detrended Fluctuation Analysis or Wavelet Leader methods yielding their focus-based variants (Mukli et al., 2015).

The NMTSs and SL(*t*)s were analyzed by FMF-SSC with the following input parameters as suggested in Mukli et al. (2015): *s*_{min} = 8 with *s* increasing dyadically until *s*_{max} = 8192 and *q* ranging from −15 to +15 with unit increments. The scaling range was defined between scales 2^{3} and N/4, as scaling function values below and over these scales become statistically unreliable (Cannon et al., 1997; Gulich and Zunino, 2014). Scaling windows based on a dyadic scale provide efficient computation. The range of *q* was selected based on Grech and Pamula (2012), Mukli et al. (2015), and Nagy et al. (2017). Global (i.e., monofractal) scaling associated with the long-term memory of the signal was captured in *H*(2), and the degree of multifractality was described by Δ*H*15 = *H*(−15)–*H*(15), as a measure of how much the scaling is different for small and large fluctuations (Figure 1). Note, that Δ*H*15 captures the distribution of local scaling in a signal equivalently with the often used multifractal/singularity spectrum that can be acquired from *S*(*q*,*s*) by Legendre-transformation (Frisch and Parisi, 1985; Chhabra et al., 1989).

### Surrogate Data Testing

Multifractal scaling in a time series can appear as a consequence of different long-range correlations present in the signal, however multifractality can also originate in the heavy-tailed probability distribution of signal values without any correlations whatsoever (Ivanov et al., 1999; Kantelhardt et al., 2002). Also, the finite size and/or simple constant linear autocorrelation of the signal can also produce a so-called “multifractal background noise” (Grech and Pamula, 2012), that has to be distinguished from true multifractality. Thus, a proper surrogate data testing is indispensable when analyzing empirical signals. Therefore, all time series (NMTSs and SL(*t*)s) underwent steps of thorough tests to verify if they indeed showed true multifractality. In all of these steps, *n* = 39 surrogate datasets of equal length were generated with the null hypothesis that the investigated time series cannot be distinguished from its surrogates based on the discriminating statistical measure (Theiler et al., 1992). The null hypothesis was rejected if the discriminating statistic was found outside the mean ± 2σ range calculated from the surrogate datasets, that with *n* = 39 surrogate datasets corresponds to a 95% confidence level (Theiler et al., 1992; Kaplan and Glass, 1995).

First, we tested the presence of *global power-law scaling*. This can best be done in the frequency domain, as a signal with global power-law scaling also has a corresponding $\frac{1}{{f}^{\beta}}$-like spectrum (Eke et al., 2000). Hence, for every time series surrogate datasets with equal β were generated with the spectral synthesis method (Saupe, 1988), and the Kolmogorov-Smirnov distances were calculated to estimate the goodness-of-fit (GoF) of the power-law function. A time series was considered scale-free, if its GoF to a power-law function was within the mean ± 2 σ range of those calculated from surrogate data with known (identical) power-law spectra (Clauset et al., 2009; He, 2011).

Second, we tested the presence of *distribution-related multifractality* by randomly shuffling the values of the time series (Ivanov et al., 1999; Kantelhardt et al., 2002). Shuffling destroys all correlations and reduces the signal to pure random noise but has no effect on the distribution of values. Consequently, if shuffling renders *H*(*q*) ≈ 0.5 (i.e., white noise) for all *q* the case of distribution-related multifractality can be excluded, otherwise presence of scaling (at least partly) could be attributed to a power-law type distribution of signal values.

Finally, we tested if the observed multifractality resulted from the presence of different *long-range correlations*. For this purpose, surrogate datasets were generated by Fourier transforming the signal with Fast Fourier Transformation, randomizing its phases and then performing inverse Fourier transformation (Theiler et al., 1992). Phase randomization leaves the amplitudes and therefore the power spectrum (hence, linear autocorrelation) unaffected (Kantz and Schreiber, 2004), while destroying all non-linear correlations in the signal (Schreiber and Schmitz, 2000). Hence it yields surrogate datasets with equal *H*(2) i.e., *monofractal* characteristic preserved (Eke et al., 2000). As the resulting time series are monofractal signals, this step simultaneously tests for multifractal background noise/true multifractality (Grech and Pamula, 2012) and the presence of non-linearity (Ivanov et al., 1999). Multifractality was considered true and as a sign of non-linear dynamics, if the Δ*H*15 value of a signal was significantly larger than those of surrogate datasets. In the following, we will refer those processes passing all tests as *true multifractals*.

### Analysis Strategy

The flowchart of the analysis steps is shown as a summary of the methods on Figure 2. The first aim of this study was to show if global DFC when investigated by EEG show multifractal properties. For this purpose, the acquired D(*t*), C(*t*), and E(*t*) time series were made subject to surrogate data tests for power-law scaling, long-range autocorrelation and true multifractality as described above, and the fraction of subjects passing each test were computed for all measures. To explore the effect of gender (F and M), network measure (D, C, and E), and state (EO and EC) on the MF properties of global DFC, two-way repeated measures ANOVA tests were performed for each frequency band separately, with gender as grouping variable, while state and network measure as repeated measure factors. Bonferroni *post-hoc* tests were performed to identify significant differences between interactions (in the following marked by ×). Results will be presented in the following manner: for each frequency band, we present the significant differences considering the effect of gender, measure and state separately, along with interaction effects of gender × measure, gender × state and measure × state. When they bear significance, results from the *post-hoc* tests are also discussed in detail.

Our second aim was to show if dynamic functional connections themselves—captured in SL(*t*) between different locations—exhibited multifractal scaling. Therefore, all SL(*t*)s were also processed by the same testing framework as the NMTSs. Fractions of edges in the network passing each test were calculated on the subject level and were averaged combining all four groups (F_{EO}, F_{EC}, M_{EO}, and M_{EC}) to obtain a general picture on the presence of true multifractality in functional connections. Subsequently, to reveal if functional connections express any particular topology in regards of their mono- and multifractal properties, *H*(2) and Δ*H*15 values of the connections were standardized (z-scored) on the subject level and averaged across subjects within all analysis groups collectively. This resulted in group-averaged networks where edge weights represent their corresponding z-scored values [z(*H*(2)) and z(Δ*H*15), respectively]. Also, to investigate the possible correlation between the long-term memory and the degree of multifractality of connections, the Pearson coefficient (*r*) was calculated between the group averaged z(*H*(2)) and z(Δ*H*15) values. In order to distinguish between a true correlation effect and pure coincidence, *r* was calculated on the shuffled data as well when *n* = 39 spatial surrogates (Aaria et al., 2013) were generated by randomly shuffling the z(*H*(2)) and z(Δ*H*15) values [thus destroying possible correlation between *H*(2) and Δ*H*15]. Again, if the *r* of original values was outside the mean ± 2 σ range acquired from surrogate data, we regarded the correlation between *H*(2) and Δ*H*15 as significant.

It should be noted, that multifractality of individual connections were investigated in each group (F_{EO}, F_{EC}, M_{EO}, and M_{EC}) individually as well, however as results were comparable between analysis groups, for the sake of simplicity we decided to present the results of all groups combined. All statistical analyses were carried out in StatSoft Statistica 13.2.

## Results

### Testing for True Multifractality

Results for surrogate data testing of NMTSs (including all groups and states) are shown in Table 2. In the vast majority of the cases, D(*t*), C(*t*), and E(*t*) were proven to have broadband power-law spectra, qualifying them as scale-free (fractal) processes. Shuffling reduced *H*(*q*) approximately to 0.5 in all cases, proving again the presence of long-range correlations in the signals, while also excluding the contribution of distribution-related multifractality. True multifractality was also present in network dynamics in most of the cases for all frequency bands as well for broadband EEG data (Table 2).

Similar results were obtained when investigating the power-law scaling and autocorrelation properties of individual connections in the dynamic functional networks on the subject level (Table 3). In that most connections in the networks yielded power-law spectra, and all connections contained LRCs while distribution-related multifractality could be excluded. Although the fraction of true multifractal connections were generally high in the delta, theta and alpha bands, slightly lower values were found in functional networks of the higher (beta and gamma) frequency bands with a slightly lower amount of connections passing the true multifractality test in broadband EEG networks. Nevertheless, these results compare well with those obtained from global network dynamics (Table 2) implying that the multifractal nature of global network dynamics emerges from multifractally fluctuating individual connection strengths in the network. Note that values in Table 3 refer to the fraction of connections in the functional networks of each subject, thus they are presented as mean ± σ across subjects (combining all four groups).

### Multifractal Nature of Global DFC

In the following, for the *H*(2) of D(*t*), C(*t*), and E(*t*) we will use the abbreviations *H*_{D}(2), *H*_{C}(2), and *H*_{E}(2), respectively. Similarly, Δ*H*15 of D(*t*), C(*t*), and E(*t*) will be referred to as Δ*H*_{D}15, Δ*H*_{C}15, and Δ*H*_{E}15. Summary of the results regarding *H*(2) of global DFC is shown in Figure 3, while those of Δ*H*15 are presented on Figure 4.

**Figure 3**. Summary of results regarding *H*(2) of global network metrics for the delta **(A)**, theta **(B)**, alpha **(C)**, beta **(D)**, and gamma **(E)** bands, as well as for broadband EEG **(F)**. Dotted lines connect the means of the corresponding analysis groups. Significant differences are marked by * for effects of gender, measure and state (displayed on the left, in front and right of each subplot, respectively). Note, that significant differences of interaction effects are not marked.

**Figure 4**. Summary of results regarding Δ*H*15 of global network metrics in for the delta **(A)**, theta **(B)**, alpha **(C)**, beta **(D)**, and gamma **(E)** bands, as well as for broadband EEG **(F)**. Dotted lines connect the means of the corresponding analysis groups. Significant differences are marked by * for effects of gender and measure (displayed on the left and in front of each subgraph, respectively). Again, significant differences of interaction effects are not marked.

A prominent measure-related effect appeared in *H*(2) which was present in all frequency bands, as *H*_{C}(2) values were found significantly lower than those of *H*_{D}(2) and *H*_{E}(2). The same difference was found regarding Δ*H*15 (Δ*H*_{D}15 > Δ*H*_{C}15 < Δ*H*_{E}15), however in the delta band, only. Moreover, state had a significant effect in the alpha and beta bands as *H*(2) values increased during EC compared to EO condition. Gender related effects were found sparsely with a tendency of higher *H*(2) and Δ*H*15 values in males. In the following, we elaborate on these results for every frequency band separately.

#### Delta Band

*H*(2) values were found significantly higher in the male groups (main effect of gender, *p* = 0.023). The three network measures also showed significant difference with *H*_{C}(2) being lower than those of *H*_{D}(2) and *H*_{E}(2) (main effect of measure, *p* < 0.0001). Interaction of measure × gender and Bonferroni *post hoc* tests showed that these differences occurred in both genders (Figure 3A).

Measure had a significant main effect on Δ*H*15 as well (*p* < 0.0001). The measure × gender interaction revealed, that in females Δ*H*_{D}15 was significantly larger than Δ*H*_{C}15 andΔ*H*_{E}15 with no difference between the latter two, while in males Δ*H*_{E}15 was also larger than Δ*H*_{C}15 (*p* < 0.01 in all cases). Neither gender nor state had significant effect on Δ*H*15 in the delta band (Figure 4A).

#### Theta Band

Measure had a similar main effect on *H*(2) as in the delta band [*H*_{D}(2) > *H*_{C}(2) < *H*_{E}(2), *p* < 0.001 in all cases], however neither gender- nor state-related differences were found (Figure 3B).

Regarding Δ*H*15, only measure-related differences were found (main effect of measure, *p* = 0.026). Bonferroni *post hoc* test indicated that Δ*H*_{D}15 was higher than Δ*H*_{C}15, however the interaction measure × gender showed that this was not significant in neither the male or female groups individually (Figure 4B).

#### Alpha Band

The previously observed difference in measure regarding *H*(2) was also found in the alpha band [*H*_{D}(2) > *H*_{C}(2) < *H*_{E}(2), *p* < 0.001 in all cases], however additionally in males *H*_{E}(2) was slightly higher than *H*_{D}(2) and *H*_{C}(2) (*p* = 0.043). Also, a significant difference related to state was revealed with *H*(2) being higher in EC than in EO state (*p* = 0.035). Interaction of measure × state verified that this difference was present in all network measures (Figure 3C).

Larger Δ*H*15 values were found in the male groups (main effect of gender, *p* = 0.014). Interestingly, the interaction gender × state revealed a trend in which state had the opposite effect in the two genders, as Δ*H*15 increased in males while decreased it in females during EC condition (*p* = 0.066) (Figure 4C).

#### Beta Band

The main effect of measure was found significant (*p* < 0.0001) with Bonferroni *post hoc* test revealing *H*_{C}(2) being smaller than *H*_{D}(2) and *H*_{E}(2) (*p* < 0.001 in all cases). Moreover, similarly to the alpha band, state had a significant effect in increasing *H*(2) when transitioning from EO to EC (*p* < 0.001) (Figure 3D).

Δ*H*_{D}15, Δ*H*_{C}15, andΔ*H*_{E}15 were found increased in EC state, although the main effect of state was not significant (*p* = 0.136). The interaction of gender × state revealed that this increase was only present in the male group, while Δ*H*15 values remained unchanged in females, although yet again only in tendency (*p* = 0.093) (Figure 4D).

#### Gamma Band

In *H*(2), only measure-related differences appeared significant, with a similar tendency as in the alpha band, with *H*_{C}(2) being lower than *H*_{D}(2) and *H*_{E}(2), and in males *H*_{E}(2) also being higher than *H*_{D}(2) too (*p* < 0.001 in all cases) (Figure 3E). No differences were found regarding Δ*H*15 (Figure 4E).

#### Broadband EEG

Measure related differences were found (main effect of measure, *p* < 0.0001) with Bonferroni *post hoc* test confirming the same tendency as in most previous cases [*H*_{D}(2) > *H*_{C}(2) < *H*_{E}(2), *p* < 0.01]. Also, higher *H*(2) values were found in EC (main effect of state, *p* = 0.014) (Figure 3F).

The main effect of gender on Δ*H*15 appeared significant (*p* = 0.012) with higher values in the male groups. Measure also had a significant effect on Δ*H*15 (*p* = 0.0006) however the measure × gender interaction and Bonferroni *post hoc* test revealed that significant differences only occurred in the male group, with Δ*H*_{D}15 being significantly higher than Δ*H*_{C}15 and Δ*H*_{E}15 (*p* < 0.001 in both cases) (Figure 4F).

In order to keep the statistical analysis simple, frequency band was not included as a factor in the statistical analysis, however it can be clearly seen on Figure 3 that *H*(2) values in the delta and alpha band are considerably higher than those of the theta, beta, and gamma bands, with *H*(2) of unfiltered signals being approximately in between. Figure 4 shows that Δ*H*15 values were the highest in the delta-, alpha-, and broad-band EEG, with slightly lower values in the theta, beta, and gamma bands.

### Individual Connection Dynamics

#### Long-Range Correlation–*H*(2)

*H*(2) of individual connections showed a characteristic topological distribution in all frequency bands as well as broadband EEG data (Figure 5). In that spatially proximal functional connections between ROIs expressed higher, while those between distant regions showed lower *H*(2) values. This pattern could be observed most prominently in connections between regions of the frontal and pre-frontal cortex that had the highest z(*H*(2)) values, while in general connections linking regions of the occipital and parietal cortex to those of the frontal- and pre-frontal cortex had the lowest z(*H*(2)) values. Also, nearby connections linking contralateral regions in the frontal and pre-frontal areas expressed higher long-range correlation, however this did not hold for connections linking areas of the occipital and parietal cortex to those of the contralateral hemisphere (Figure 5).

**Figure 5**. Topology of global long-range correlations in dynamic functional connections across the brain cortex. In all frequency bands as well as in broadband EEG, connections between frontal and prefrontal regions exhibit higher, while occipito-frontal connections lower *H*(2) values.

#### Degree of Multifractality–Δ*H*15

Topology of Δ*H*15 of individual connections in the delta band were markedly different from those of other frequency bands (Figure 6). While dynamic connections in the theta, alpha, beta and gamma bands resembled a topology similar to that observed for *H*(2) (with higher values for short- and lower-values for long-distance connections), a markedly opposite distribution appeared in the delta band. Within this range, stronger multifractality was found in connections linking the occipital and parietal regions to pre-frontal and frontal regions while lower Δ*H*15 values appeared in connections between nearby regions. Connections estimated on broadband EEG appeared to exhibit a homogenous distribution of Δ*H*15 without any particular topology (Figure 6).

**Figure 6**. Topology of the degree of multifractality in dynamic functional connections across the brain cortex. Connections of the delta band showed a markedly different topology from the rest of the frequency bands. Also, no particular topological pattern was apparent in connections estimated on broadband EEG.

#### Relationship Between z(*H*(2)) and z(Δ*H*15)

The relationship between the multifractal properties of dynamic functional connections was captured in the Pearson correlation coefficient between z(*H*(2)) and z(Δ*H*15) values with their scatter plots shown in Figure 7. In the theta, alpha, beta, and gamma bands, significant positive correlations were found indicating that connections with higher long-term autocorrelation was associated with higher degree of multifractality. The opposite topology of z(*H*(2)) and z(Δ*H*15) in the delta band however resulted in a strong negative correlation as connections with high *H*(2) were associated with lower Δ*H*15 values. In broadband EEG no correlation was found concomitant with the absence of topology in Δ*H*15 of individual connections (Figure 5).

**Figure 7**. Cross-correlation between *z(H*(2)) and *z*(Δ*H*15) of individual connections. Panels are scatter plots of all connections (14*13/2) of all subjects (48). The Pearson correlation coefficient (*r*) is displayed in the lower right corner of each panel. Significant correlation is indicated by red regression line, while the black regression line for broadband EEG indicates that the correlation was indistinguishable from spatial surrogates.

## Discussion

### Multifractal Nature of Global DFC

In this study, we reported that dynamic functional connectivity of the brain—as reconstructed from 14 channel whole-head EEG measurements and captured with dynamic graph theoretical analysis—fluctuates according to multifractal dynamics. Surrogate data tests proved that in majority, this temporal structuring was of true multifractal nature in all frequency bands as well as in broadband EEG. We also identified several significant differences in the MF characteristics of global DFC related to network measure, gender, and state.

We found that the degree of long-range temporal correlation could be attributed to specific topological aspects of the dynamic functional networks in that *H*_{C}(2) was found significantly lower than *H*_{D}(2) and *H*_{E}(2). This pattern was present almost universally in all frequency bands as well as in broadband EEG data. A very similar pattern regarding the *H*(2) of the same dynamic graph theoretical measures (D, C, and E) was reported previously in Racz et al. (2018), where multifractal nature of DFC in the pre-frontal cortex was investigated using fNIRS imaging. Findings of the present study and those of Racz et al. (2018) suggest therefore that this pattern is a genuine feature of FC dynamics as it can be captured in different imaging modalities across a broad range of spatio-temporal scales. In Δ*H*15, significant differences related to network measure appeared prominently only in the delta band, where Δ*H*_{C}15 was significantly lower than Δ*H*_{D}15 and Δ*H*_{E}15. These findings are also consistent with the similar results of Racz et al. (2018), where Δ*H*15 showed comparable differences regarding the applied network measures (Δ*H*_{D}15 > Δ*H*_{C}15 < Δ*H*_{E}15). As in this study, this measure-related pattern was found only in the delta band with basically no significant measure-related differences in higher frequency bands. As Racz et al. (2018) investigated spontaneous brain activity in the 0.01–0.1 Hz range, these results may apply for lower-frequency brain activity, only.

Multifractal time series analysis characterizes the time-dependent scaling in a temporal process, that can emerge from intermittent periods with high variance (Ihlen and Vereijken, 2010). In this study, employing three graph theoretical measures–each capturing different aspects of network structure and topology–all were shown fluctuating in a multifractal manner, although to a different extent. The true multifractal scaling of all three network topological measures indicates that their scaling is in fact a local, time-dependent property. This implies that the temporal evolution of dynamic resting-state connection networks was interspersed with short periods of high variance in their segregation and integration, suggesting that these reorganization events leave their impact on localized and global information transfer alike (as captured in C and E, respectively). Moreover, the differences observed in the *H*(2) and Δ*H*15 values of these measures—at least in the delta band—imply that local and global information processing is affected differently: the lower Δ*H*_{C}15 suggesting a more “balanced” temporal structuring of localized activity, while global network integration (related to fast information transfer between distant network sites) and overall “wiring cost” is associated with larger variability in scaling.

The scale-free nature of global DFC has been demonstrated earlier by Stam and de Bruin (2004). In that study, DFC was captured in global synchronization of EEG signals (as acquired by averaging SL time series) and fractal dynamics were characterized by the monofractal Hurst exponent estimated with Detrended Fluctuation Analysis. The authors found an increased Hurst exponent in eyes closed condition in the alpha and beta bands (Stam and de Bruin, 2004). As weighted Density is in fact the averaged connectivity strength (captured in SL), our current observations regarding the same difference–i.e., increased *H*(2) in the alpha and beta bands–are in agreement with those of Stam and de Bruin (2004) thereby extending their validity to Clustering Coefficient and Efficiency. The transition from eyes closed to eyes open state is often associated with desynchronization in cortical activity [Berger (1929), for a review see Barry et al. (2007)], with lower *H*(2) values reflecting on the dominance of rapid fluctuations (Stam and de Bruin, 2004). The lower *H*(2) values found in EO condition for all three network measures may indeed reflect the large-scale network reorganization affecting both local and global information flow (Rubinov and Sporns, 2010; Preti et al., 2017).

In a more recent study, Van de Ville et al. (2010) reported that EEG microstates also exhibit scale-free dynamics. Nevertheless, the authors did not find any difference in the scaling between various statistical moments as EEG microstates were found fluctuating in a strictly monofractal manner. However, this absence of MF dynamics could be the consequence of the applied analysis method. Specifically, it has been shown previously, that brain electrical activity fluctuates between only four different microstates in resting state (Lehmann et al., 1998). Van de Ville et al. (2010) used bipartitioning between these four microstates to capture microstate transitions, and confirmed fractal dynamics in all possible bipartitioning scenarios independently from the applied partition itself (Van de Ville et al., 2010). The random walk time series analyzed by Van de Ville et al. (2010) were acquired by cumulatively summing the bipartition label sequences (consisting of −1 and +1 values) of the EEG microstates. These are reminiscent of the sign time series obtained from an increment series (Ashkenazy et al., 2001, 2003), that can be obtained as the sign sequence of the local derivatives of a fluctuating signal. Sign time series were indeed shown to be related to the monofractal character (Gomez-Extremera et al., 2016), while scaling in the magnitude time series acquired as the absolute value sequence of local derivatives was reported to be related to multifractality of a dynamic process (Gomez-Extremera et al., 2016; Bernaola-Galvan et al., 2017). Our results show that dynamic graph theoretical analysis can capture the dynamics of FC in a more detailed fashion than EEG microstates, allowing for the unfolding of finer temporal structuring such as multifractal scaling.

Gender-related differences in DFC has also been reported previously (Yaesoubi et al., 2015a,b). Yaesoubi et al. (2015b) investigated the simultaneous occupation of different FC states (i.e., state combinations termed as “combo states”) and showed that males occupy a larger fraction of all possible combo states than females. Since spatio-temporally overlapping events were suggested as a possible origin of multifractal dynamics in a complex system (Lima et al., 2017), this may indeed be the underlying reason of the higher degree of multifractality of DFC observed in males in the alpha band. We also found higher autocorrelation i.e., higher *H*(2) values in male subjects in the delta band. It is well-known, that functional connectivity correlates well with structural (anatomical) connectivity (SC) (Greicius et al., 2009; Honey et al., 2009) even when investigated on multiple time scales (Honey et al., 2007). Our present findings showed that in some cases gender can have an influence on FC dynamics, which could well be–at least in part–to gender-related differences in brain anatomical connectivity (Ingalhalikar et al., 2014).

Finally, 1/*f* ^{β} i.e., scale-free neurodynamics were suggested in numerous studies as an indication of an underlying self-organized critical state of the brain (Linkenkaer-Hansen et al., 2001; Stam and de Bruin, 2004; Stam, 2005; Kitzbichler et al., 2009; Chialvo, 2010; Van de Ville et al., 2010; Tagliazucchi et al., 2012a; Racz et al., 2018). Moreover, it has been reported that not only mono- but indeed multifractal scaling could also occur in a three dimensional system at a critical state (Lima et al., 2017). The concept of self-organized criticality (SOC) as introduced by Bak et al. (1987) refers to the intrinsic property of a dynamic system constantly approaching a critical state without the fine tuning of an external “control parameter.” A critical state—where local perturbations are allowed to dissipate across all spatial and temporal scales in the system—usually appears near (first- or secondorder) phase transitions. Therefore, SOC is often considered as an ideal state of the brain in which fast adaptation to changes in the external or internal environment can be effectively achieved by rapid large-scale reorganization (Van de Ville et al., 2010). The scale-free–and indeed, multifractal–nature of the DFC measures reported in this study may therefore be considered evidencing yet another implication of the self-organized critical nature of resting-state brain activity.

### Multifractal Dynamics of Individual Functional Connections

We showed that DFC networks exhibited multifractal dynamics not only in their global graph theoretical parameters, but in their dynamic functional connections, too, as captured in the varying connectivity strengths between the nodes. This property was most prevalent in the delta, theta and alpha band connections, while slightly lower fraction of beta and gamma band as well as broadband EEG connections was proven as true multifractals.

Given the moderate spatial resolution in this study, instead of focusing on individual differences between inter-regional connections, we rather investigated if dynamic functional connections showed any global topology in their *H*(2) and Δ*H*15 values and if there was a relationship between them. In all frequency bands as well as in broadband EEG data we found a characteristic topology in the long-term autocorrelation of FC strength, as short-range connections between the frontal and pre-frontal areas tend to have higher *H*(2) values than long-range connections linking the same regions to the occipital and parietal cortex. A very similar topology regarding the degree of multifractality captured in Δ*H*15 was found in the theta, alpha, beta, and gamma bands. On the contrary, the topological pattern was the opposite in the delta band: long-range connections linking occipital regions with mostly the frontal and pre-frontal regions showed higher Δ*H*15 values than those linking nearby regions prominently again the frontal and pre-frontal areas. This inverse relationship in the delta band could be captured in a strong negative correlation between *H*(2) and Δ*H*15 values. However, in the higher frequency bands the topology of Δ*H*15 values were similar to those of *H*(2) thus resulting in strong positive correlation between the two. No correlation was found between the two MF measures in broadband EEG connections. It should be emphasized that a relationship between *H*(2) and Δ*H*15 is indeed non-trivial, as these measures capture two separate properties of the multifractal spectrum (i.e., spectral center and spectral width) (Theiler, 1990; Kantelhardt et al., 2002; Kantelhardt, 2009; Mukli et al., 2015).

It has been shown in dynamic processes, that monofractal scaling [captured in *H*(2) of the scaling function or the center of the multifractal spectrum] is attributed to linear properties, while the degree of multifractality [as assessed equivalently either by Δ*H*15 or the multifractal spectrum width] correlates well with the degree of non-linearity (Gomez-Extremera et al., 2016; Bernaola-Galvan et al., 2017). In addition non-linearity in a power-law correlated signal can also be well estimated by the long-term autocorrelation of its magnitude time series (Ashkenazy et al., 2003; Schmitt et al., 2009). Multifractality often emerges from the presence of intermittent periods of higher variance (Ihlen and Vereijken, 2010). Accordingly, from the topological pattern observed in broadband EEG and in the lower frequency range (i.e., delta band), it is apparent that connections between nearby regions express a lower degree of non-linearity and they are more linearly autocorrelated (however in most cases still non-linear as well), while on the contrary, the opposite is true for long-distance connections where linear autocorrelation is weaker and the dynamics appear more non-linear. Therefore our results are in good agreement with—and presumably reflect the same phenomenon as—previous findings demonstrating the non-linear nature of neuronal synchronization (Stam et al., 2003) and intermittent periods of high synchronization between neuronal cell assemblies (Friston, 2000). Also, this is in line with the fact that brain activity is generally scale-free on the large-scale level while becoming more synchronized when observed smaller spatial scales (Buzsaki, 2006).

Dynamic synchronization levels mainly capture the joint activity of neuron populations of the regions of interest. The fractal and multifractal nature of individual connections therefore may also indicate that a possible critical state of brain activity may not only be present on a global level, but also on smaller spatial scales, too. In a SOC simulation using a modified version of the classic sand pile model, Mukli et al. (2018) investigated the effect of the size and connection density of the cellular automata on the multifractal properties of their dynamics, which was captured in the total number of sand grains resident in the system at a given moment. Their results demonstrated, that with (independently) increasing system size or connection density, both *H*(2) and Δ*H*15 of the system dynamics increased (Mukli et al., 2018). This phenomenon is very reminiscent of the positive correlation found in this study between *H*(2) and Δ*H*15 of individual connections in the higher frequency bands. Therefore, higher *H*(2) and Δ*H*15 values of a dynamic connection may well reflect that the joint activity between the two regions involved larger and/or more densely interconnected neuronal populations. Also, the frequency range of network oscillations is constrained by the networks size because most neuronal connections are local and thus yield high-frequency oscillations (Braitenberg and Schüz, 2013). Accordingly, slower oscillations can only be produced when larger neuronal cell assemblies are involved (Buzsaki and Draguhn, 2004). Higher *H*(2) values indicate the dominance of slow fluctuations thus suggesting the involvement of large neural populations. These considerations are well in line with the results of this study regarding that the highest *H*(2) and Δ*H*15 values were found in connections of the frontal and prefrontal cortex; regions that are both anatomically and functionally well connected and form a high-level association cortex with diverse functionality (Kandel, 2013). On the contrary, lower values were found between both anatomically and functionally distant regions.

Nevertheless the inverse relationship between *H*(2) and Δ*H*15 found in the delta band raises some questions. Delta band activity is not prominent in awake state during physiological conditions, however several studies demonstrated the critical nature of brain dynamics during sleep stages (Lo et al., 2002, 2004, 2013). It also has been shown recently (using an alternate DFC approach termed time delay stability, TDS), that different sleep stages–that are known to be dominated by specific EEG frequency bands, as α oscillations are most prominent during quiet wake and rapid eye movement, while θ and δ fluctuations usually characterize light- and deep sleep, respectively–exhibit FC topologies characteristic to sleep stage and frequency band as well (Bartsch et al., 2015; Liu et al., 2015). Therefore, sleep would probably be a better setting for further investigating the multifractal and critical nature of delta band connections—and their relation to other frequency bands—that is evidently beyond the scope of this present study.

### Limitations, Future Perspectives

Multifractal analysis could reveal relevant information at the global (network) and the local (connection) levels, that otherwise may remain hidden in static and also in most dynamic functional connectivity approaches. MF analysis of DFC carries potential for future applications both in basic science and clinical fields, however one has to consider some methodological difficulties with this approach. For a reliable numerical estimation of fractal and/or multifractal parameters, a signal length of at least a few thousand data points is desired, as well as long-enough measurement time and high temporal resolution so that the signal could represent a sufficiently broad range of temporal scales (Eke et al., 2002; Ihlen, 2013). While these prerequisites can be readily met with imaging modalities such as EEG, fNIRS, or MEG, still major drawbacks remain owing to their lower spatial resolution, their lack of exact source localization and the fact that these techniques cannot access subcortical regions. These limitations can be partly overcome by using fMRI—with spatial normalization even allowing for exact comparison between separate studies—however at the expense of lower temporal resolution and limited signal length, both considerably affecting the applicability of multifractal analysis (Eke et al., 2012). As these methods–e.g., EEG and fMRI–can provide complementary information, the importance of simultaneous EEG-fMRI (and/or fNIRS-fMRI) measurements is indeed crucial in revealing the relationship between EEG-DFC and fMRI-DFC so that DFC could be investigated with high temporal resolution and exact spatial localization, alike.

Although in this study the whole brain cortex is sampled, the spatial resolution is still fairly limited. Using a higher spatial resolution method would not only benefit from a more detailed sampling of brain activity, but would also allow for calculating more complex network measures—i.e., those related with modularity, centrality or network motifs being more complex than triangles—that could reveal even more details on functional brain organization. Also, from Table 1 it is apparent, that gamma band activity may not be well represented in the state space during SL calculation. This is a limitation brought about by the sampling frequency that calls for greater caution when results regarding the gamma band are evaluated.

DFC analyses carry great potentials not only in basic research leading to a better understanding of brain functions, but also in the clinical field, as several studies already demonstrated their applicability in neuropsychiatric diseases such as schizophrenia (SZ) (Sakoglu et al., 2010; Calhoun et al., 2014; Damaraju et al., 2014) or autism (Price et al., 2014). Thus dynamic graph theoretical analyses similar to the one presented in this study could prove a useful and potentially powerful tool when investigating DFC in clinical settings. When analyzing the DFC of the default mode network with fMRI imaging in SZ patients, the same dynamic graph theoretical measures as used in this study were found to be fluctuating around a lower average value than in healthy controls (Du et al., 2016). Moreover in a whole brain DFC study on SZ patients, in addition to the same results it was also demonstrated that D(*t*), C(*t*), and E(*t*) showed less variance than in healthy controls (Yu et al., 2015). Results of the present study clearly show, that the standard deviation (i.e., the square root of variance) of graph theory metrics depends on the observation scale and that second order statistics alone are insufficient to fully characterize the network dynamics. Multifractal analysis however is suitable to capture such features, therefore may serve as a more sensitive tool in distinguishing physiological states from pathological conditions based on their dynamic graph theoretical measures.

Finally, in a recently established field of biological systems science termed network physiology (Bashan et al., 2012; Ivanov and Bartsch, 2014), the dynamic interactions of local neural activity with several other physiological subsystems (e.g., the cardiac and respiratory system) were analyzed (Bartsch and Ivanov, 2014; Bartsch et al., 2015; Lin et al., 2016). These studies showed that during different physiological states such as sleep stages, the interactions between the elements of this physiological network change significantly. As functional connectivity was also shown to alter during different physiological conditions (Horovitz et al., 2009; Liu et al., 2015; Racz et al., 2017), a dynamic graph theoretical approach may contribute to this emerging field by providing a way for capturing coordinated states of neural activity so that its interactions with other functions of the human body could be further analyzed.

## Conclusions

In this study, we showed that dynamic global functional connectivity of the brain—when investigated by EEG mapping and captured in dynamic graph theoretical measures—fluctuates according to a multifractal temporal pattern. Our results revealed that several network topological aspects exhibit different characteristics. Moreover, the dynamic functional connections assembling these networks showed multifractal dynamics themselves. We found a characteristic topology in both mono- and multi-fractal measures with a positive correlation between them in the higher frequency bands, while anticorrelation in the delta band. Our results suggest that multifractality is indeed a fundamental property of both global and local (i.e., individual) DFC with specific global and local attributes to network topology and anatomical localization, respectively. Our findings are in support of a possible self-organized critical nature of resting-state brain activity. We propose that multifractal analysis can provide a more detailed description of global and local connectivity dynamics than most methods applied in the field, and it could serve as a valuable tool for a better characterization of healthy and pathological brain functions.

## Author Contributions

FR carried out the EEG measurements, developed the analysis framework, performed the analysis on measured data, and wrote the manuscript. OS helped carrying out the measurements and did the data pre-processing. PM helped conducting the measurements and writing the manuscript. AE helped developing and writing the manuscript and provided conceptual guidance in the study.

## Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

The handling editor is currently co-organizing a Research Topic with one of the authors, AE and confirms the absence of any other collaboration.

## References

Aaria, A. T., Mansour, A., Huang, L., Baliki, M. N., Cecchi, G. A., Mesulam, M. M., et al. (2013). Linking human brain local activity fluctuations to structural and functional network architectures. *Neuroimage* 73, 144–155. doi: 10.1016/j.neuroimage.2013.01.072

Achard, S., and Bullmore, E. (2007). Efficiency and cost of economical brain functional networks. *PLoS Comput. Biol.* 3:e17. doi: 10.1371/journal.pcbi.0030017

Allen, E. A., Damaraju, E., Plis, S. M., Erhardt, E. B., Eichele, T., and Calhoun, V. D. (2014). Tracking whole-brain connectivity dynamics in the resting state. *Cerebral Cortex* 24, 663–676. doi: 10.1093/cercor/bhs352

Ashkenazy, Y., Hausdorff, J. M., Ivanov, P. C., and Stanley, H. (2002). A stochastic model of human gait dynamics. *Physica A* 316, 662–670. doi: 10.1016/S0378-4371(02)01453-X

Ashkenazy, Y., Havlin, S., Ivanov, P. C., Peng, C. K., Schulte-Frohlinde, V., and Stanley, H. E. (2003). Magnitude and sign scaling in power-law correlated time series. *Physica A* 323, 19–41. doi: 10.1016/S0378-4371(03)00008-6

Ashkenazy, Y., Ivanov, P. C., Havlin, S., Peng, C. K., Goldberger, A. L., and Stanley, H. E. (2001). Magnitude and sign correlations in heartbeat fluctuations. *Phys. Rev. Lett.* 86, 1900–1903. doi: 10.1103/PhysRevLett.86.1900

Bak, P., Tang, C., and Wiesenfeld, K. (1987). Self-organized criticality: an explanation of the 1/f noise. *Phys. Rev. Lett.* 59, 381–384. doi: 10.1103/PhysRevLett.59.381

Barry, R. J., Clarke, A. R., Johnstone, S. J., and Magee, C. (2007). EEG differences between eyes-closed and eyes-open resting conditions. *Clin. Neurophysiol.* 118, 2765–2773. doi: 10.1016/j.clinph.2007.07.028

Bartsch, R. P., and Ivanov, P. C. (2014). Coexisting forms of coupling and phase-transitions in physiological networks. *Nonlin. Dynam. Electron. Syst.* 438, 270–287. doi: 10.1007/978-3-319-08672-9_33

Bartsch, R. P., Liu, K. K. L., Bashan, A., and Ivanov, P. C. (2015). Network physiology: how organ systems dynamically interact. *PLoS ONE* 10:e0142143. doi: 10.1371/journal.pone.0142143

Barunik, J., and Kristoufek, L. (2010). On hurst exponent estimation under heavy-tailed distributions. *Physica A* 389, 3844–3855. doi: 10.1016/j.physa.2010.05.025

Bashan, A., Bartsch, R. P., Kantelhardt, J. W., and Havlin, S. (2012). Network physiology reveals relations between network topology and physiological function. *Nat. Commun.* 3:702. doi: 10.1038/ncomms1705

Berger, H. (1929). Über das elektrenkephalogramm des menschen. *Archiv für Psychiatrie Nervenkrankheiten* 87, 527–570.

Bernaola-Galvan, P. A., Gomez-Extremera, M., Romance, A. R., and Carpena, P. (2017). Correlations in magnitude series to assess nonlinearities: application to multifractal models and heartbeat fluctuations. *Phys. Rev. E* 96:032218. doi: 10.1103/PhysRevE.96.032218

Biswal, B., Yetkin, F. Z., Haughton, V. M., and Hyde, J. (1995). Functional connectivity in the motor cortex of resting human brain using echo-planar MRI. *Magn. Reson. Med.* 34, 537–541. doi: 10.1002/mrm.1910340409

Braitenberg, V., and Schüz, A. (2013). *Cortex: Statistics and Geometry of Neuronal Connectivity.* Berlin: Springer Science & Business Media.

Bullmore, E., and Sporns, O. (2009). Complex brain networks: graph theoretical analysis of structural and functional systems. *Nat. Rev. Neurosci.* 10, 186–198. doi: 10.1038/nrn2575

Buzsaki, G., and Draguhn, A. (2004). Neuronal oscillations in cortical networks. *Science* 304, 1926–1929. doi: 10.1126/science.1099745

Cader, S., Cifelli, A., Abu-Omar, Y., Palace, J., and Matthews, P. M. (2006). Reduced brain functional reserve and altered functional connectivity in patients with multiple sclerosis. *Brain* 129, 527–537. doi: 10.1093/brain/awh670

Calhoun, V. D., Miller, R., Pearlson, G., and Adali, T. (2014). The chronnectome: time-varying connectivity networks as the next frontier in fMRI data discovery. *Neuron* 84, 262–274. doi: 10.1016/j.neuron.2014.10.015

Cannon, M. J., Percival, D. B., Caccia, D. C., and Raymond, G. (1997). Evaluating scaled windowed variance methods for estimating the Hurst coefficient of time series. *Physica A* 241, 606–626. doi: 10.1016/S0378-4371(97)00252-5

Chang, C., and Glover, G. H. (2010). Time-frequency dynamics of resting-state brain connectivity measured with fMRI. *Neuroimage* 50, 81–98. doi: 10.1016/j.neuroimage.2009.12.011

Chhabra, A. B., Meneveau, C., Jensen, R. V., and Sreenivasan, K. R. (1989). Direct determination of the f(alpha) singularity spectrum and its application to fully developed turbulence. *Phys. Rev. A Gen. Phys.* 40, 5284–5294. doi: 10.1103/PhysRevA.40.5284

Chialvo, D. R. (2010). Emergent complex neural dynamics. *Nat. Phys.* 6, 744–750. doi: 10.1038/Nphys1803

Clauset, A., Shalizi, C. R., and Newman, M. E. J. (2009). Power-law distributions in empirical data. *Siam Rev.* 51, 661–703. doi: 10.1137/070710111

Damaraju, E., Allen, E. A., Belger, A., Ford, J. M., McEwen, S., Mathalon, D. H., et al. (2014). Dynamic functional connectivity analysis reveals transient states of dysconnectivity in schizophrenia. *Neuroimage Clin.* 5, 298–308. doi: 10.1016/j.nicl.2014.07.003

Delorme, A., and Makeig, S. (2004). EEGLAB: an open source toolbox for analysis of single-trial EEG dynamics including independent component analysis. *J. Neurosci. Methods* 134, 9–21. doi: 10.1016/j.jneumeth.2003.10.009

Du, Y., Pearlson, G. D., Yu, Q., He, H., Lin, D., Sui, J., et al. (2016). Interaction among subsystems within default mode network diminished in schizophrenia patients: a dynamic connectivity approach. *Schizophr. Res.* 170, 55–65. doi: 10.1016/j.schres.2015.11.021

Eke, A., Hermán, P., Bassingthwaighte, J. B., Raymond, G. M., Percival, D. B., Cannon, M., et al. (2000). Physiological time series: distinguishing fractal noises from motions. *Pflugers Arch.* 439, 403–415. doi: 10.1007/s004249900135

Eke, A., Herman, P., Kocsis, L., and Kozak, L. R. (2002). Fractal characterization of complexity in temporal physiological signals. *Physiol. Meas.* 23, R1–R38. doi: 10.1088/0967-3334/23/1/201

Eke, A., Herman, P., Sanganahalli, B. G., Hyder, F., Mukli, P., and Nagy, Z. (2012). Pitfalls in fractal time series analysis: fmri bold as an exemplary case. *Front. Physiol.* 3:417. doi: 10.3389/fphys.2012.00417

Esposito, F., Bertolino, A., Scarabino, T., Latorre, V., Blasi, G., Popolizio, T., et al. (2006). Independent component model of the default-mode brain function: Assessing the impact of active thinking. *Brain Res. Bull.* 70, 263–269. doi: 10.1016/j.brainresbull.2006.06.012

Falahpour, M., Thompson, W. K., Abbott, A. E., Jahedi, A., Mulvey, M. E., Datko, M., et al. (2016). Underconnected, but not broken? dynamic functional connectivity mri shows underconnectivity in autism is linked to increased intra-individual variability across time. *Brain Connect* 6, 403–414. doi: 10.1089/brain.2015.0389

Fox, M. D., Snyder, A. Z., Vincent, J. L., Corbetta, M., Van Essen, D. C., and Raichle, M. E. (2005). The human brain is intrinsically organized into dynamic, anticorrelated functional networks. *Proc. Natl. Acad. Sci. U.S.A.* 102, 9673–9678. doi: 10.1073/pnas.0504136102

Freeman, W. J., and Quian Quiroga, R. (2013). *Imaging Brain Function With EEG : Advanced Temporal and Spatial Analysis of Electroencephalographic Signals*. New York, NY: Springer.

Frisch, U., and Parisi, G. (1985). *On the Singularity Structure of Fully Developed Turbulence*. Amsterdam: Turbulence and Predictability in Geophysical Fluid Dynamics and Climate Dynamics, Elsevier.

Friston, K. J. (2000). The labile brain. I. Neuronal transients and nonlinear coupling. *Philos. Trans. R Soc. Lond. B Biol. Sci.* 355, 215–236. doi: 10.1098/rstb.2000.0560

Friston, K. J. (2011). Functional and effective connectivity: a review. *Brain Connect.* 1, 13–36. doi: 10.1089/brain.2011.0008

Friston, K. J., Frith, C. D., Liddle, P. F., and Frackowiak, R. (1993). Functional connectivity: the principal-component analysis of large (PET) data sets. *J. Cereb. Blood Flow Metab.* 13, 5–14. doi: 10.1038/jcbfm.1993.4

Gaudes, C. C., Petridou, N., Francis, S. T., and Dryden, I. L. (2013). Paradigm free mapping with sparse regression automatically detects single-trial functional magnetic resonance imaging blood oxygenation level dependent responses. *Hum. Brain Mapp.* 34, 501–518. doi: 10.1002/hbm.21452

Gomez-Extremera, M., Carpena, P., Ivanov, P., and Bernaola-Galvan, P. A. (2016). Magnitude and sign of long-range correlated time series: decomposition and surrogate signal generation. *Phys. Rev. E* 93:042201. doi: 10.1103/PhysRevE.93.042201

Gong, P., Nikolaev, A. R., and van Leeuwen, C. (2003). Scale-invariant fluctuations of the dynamical synchronization in human brain electrical activity. *Neurosci. Lett.* 336, 33–36. doi: 10.1016/S0304-3940(02)01247-8

Gonzalez-Castillo, J., Handwerker, D. A., Robinson, M. E., Hoy, C. W., Buchanan, L. C., Saad, Z. S., et al. (2014). The spatial structure of resting state connectivity stability on the scale of minutes. *Front Neurosci* 8:138. doi: 10.3389/fnins.2014.00138

Grassberger, P., and Procaccia, I. (1983). Characterization of strange attractors. *Phys. Rev. Lett.* 50, 346–349. doi: 10.1103/PhysRevLett.50.346

Grech, D., and Pamula, G. (2012). Multifractal Background Noise of Monofractal Signals. *Acta Phys Pol A* 121, B34–B39.

Greicius, M. D., Krasnow, B., Reiss, A. L., and Menon, V. (2003). Functional connectivity in the resting brain: a network analysis of the default mode hypothesis. *Proc. Natl. Acad. Sci. U.S.A.* 100, 253–258. doi: 10.1073/pnas.0135058100

Greicius, M. D., Supekar, K., Menon, V., and Dougherty, R. F. (2009). Resting-state functional connectivity reflects structural connectivity in the default mode network. *Cereb. Cortex* 19, 72–78. doi: 10.1093/cercor/bhn059

Gulich, D., and Zunino, L. (2014). A criterion for the determination of optimal scaling ranges in DFA and MF-DFA. *Phys. Stat. Mech. Appl.* 397, 17–30. doi: 10.1016/j.physa.2013.11.029

He, B. J. (2011). Scale-free properties of the functional magnetic resonance imaging signal during rest and task. *J. Neurosci.* 31, 13786–13795. doi: 10.1523/JNEUROSCI.2111-11.2011

Honey, C. J., Kotter, R., Breakspear, M., and Sporns, O. (2007). Network structure of cerebral cortex shapes functional connectivity on multiple time scales. *Proc. Natl. Acad. Sci. U.S.A.* 104, 10240–10245. doi: 10.1073/pnas.0701519104

Honey, C. J., Sporns, O., Cammoun, L., Gigandet, X., Thiran, J. P., Meuli, R., et al. (2009). Predicting human resting-state functional connectivity from structural connectivity. *Proc. Natl. Acad. Sci. U.S.A.* 106, 2035–2040. doi: 10.1073/pnas.0811168106

Horovitz, S. G., Braun, A. R., Carr, W. S., Picchioni, D., Balkin, T. J., Fukunaga, M., et al. (2009). Decoupling of the brain's default mode network during deep sleep. *Proc. Natl. Acad. Sci. U.S.A.* 106, 11376–11381. doi: 10.1073/pnas.0901435106

Hutchison, R. M., Womelsdorf, T., Allen, E. A., Bandettini, P. A., Calhoun, V. D., Corbetta, M., et al. (2013a). Dynamic functional connectivity: promise, issues, and interpretations. *Neuroimage* 80, 360–378. doi: 10.1016/j.neuroimage.2013.05.079

Hutchison, R. M., Womelsdorf, T., Gati, J. S., Everling, S., Menon, R., et al. (2013b). Resting-state networks show dynamic functional connectivity in awake humans and anesthetized macaques. *Hum. Brain Mapp.* 34, 2154–2177. doi: 10.1002/hbm.22058

Hyvarinen, A., and Oja, E. (2000). Independent component analysis: algorithms and applications. *Neural Netw.* 13, 411–430. doi: 10.1016/S0893-6080(00)00026-5

Ihlen, E. A., and Vereijken, B. (2010). Interaction-dominant dynamics in human cognition: beyond 1/f(alpha) fluctuation. *J. Exp. Psychol. Gen.* 139, 436–463. doi: 10.1037/a0019098

Ihlen, E. A., and Vereijken, B. (2013). Identifying multiplicative interactions between temporal scales of human movement variability. *Ann. Biomed. Eng.* 41, 1635–1645. doi: 10.1007/s10439-012-0724-z

Ihlen, E. A. F. (2013). Multifractal analyses of response time series: a comparative study. *Behav. Res. Methods* 45, 928–945. doi: 10.3758/s13428-013-0317-2

Ingalhalikar, M., Smith, A., Parker, D., Satterthwaite, T. D., Elliott, M. A., Ruparel, K., et al. (2014). Sex differences in the structural connectome of the human brain. *Proc. Natl. Acad. Sci. U.S.A.* 111, 823–828. doi: 10.1073/pnas.1316909110

Ivanov, P. C., Amaral, L. A., Goldberger, A. L., Havlin, S., Rosenblum, M. G., Struzik, Z. R., et al. (1999). Multifractality in human heartbeat dynamics. *Nature* 399, 461–465. doi: 10.1038/20924

Ivanov, P. C., and Bartsch, R. P. (2014). Network Physiology: Mapping Interactions Between Networks of Physiologic Networks. *Networks of Networks: the Last Frontier of Complexity.* Berlin: Springer, 203–222.

Ivanov, P. C., Chen, Z., Hu, K., and Stanley, H. E. (2004). Multiscale aspects of cardiac control. *Phys. Stat. Mech. Appl.* 344, 685–704. doi: 10.1016/j.physa.2004.08.016

Ivanov, P. C., Nunes Amaral, L. A., Goldberger, A. L., Rosenblum, M. G., Stanley, H. G., and Struzik, Z. R. (2001). From 1/f noise to multifractal cascades in heartbeat dynamics. *Chaos* 11, 641–652. doi: 10.1063/1.1395631

Kantelhardt, J. W. (2009). *Fractal and Multifractal Time Series. Encyclopedia of Complexity and Systems Science.* Berlin: Springer, 3754–3779.

Kantelhardt, J. W., Zschiegner, S. A., Koscielny-Bunde, E., Bunde, A., Havlin, H., Stanley, H. E., et al. (2002). Multifractal detrended fluctuation analysis of nonstationary time series. *Phys. A* 316, 87–114. doi: 10.1016/S0378-4371(02)01383-3

Kantz, H., and Schreiber, T. (2004). *Nonlinear Time Series Analysis*. Cambridge, NY: Cambridge University Press.

Kaplan, A. Y., Fingelkurts, A. A., Fingelkurts, A. A., and Borisov, S. (2005). Nonstationary nature of the brain activity as revealed by EEG/MEG: Methodological, practical and conceptual challenges. *Signal Process* 85, 2190–2212. doi: 10.1016/j.sigpro.2005.07.010

Kitzbichler, M. G., Smith, M. L., Christensen, S. R., and Bullmore, E. (2009). Broadband criticality of human brain network synchronization. *Plos Comput. Biol.* 5:e1000314. doi: 10.1371/journal.pcbi.1000314

Kucyi, A., and Davis, K. D. (2014). Dynamic functional connectivity of the default mode network tracks daydreaming. *Neuroimage* 100, 471–480. doi: 10.1016/j.neuroimage.2014.06.044

Latora, V., and Marchiori, M. (2001). Efficient behavior of small-world networks. *Phys. Rev. Lett.* 87:198701. doi: 10.1103/PhysRevLett.87.198701

Lehmann, D., Ozaki, H., and Pal, I. (1987). Eeg alpha-map series - brain micro-states by space-oriented adaptive segmentation. *Electroen. Clin. Neuro* 67, 271–288. doi: 10.1016/0013-4694(87)90025-3

Lehmann, D., Strik, W. K., Henggeler, B., Koenig, T., and Koukkou, M. (1998). Brain electric microstates and momentary conscious mind states as building blocks of spontaneous thinking: i. Visual imagery and abstract thoughts. *Int. J. Psychophysiol.* 29, 1–11.

Lima, G. Z. D. S., Corso, G., Correa, M. A., Sommer, R. L., Ivanov, P. C., and Bohn, F. (2017). Universal temporal characteristics and vanishing of multifractality in Barkhausen avalanches. *Phys. Rev. E* 96:022159. doi: 10.1103/PhysRevE.96.022159

Lin, A. J., Liu, K. K. L., Bartsch, R. P., and Ivanov, P (2016). Delay-correlation landscape reveals characteristic time delays of brain rhythms and heart interactions. *Philos. Trans. R. Soc. Math. Phys. Eng. Sci.* 374:20150182. doi: 10.1098/rsta.2015.0182

Linkenkaer-Hansen, K., Nikouline, V. V., Palva, J. M., and Ilmoniemi, R. (2001). Long-range temporal correlations and scaling behavior in human brain oscillations. *J. Neurosci.* 21, 1370–1377. doi: 10.1523/JNEUROSCI.21-04-01370.2001

Liu, K. K., Bartsch, R. P., Lin, A., and Mantegna, R. N (2015). Plasticity of brain wave network interactions and evolution across physiologic states. *Front. Neural Circuits* 9:62. doi: 10.3389/fncir.2015.00062

Liu, Y., Liang, M., Zhou, Y., He, Y., Hao, Y., Song, M., et al. (2008). Disrupted small-world networks in schizophrenia. *Brain* 131, 945–961. doi: 10.1093/brain/awn018

Lo, C.-C., Bartsch, R. P., and Ivanov, P. C. (2013). Asymmetry and basic pathways in sleep-stage transitions. *EPL* 102:10008. doi: 10.1209/0295-5075/102/10008

Lo, C.-C., Nunes Amaral, L. A., Havlin, S., Ivanov, P. C., Penzel, T., Peter, J.-H., et al. (2002). Dynamics of sleep-wake transitions during sleep. *EPL* 57:625. doi: 10.1209/epl/i2002-00508-7

Lo, C. C., Chou, T., Penzel, T., Scammell, T. E., Strecker, R. E., Stanley, H. E., et al. (2004). Common scale-invariant patterns of sleep–wake transitions across mammalian species. *Proc. Natl. Acad. Sci. U.S.A.* 101, 17545–17548. doi: 10.1073/pnas.0408242101

Montez, T., Linkenkaer-Hansen, K., van Dijk, B. W., and Stam, C. J. (2006). Synchronization likelihood with explicit time-frequency priors. *Neuroimage* 33, 1117–1125. doi: 10.1016/j.neuroimage.2006.06.066

Mukli, P., Nagy, Z., and Eke, A. (2015). Multifractal formalism by enforcing the universal behavior of scaling functions. *Phys. A* 417, 150–167. doi: 10.1016/j.physa.2014.09.002

Mukli, P., Nagy, Z., Racz, F. S., Herman, P., and Eke, A. (2018). Impact of healthy aging on multifractal hemodynamic fluctuations in the human prefrontal cortex. *Front Physiol.* 9:1072. doi: 10.3389/fphys.2018.01072

Nagy, Z., Mukli, P., Herman, P., and Eke, A. (2017). Decomposing multifractal crossovers. *Front. Physiol.* 8:533. doi: 10.3389/fphys.2017.00533

Nunes Amaral, L. A., Ivanov, P. C., Aoyagi, N., Hidaka, I., Tomono, S., Goldberger, A. L., et al. (2001). Behavioral-independent features of complex heartbeat dynamics. *Phys. Rev. Lett.* 86, 6026–6029. doi: 10.1103/PhysRevLett.86.6026

Onnela, J. P., Saramaki, J., Kertesz, J., and Kaski, K. (2005). Intensity and coherence of motifs in weighted complex networks. *Phys. Rev. E Stat. Nonlin. Soft Matter Phys.* 71:065103. doi: 10.1103/PhysRevE.71.065103

Pievani, M., de Haan, W., Wu, T., Seeley, W. W., and Frisoni, G. (2011). Functional network disruption in the degenerative dementias. *Lancet Neurol.* 10, 829–843. doi: 10.1016/S1474-4422(11)70158-2

Preti, M. G., Bolton, T. A., and Van De Ville, D. (2017). The dynamic functional connectome: state-of-the-art and perspectives. *Neuroimage* 160, 41–54. doi: 10.1016/j.neuroimage.2016.12.061

Price, T., Wee, C. Y., Gao, W., and Shen, D. G. (2014). Multiple-network classification of childhood autism using functional connectivity dynamics. *Lect. Notes Comput. Sci.* 8675, 177–184. doi: 10.1007/978-3-319-10443-0_23

Rack-Gomer, A. L., and Liu, T. T. (2012). Caffeine increases the temporal variability of resting-state BOLD connectivity in the motor cortex. *Neuroimage* 59, 2994–3002. doi: 10.1016/j.neuroimage.2011.10.001

Racz, F. S., Mukli, P., Nagy, Z., and Eke, A. (2017). Increased prefrontal cortex connectivity during cognitive challenge assessed by fNIRS imaging. *Biomed. Opt. Express* 8, 3842–3855. doi: 10.1364/BOE.8.003842

Racz, F. S., Mukli, P., Nagy, Z., and Eke, A. (2018). Multifractal dynamics of resting-state functional connectivity in the prefrontal cortex. *Physiol. Meas.* 39:024003. doi: 10.1088/1361-6579/aaa916

Raichle, M. E., MacLeod, A. M., Snyde, A. Z., Powers, W. J., Gusnard, D. A., Shulman, G. L., et al. (2001). A default mode of brain function. *Proc. Natl. Acad. Sci. U.S.A.* 98, 676–682. doi: 10.1073/pnas.98.2.676

Rubinov, M., and Sporns, O. (2010). Complex network measures of brain connectivity: uses and interpretations. *Neuroimage* 52, 1059–1069. doi: 10.1016/j.neuroimage.2009.10.003

Sakoglu, U., Pearlson, G. D., Kiehl, K. A., Wang, Y. M., Michael, A. M., and Calhoun, V. D. (2010). A method for evaluating dynamic functional network connectivity and task-modulation: application to schizophrenia. *Magn. Reson. Mater. Phy.* 23, 351–366. doi: 10.1007/s10334-010-0197-8

Saupe, D. (1988). “Algorithms for random fractals,” in *The Science of Fractal Images*, eds H. O. Peitgen and D. Saupe (Springer), 71–136.

Schmitt, D. T., Stein, P. K., and Ivanov, P. C. (2009). Stratification pattern of static and scale-invariant dynamic measures of heartbeat fluctuations across sleep stages in young and elderly. *Ieee Trans. Biomed. Eng.* 56, 1564–1573. doi: 10.1109/Tbme.2009.2014819

Schreiber, T., and Schmitz, A. (2000). Surrogate time series. *Phys. D* 142, 346–382. doi: 10.1016/S0167-2789(00)00043-9

Shimizu, Y., Barth, M., Windischberger, C., Moser, E., and Thurner, S. (2004). Wavelet-based multifractal analysis of fMRI time series. *Neuroimage* 22, 1195–1202. doi: 10.1016/j.neuroimage.2004.03.007

Sporns, O. (2011). The human connectome: a complex network. *Ann. N. Y. Acad. Sci.* 1224, 109–125. doi: 10.1111/j.1749-6632.2010.05888.x

Sporns, O., Tononi, G., and Kotter, R. (2005). The human connectome: a structural description of the human brain. *PLoS Comput. Biol.* 1:e42. doi: 10.1371/journal.pcbi.0010042

Stam, C. J. (2005). Nonlinear dynamical analysis of EEG and MEG: review of an emerging field. *Clin. Neurophysiol.* 116, 2266–2301. doi: 10.1016/j.clinph.2005.06.011

Stam, C. J., Breakspear, M. A. M., van Cappellen, W., and van Dijk, B. W. (2003). Nonlinear synchronization in EEG and whole-head MEG recordings of healthy subjects. *Hum. Brain Mapp.* 19, 63–78. doi: 10.1002/hbm.10106

Stam, C. J., and de Bruin, E. A. (2004). Scale-free dynamics of global functional connectivity in the human brain. *Hum. Brain Mapp.* 22, 97–109. doi: 10.1002/hbm.20016

Stam, C. J., and van Dijk, B. W. (2002). Synchronization likelihood: an unbiased measure of generalized synchronization in multivariate data sets. *Phys. D* 163, 236–251. doi: 10.1016/S0167-2789(01)00386-4

Tagliazucchi, E., Balenzuela, P., Fraiman, D., and Chialvo, D. R. (2012a). Criticality in large-scale brain FMRI dynamics unveiled by a novel point process analysis. *Front. Physiol.* 3:15. doi: 10.3389/fphys.2012.00015

Tagliazucchi, E., von Wegner, F., Morzelewski, A., Brodbeck, V., and Laufs, H. (2012b). Dynamic BOLD functional connectivity in humans and its electrophysiological correlates. *Front. Hum. Neurosci.* 6:339. doi: 10.3389/fnhum.2012.00339

Takens, F. (1981). “Detecting strange attractors in turbulence,” in *Dynamical Systems and Turbulence, Warwick 1980*, eds D. Rand and L.S. Young (Berlin: Springer), 366–381.

Tel, T. (1988). Fractals, multifractals, and thermodynamics - an introductory review. *Z Naturforsch A* 43, 1154–1174. doi: 10.1515/zna-1988-1221

Theiler, J. (1986). spurious dimension from correlation algorithms applied to limited time-series data. *Phys. Rev. A* 34, 2427–2432. doi: 10.1103/PhysRevA.34.2427

Theiler, J. (1990). Estimating fractal dimension. *J. Opt. Soc. Am. A* 7, 1055–1073. doi: 10.1364/Josaa.7.001055

Theiler, J., Eubank, S., Longtin, A., Galdrikian, B., and Farmer, J. D. (1992). Testing for nonlinearity in time-series - the method of surrogate data. *Phys. D* 58, 77–94. doi: 10.1016/0167-2789(92)90102-S

Van de Ville, D., Britz, J., and Michel, C. M. (2010). EEG microstate sequences in healthy humans at rest reveal scale-free dynamics. *Proc. Natl. Acad. Sci. U.S.A.* 107, 18179–18184. doi: 10.1073/pnas.1007841107

van den Heuvel, M. P., and Hulshoff Pol, H. E. (2010). Exploring the brain network: a review on resting-state fMRI functional connectivity. *Eur. Neuropsychopharmacol.* 20, 519–534. doi: 10.1016/j.euroneuro.2010.03.008

Watts, D. J., and Strogatz, S. H. (1998). Collective dynamics of ‘small-world’ networks. *Nature* 393, 440–442. doi: 10.1038/30918

Werner, G. (2010). Fractals in the nervous system: conceptual implications for theoretical neuroscience. *Front. Physiol.* 1:15. doi: 10.3389/Fphys.2010.00015

Wink, A. M., Bullmore, E., Barnes, A., Bernard, F., and Suckling, J. (2008). Monofractal and multifractal dynamics of low frequency endogenous brain oscillations in functional MRI. *Hum. Brain Mapp.* 29, 791–801. doi: 10.1002/hbm.20593

Yaesoubi, M., Allen, E. A., Miller, R. L., and Calhoun, V. (2015a). Dynamic coherence analysis of resting fMRI data to jointly capture state-based phase, frequency, and time-domain information. *Neuroimage* 120, 133–142. doi: 10.1016/j.neuroimage.2015.07.002

Yaesoubi, M., Miller, R. L., and Calhoun, V. D. (2015b). Mutually temporally independent connectivity patterns: A new framework to study the dynamics of brain connectivity at rest with application to explain group difference based on gender. *Neuroimage* 107, 85–94. doi: 10.1016/j.neuroimage.2014.11.054

Keywords: functional connectivity, dynamic functional connectivity, multifractal analyses, brain, synchronization likelihood (SL), self-organized criticality (SOC), electroencephalography (EEG)

Citation: Racz FS, Stylianou O, Mukli P and Eke A (2018) Multifractal Dynamic Functional Connectivity in the Resting-State Brain. *Front. Physiol*. 9:1704. doi: 10.3389/fphys.2018.01704

Received: 19 July 2018; Accepted: 12 November 2018;

Published: 30 November 2018.

Edited by:

Plamen Ch. Ivanov, Boston University, United StatesReviewed by:

Chengyu Huo, Changshu Institute of Technology, ChinaSebastiano Stramaglia, Università degli Studi di Bari, Italy

Copyright © 2018 Racz, Stylianou, Mukli and Eke. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Andras Eke, eke.andras@med.semmelweis-univ.hu