Scaling of whole-brain dynamics reproduced by high-order moments of turbulence indicators

We investigate how brain activity can be supported by a turbulent regime based on the deviations of a self-similar scaling of high-order structure functions within the phenomenological Kolmogorov’s theory. By analyzing a large neuroimaging data set, we establish the relationship between scaling exponents and their order, showing that brain activity has more than one invariant scale, and thus orders higher than 2 are needed to accurately describe its underlying statistical properties. Furthermore, we build whole-brain models of coupled oscillators to show that high-order information allows for a better description of the brain’s empirical information transmission and reactivity.

scaling law (or, equivalently, one critical exponent).Instead, many critical exponents are needed to identify this behavior.
A key outstanding question is, thus, whether scaling laws of higher-order structure functions akin to those found in turbulence in fluid dynamics are also found in neuronal dynamics.If so, this would provide fundamental information about the nature of information transfer in the brain and how it changes in disease or in altered states of consciousness [10,11].In turn, this could help more accurate whole-brain models [12], which could have important translational consequences for providing more precise biomarkers for neuropsychiatric disorders that may eventually lead to better treatments [13].
Kolmogorov's scaling theory for the inertial range in turbulence makes a claim for the power-law scaling of the second-order structure function where x and x + r are points in the fluid, and u is the component of the fluid velocity along r.The range where the powerlaw scaling holds (i.e., the inertial range) is bounded between the dissipation and the integral scales of motion [14,15].Deco and Kringelbach adapted Kolmogorov's second-order structure function (termed the "turbulent core" by the authors) by defining u as the blood-oxygen-level-dependent (BOLD) signal and found an inertial range and its corresponding power law in brain dynamics data [16].Nevertheless, the statistical properties of human brain dynamics in terms of higher-order structure functions have not been investigated.These functions are important as they provide information on intermittency, fractality, and out-of-equilibrium dynamics and are thus required to fully characterize turbulence in any physical system.The extended Kolmogorov theory for higher orders (p > 2) claims that the scaling results generalize to structure functions of any order as with p > 0, i.e., they present a power-law scaling with exponent ζ changing with the order as ζ (p) = p/3 [3].This corresponds to perfect (Gaussian) scale invariance.However, extensive numerical and experimental investigations have found deviations from Kolmogorov's prediction, which are due to intermittency (i.e., localized spatial inhomogeneities) of the energy dissipation rate that cannot be treated as uniform [17].In spite of these deviations, the existence of universal scaling laws for higher-order structure functions has been verified for a variety of different turbulent flows [18][19][20].
Here we study these out-of-equilibrium indicators in brain dynamics, by investigating higher-order structure functions within the inertial range determined from the second-order moment [16].We show that brain signals present more than one invariant scale, and thus high orders are necessary to describe the statistics, i.e., second-order structure functions are not sufficient for describing underlying statistical properties.To this end we compute the power-law scaling ranges for structure functions up to the eighth order, by considering u as the BOLD signal in Eq. ( 2) for the empirical brain resting from the 1003 Human Connectome Project (HCP) participants.The computation of scaling exponents using extended self-similarity (ESS) [19] allows us to confirm the scaling.We then examine the relationship between these exponents and the order of the structure function to investigate the relevant scales present in the system.The nontrivial dependence of the exponents with the order confirms that the analogy with turbulent systems includes important features such as deviations from perfect scale invariance and intermittency.Finally, we ask whether higher-order structure functions provide useful information for whole-brain computational models built upon simple local dynamical rules coupled according to empirical measurements of anatomical connectivity [21], i.e., white matter tracts inferred from diffusion tensor imaging [22].Typically, a parameter optimization procedure is conducted to reproduce an empirical observable [23][24][25], perhaps the most common being the functional connectivity, identified with the second-order structure function.Implementing local dynamics with Stuart-Landau oscillators and considering the second-order structure function as an empirical target, we examine the effect of including high-order structure functions on the model performance and its capacity to reproduce information processing properties found in empirical data such as information cascades and fast reactiveness to external stimuli [10,11,16].

II. EMPIRICAL HIGH-ORDER STRUCTURE FUNCTIONS COMPUTED FROM RESTING-STATE FUNCTIONAL MRI DATA SHOW POWER LAWS
First, we extracted and bandpass filtered (0.008-0.08 Hz) the time series from the 1003 Human Connectome Project (HCP) resting-state functional magnetic resonance imaging (fMRI) recordings, each time series corresponding to the average of a parcel within the Schaefer 1000 parcellation.We adapted the increment of velocity in fluid dynamics, u(r), by computing the difference between BOLD signals at pairs of nodes separated by distance r ± r, where r equals one bin of the spatial grid.To do so, we computed the distances between the geometric centers corresponding to all parcels and then averaged the BOLD signal across all pairs separated by r ± r.For each participant, we calculated the p-order structure functions following Eq.( 2) from p = 1 to p = 8, and then we estimated within the inertial subrange (previously determined by the scaling of the second-order structure function) the power-law scaling exponents corresponding to each order.In Fig. 1(a) we present the structure functions as a function of the distance in their brain r for each of the 1003 participants of the HCP data set (see Appendix A for detailed explanation on the participants, data acquisition, and preprocessing).We observed that the structure functions are consistent across all participants, with odd structure functions showing a mean around 0. Consequently, we considered only the even functions within the inertial range, with r = 8.1 to r = 33.8mm [highlighted in Fig. 1(a)] to compute the scaling exponent ζ (n), as is the standard practice in turbulence studies [18,20].We plotted in a log-log scale the even structure functions, fitted a power law within the inertial range [Fig.1(b)], and found the following spatial power-law scaling exponents: ζ (2) = 0.38 ± 0.06, ζ (4) = 0.59 ± 0.09, ζ (6) = 0.57 ± 0.09, and ζ (8) = 0.44 ± 0.08, with regression coefficients of R 2 = 0.97 and p < 0.001, R 4 = 0.97 and p < 0.001, R 6 = 0.97 and p < 0.001, and R 8 = 0.95 and p < 0.001, respectively.

III. EXTENDED SELF-SIMILARITY VALIDATES THE INERTIAL RANGE FOR HIGHER-ORDER STRUCTURE FUNCTIONS
We leveraged the extended self-similarity approach proposed by Benzi and colleagues to validate the inertial range in turbulent flows, defined by the scaling of the second-order structure function and allowing better estimation of the other scaling exponents [19].Figure 2(a) shows in log-log scale the relationship between high-order structure functions from 3 to 8 and the second-order structure function.We were able to estimate the odd scaling exponents through ESS, including the corresponding structure functions, by considering the absolute value of the increment, a practice also often used in the study of turbulence [18].We investigated the convergence of the accumulated moments of the different orders in our data following the work of Gotoh and colleagues [26].We found that the odd moments were less converged than the even orders and that in both cases the convergence degraded as a function of the increasing order.In particular, we found that the seventh order did not properly converge and that it violated Holder's inequalities.As a consequence, we decided to exclude this order from the analysis.Also, the eighth order did not properly converge although it respected Holder's inequalities.Thus, we considered that the value that we obtained corresponds to a lower bound of the scaling exponent (see Appendix B).After this consideration, and following ESS, we computed the scaling exponents for highorder structure functions by fitting the linear relation between the logarithm of the structure functions and the logarithm of the second-order structure function.We thus obtained the following power-

IV. MULTIFRACTALITY ON WHOLE BRAIN DYNAMICS BASED ON fMRI RECORDINGS OF HEALTHY PARTICIPANTS
We investigated the scale invariance of brain dynamics through the relationship between the scaling exponents and the corresponding order shown in Fig. 2(b).Notably, the exponents could be linearly organized up to order 4 as a signature of one-scale invariance.Nevertheless, for higher orders a change in slope was observed, and at least one different scale could be identified.This could be interpreted as evidence that the system presents at least a bifractal structure, as in the case of fluids following the Burgers equation, where the intermittency is associated with large-scale singularities resulting from a broken statistical symmetry of the disordered system [27].We then investigated the distribution of the first-order structure function, S 1 (r), at different scales within the inertial range (i.e., the probability distribution function of the BOLD spatial fluctuations).We found that, while for small scales the distributions were not Gaussian, as the scale r increased the distributions of S 1 (r) tended monotonically towards Gaussian distributions [Fig.2(c)].Thus, the fluctuations were not perfectly self-similar (in which case all normalized distributions for different r should have collapsed to the same shape), and the analysis further confirmed the presence of deviations of the scaling exponents characteristic of turbulencelike phenomena.

V. WHOLE-BRAIN MODEL
We implemented a whole-brain model based on nonlinear Stuart-Landau oscillators coupled by the structural connectivity (SC), which can be globally scaled by a coupling strength parameter (G).The coupled dynamics are given by dx FIG. 2. Multifractility on whole brain dynamics.(a) Extended self-similarity confirms the inertial range for higher-order structure functions.We investigated the selected inertial range by computing the extended self-similarity based on studying the higher-order structure functions as a function of that of order 2. We found that the fit of the ESS is excellent for the odd structure functions, providing evidence that the selected inertial range is suitable to assess the spatial scaling properties at the different orders considered.(b) We computed scaling power laws based on high-order structure functions.We obtained the exponents of the even structure functions through the direct computation of the log-log linear fitting (green circles) and also through ESS (blue triangles).We found a high level of agreement between both measures, which confirms inertial range scaling.Interestingly, the exponents can be organized linearly up to order 4 and then the slope changes and the exponent saturates.This is compatible with at least a bifractal structure with two characterstic scales.(c) We investigated the distribution of the structure function of order 1 [S 1 (r)] at different scales within the inertial range.We found that, whereas the scale increases, the distributions of S 1 (r) tend monotonically to a Gaussian distribution, indicating small scale fluctuations are less common.
where x j is the dynamical variable that simulates the fMRI signal of the region j, ω j is the corresponding empirical frequency, and SC i j represents the symmetrical coupling matrix that weights the connectivity between regions i and j, (Appendix C contains a detailed explanation of the model, and how the the structural connectivity matrix is obtained is detailed in Appendix A).In the model, ν j stands for an additive Gaussian noise to node j with amplitude β.We fixed the bifurcation parameter equally for all brain regions close to the bifurcation point (a = −0.02),and we explored G from 0 to 3 in 0.1 steps to find the optimal working point of the model in two scenarios: only considering the second-order structure function, and including information from higherorder structure functions.Notably, the second-order structure function is, essentially, the functional connectivity (probably the empirical observable most frequently reproduced in the whole-brain modeling literature [23,28]), while the inclusion of higher-order information is one of the important aspects of this work.We computed the fitting of the model as the Euclidean difference between the empirical and simulated second-order structure functions within the inertial range, as done previously [9][10][11]16].The fitting of the higher-order functions was computed as the distance between the empirical and simulated scaling exponents for the even functions.Note that for this analysis, we only considered the even exponents, based on the fact that we could compute them directly from the structure function (see Fig. 1) and the previous convergence analysis (see Appendix B).Each optimal working point was determined as the value of G that minimizes the fitting function in each case.Appendix C contains more details of the fitting for each model.
We simulated 100 trials with the optimal working point for the second order (G 02 op = 1.26) and 100 trials for the optimal working point for higher orders (G HO op = 1.44).We found statistically significant differences when fitting the model to the second-order structure function in both working points (Wilcoxon signed-rank test p = 0.049 and effect size rank-biserial correlation = 0.14) [see Fig. 3(a)].We then assessed which model more faithfully reproduced the information transmission properties of the empirical data in terms of an information cascade measure.This measure was introduced in previous works to quantify information transmission across scales in terms of the correlation between the local level of synchronization (see Refs. [9,10]).We found that the information cascades obtained when considering G HO op were higher compared to those obtained when considering G 02 op (Wilcoxon signed-rank test p < 0.001, rank-biserial correlation = 0.35) and also closer to the empirical values (G 02 op compared with empirical: Wilcoxon rank-sum test p < 0.001, effect size computed by the rank-biserial correlation = 0.18; G HO op compared with empirical: Wilcoxon rank-sum test p = 0.005, rank-biserial correlation = 0.08) [see Fig. 3(b)].
We then evaluated the model responses to external perturbations when set up in both optimal working points.We modeled the external perturbation by randomly changing the local bifurcation parameter of each region in the range [−0.02 : 0].Note that this perturbation is carefully defined to keep the dynamical scenario in the subcritical regime of each oscillator [29].We quantified the response to the in silico perturbation by measuring the susceptibility and information-encoding capabilities.To quantify the susceptibility, we adapted the proposal of Daido, who defines the FIG. 3. Higher-order information improves model capabilities.We found for each case the optimal working point as the value of G that minimizes the fitting function in each case.We found statistical difference between both working points by comparing the fitting of the models (a).We found that the information cascade obtained when the brain is modeled with G HOop is closer to the empirical, while the information cascade modeled with G 02op is significantly lower (b).We quantified the responses when models are perturbed through the susceptibility and the information-encoding capability.We found that both measures are significantly higher in 50 trials of the perturbation when the model includes the information of higher-order structure functions.*** means p < 0.001; * means p < 0.05 and p > 0.01 Wilcoxon rank-sum test.susceptibility of a large population of coupled oscillators as the variation of the Kuramoto order parameter under external perturbation [30].The information-encoding capability measures the ability of the system to encode external inputs; in this sense, it can be related to complexity measures such as Lempel-Ziv [31] (see Appendix D for a detailed explanation, and see Refs.[9,10,16]).We show in Figs.3(c) and 3(d) the computation of susceptibility and information-encoding capability when the model is perturbed 100 times for each optimal working point (G 02 op and G HO op ).We found that, when the model is built using higher-order information, both measures increase compared with the second-order working point (susceptibility: Wilcoxon signed-rank test p = 0.02, rank-biserial correlation = 0.23; information-encoding capability: Wilcoxon signed-rank test p = 0.005, rank-biserial correlation = 0.28).These results provide insight into the relevance of the higher-order structure function information to improve the models performance.

VI. DISCUSSION
We have shown that the whole brain dynamics can be reproduced by scaling laws akin to those in Kolmogorov's phenomenological scaling theory for an inertial range in turbulence.We extended previous results, which only considered second-order structure functions, to higher-order structure functions, and we noticed that the relationship between the order (p) and the scaling exponents [ζ (p)] depart from linear dependence [e.g., ζ (p) = p/3 in Kolmogorov's case].This result suggests that the brain dynamics can display at least a bifractal structure.The behavior is reminiscent of the case of fluids following the Burgers equation, where the departure from monofractality and saturation of the scaling exponents is associated with large-scale intermittency in a disordered system [27].This observation is also aligned with findings of multifractality in brain dynamics.Moreover, our analysis shows that including higher-order statistics indicators in whole-brain models of coupled oscillators allows for a better description of information transmission and reactivity.In particular, we showed that the inclusion of higher-order information narrows the configuration of whole-brain models and allows them to reproduce more faithfully important brain functions such as information processing and brain responsiveness.
It is important to note that this work does not propose a one-to-one mapping between the physical variables of a turbulent system and brain activity measurements, and one should be careful when conceptualizing brain dynamics as a system that exhibits turbulence in a classic physical way.Instead, what the results show is that the mathematical formalism developed to capture the phenomenology of turbulence can also be used to represent important aspects of whole-brain dynamics.
Turbulence is an example of an out-of-equilibrium system that displays universal scaling behavior.Recent studies have considered the possibility that the scale invariance observed in turbulence is the result of the existence of a critical point at infinite Reynolds number [32] or have found evidence of critical points at finite Reynolds numbers as other control parameters are changed [33,34].Also, the dissipative scales in turbulent flows have been found to be compatible, under some conditions, with self-organized criticality [35].However, and as previously mentioned, it is important to note that perfect self-similarity in turbulence is broken as a result of a a phenomenon known as intermittency.This phenomenon, which is quantified by the nontrivial dependence of the scaling exponents with their order, by measurements of multifractality, or equivalently, by deviations from Gaussian statistics, is associated to the forced-dissipative out-of-equilibrium nature of the system.In other words, it is caused by the irreversibility of turbulence and by the existence of a nonzero flux of energy across scales [3,33].Our proposed framework can be used to distinguish in brain dynamics measurements between criticality as found in equilibrium systems (which are, in most cases, scale invariant) and out-of-equilibrium systems (which, as in the paradigmatic case of turbulence, display deviations from scale invariance associated with their forced and dissipative nature).
The extension of Kolmogorov's theory in fluids to higherorder structure functions, and the deviations from perfect self-similarity, has been extensively challenged by experiments and numerical investigations based on the rate of energy dissipation being intermittent, i.e., spatially inhomogeneous, and thus that it cannot be treated as a constant [17,36,37].
Our results can, therefore, be interpreted as signatures that the turbulentlike regime of brain dynamics also presents intermittency or multifractality.Another possible interpretation is that the existence of such power laws does not demonstrate the existence of turbulence in the brain (understanding turbulence as a general out-of-equilibrium multiscale nonlinear process).However, for many natural systems the existence of universal scaling laws for different orders of the statistical moments, and the usage of turbulence indicators to characterize them, has been extensively verified and used as a way to characterize nonequilibrium dynamics [18][19][20].It is worth noting that, in recordings of brain activity data, several studies have also found power laws in the context of criticality [4][5][6]38], which also could be consistent with such a class of processes.Studies of the famous sandpile model proposed by Bak et al. [39] to capture the behavior of self-organized criticality have shown that the moments of wave size differences present similar scaling to those holding for velocity structure functions in fluid turbulence [8].The results presented here push this study to include many different statistical orders, allowing for simultaneous characterization of scale invariance as well as of the deviations from this behavior.
The preprocessing of the HCP resting state is described in details on the HCP website.In short, the data are preprocessed using the HCP pipeline which is using standardized methods using FSL (FMRIB Software Library), FREESURFER, and the CONNECTOME WORKBENCH software [42,43].This preprocessing included correction for spatial and gradient distortions and head motion, intensity normalization and bias field removal, registration to the T1-weighted structural image, transformation to the 2-mm Montreal Neurological Institute (MNI) space, and using the FIX artifact removal procedure [43].The head motion parameters were regressed out and structured artifacts were removed by ICA + FIX processing (independent component analysis followed by FMRIB's ICA-based X-noiseifier [44]).Preprocessed time series of all grayordinates are in HCP CIFTI grayordinates' standard space and available in the surface-based CIFTI file for each participant for the resting state.
We used a custom-made MATLAB script using the ft-readcifti function (FIELDTRIP toolbox [45] to extract the average time series of all the grayordinates in each region of the Schaefer parcellation, which are defined in the HCP CIFTI grayordinates' standard space.Furthermore, the BOLD time series were transformed to phase space by filtering the signals in the range between 0.008 and 0.08 Hz and the lowpass cutoff to filter the physiological noise, which tends to dominate the higher frequencies [46].

Structural connectivity using diffusion MRI
To obtain the structural connectivity we use the HCP database which contains diffusion spectrum and T2-weighted imaging data from 32 participants (acquisition parameters are described in detail on the HCP website).Briefly, the data were processed using a generalized q-sampling imaging algorithm implemented in DSI studio [47].The white-matter mask is produced by the segmentation of the T2-weighted anatomical images.The co-register of the images to the b0 image of the diffusion data was carried out using Statistical Parametric Mapping tool (SPM) [48].In each HCP participant, 200 000 fibers were sampled within the white-matter mask.Fibers were transformed into MNI space using LEAD-DBS [49].We used the standardized methods in LEAD-DBS to produce the structural connectomes for Schaefer 1000 parcellation [40] where the connectivity has been normalized to a maximum of 0.2.The freely available LEAD-DBS software package [50] provides the preprocessing that we implemented and is described in detail in Horn and colleagues [51].

APPENDIX B: CONVERGENCE OF THE ACCUMULATED MOMENTS
We have analyzed the convergence of the different accumulated moments following the work of Gotoh and colleagues [26].We displayed the convergence of the orders 2 to 8 accumulated moments Cn (n = 2 to 8) for various separations r p = 22 + 2p, p = 1, 2, 3, 4, 5, and 6.Curves are for p = 1, . . ., 6 and all values are within the inertial range (see Fig. 4).In Fig. 4, the insets in the odd orders are the accumulated moments computed with the absolute value.

APPENDIX C: WHOLE-BRAIN MODEL FITTING PROCEDURE
1. Second-order and higher-order structure functions model fitting The second-order structure function [S 2 (r)] is defined as the N×r matrix of S 2 , where u is the BOLD signal and r is the distance between brain regions (N).We computed the empirical S 2 (r) for each human participant and for each simulated trial (the total number of trials matched the number of participants) and computed the Euclidean distance within the inertial range: where A and B are the lower and upper limits of the considered inertial range.On the other hand, the higher-order structure functions were computed as S p (r) = U (r) p = [u(x + r) − u(x)] p , for p = 1 to p = 8.We computed the empirical S p (r) for each human participant and extracted the spatial scaling power law for each even order.We replicated the same computation for the model signals and then compared the empirical and modeled scaling exponents.We quantified the fitting by computing S fitt HO (r) = (C2) We explored the coupling strength parameters from 0 to 3 in 0.1 steps to find the optimal working point of the model when fitting the second-order structure function and when the higher-order structure functions were considered.We computed the fitting of the model as the Euclidean difference between the empirical and simulated second-order structure function within the inertial range.The higher-orders' fitting was computed as the distance between the even empirical FIG. 5. Whole-brain model fitting to second-order structure function and higher-order structure functions.
scaling exponents and those obtained with the model.We found for each case the optimal working point as the value of G that minimizes the fitting function in each case (Fig. 5).

Susceptibility
We computed the sensitivity of the perturbations on the spatiotemporal dynamics by extending the definition of previous work, which determines the susceptibility in a system of coupled oscillators based on the response of the Kuramoto order parameter [30].The Hopf model was perturbed for each G by randomly changing the local bifurcation parameter (a) in the range [−0.02 : 0].The sensitivity of the perturbations on the spatiotemporal dynamics was calculated by measuring the modulus of the local Kuramoto order parameter as where R(m) n ( x, t ) corresponds to the perturbed case; R (m) n ( x, t ) corresponds to the unperturbed case; and t , trials , and x correspond to the average across time, trials, and space, respectively.

Information-encoding capability
The information-encoding capability measures the ability of the system to encode external inputs, and as such is more closely related to complexity measures such as Lempel-Ziv or the automatic complexity evaluator, and synchrony coalition entropy (used and defined in Ref. [3]).The information capability I was defined as the standard deviation across trials of the difference between the perturbed R(m) n ( x, t ) and the unperturbed R (m)  n ( x, t ) mean of the modulus of the local Kuramoto order parameter across time t, averaged across all brain areas n as where the brackets t , trials , and x denote the averages defined as above.

FIG. 1 .
FIG. 1. High-order structure functions and their spatial power scaling laws.(a) We computed the structure functions from 1 to 8 by replacing the longitudinal component of the fluid velocity in Eq. (2) by the BOLD signal from fMRI recordings of 1003 healthy participants.We highlight even structure functions in the inertial subrange.(b) We found the spatial power-law scaling of even structure functions S(r) as a function of log(r) within the inertial subrange.

FIG. 4 .
FIG. 4. Convergence of the accumulated moments for various separations.