University of Birmingham Characterisation of haemodynamic activity in resting state networks by fractal analysis

• Users may freely distribute the URL that is used to identify this publication. • Users may download and/or print one copy of the publication from the University of Birmingham research portal for the purpose of private study or non-commercial research. • User may use extracts from the document in line with the concept of ‘fair dealing’ under the Copyright, Designs and Patents Act 1988 (?) • Users may not further distribute the material nor use it for the purposes of commercial gain.

Intrinsic brain activity is organized into large-scale networks displaying specific structural-functional architecture, known as resting-state networks (RSNs). RSNs reflect complex neurophysiological processes and interactions, and have a central role in distinct sensory and cognitive functions, making it crucial to understand and quantify their anatomical and functional properties. Fractal dimension (FD) provides a parsimonious way of summarizing self-similarity over different spatial and temporal scales but despite its suitability for functional magnetic resonance imaging (fMRI) signal analysis its ability to characterize and investigate RSNs is poorly understood. We used FD in a large sample of healthy participants to differentiate fMRI RSNs and examine how the FD property of RSNs is linked with their functional roles. We identified two clusters of RSNs, one mainly consisting of sensory networks (C1, including auditory, sensorimotor and visual networks) and the other more related to higher cognitive (HCN) functions (C2, including dorsal default mode network and fronto-parietal networks). These clusters were defined in a completely data-driven manner using hierarchical clustering, suggesting that quantification of Blood Oxygen Level Dependent (BOLD) signal complexity with FD is able to characterize meaningful physiological and functional variability. Understanding the mechanisms underlying functional RSNs, and developing tools to study their signal properties, is essential for assessing specific brain alterations and FD could potentially be used for the early detection and treatment of neurological disorders.

Introduction
In the last two decades, ongoing brain fluctuations at rest (i.e. in the absence of external stimulation and response demands) have consistently been documented to be organized into large-scale networks, 1-4 each of them characterised by specific structural and functional architectures. [5][6][7][8] These are known as resting-state networks (RSNs), and have been widely reported in numerous neuroimaging 5,6 and electrophysiological 9-12 studies. The spatial pattern of RSNs is stable within individuals, consistent across healthy subjects and yet also contains considerable characteristic information about individuals' brain function [13][14][15] making them a powerful tool for neuroscientific enquiry. Furthermore, alterations to RSNs have been observed during healthy aging 16 as well as in many neuropsychiatric 17,18 and neurological disorders. [19][20][21] Understanding RSNs as an organizational principle of brain function is therefore of crucial importance. RSNs are generally divided into sensory (e.g. visual, auditory) and cognitive (e.g. default mode, fronto-parietal) networks based on their anatomical and functional properties. In this study, we hypothesise that different RSNs studied using the Blood Oxygen Level Dependent (BOLD) technique might display different brain dynamics, as already shown using electroencephalography (EEG) data by our group. 9 Specifically, we aimed to investigate how the dynamics of BOLD activity at rest can be used to differentiate RSNs and link them to behavioural and perceptual states.
However, despite the fact that linear methods are predominantly used in characterizing brain oscillations in both healthy and pathological conditions, linear analysis may not be suitable to describe the irregular and nonperiodic patterns recorded by electrophysiological and neuroimaging techniques. 22,23 To this end, we characterised the specific BOLD signature of each RSN, using a complexity measure called Fractal Dimension (FD) 24 that has advantages over classical linear methods such as the wellknown fast Fourier transformation (FFT) that are best suited to conditions where the analyzed signals are stationary. FD is a general measure of complexity derived from chaos theory, based on the fact that a simple process that is repeated endlessly becomes a very complex process, which is the basis for the description of fractals in nature. 25,26 These complex processes of interactions cause a pattern in brain activity that is self-similar over different spatial and temporal scales. In other words, neural activity shows similar features over and over again on a scale-free basis. 27,28 Knowing that FD is an accurate numerical measure no matter what the properties (stationary, nonstationary, deterministic or stochastic) of the analyzed signal, it is reasonable to accept this advantage over widely used FFT based or other linear methods. 29 In addition, recent evidence has demonstrated that in many cases brain signals considered as belonging to a frequency-defined class of brain rhythms do not represent sustained oscillations, but rather brief bouts of activity that are repeated intermittently (i.e. nonrhythmic). 30,31 Recognition that physiological time series contain "hidden information" that might be captured by nonlinear methods such as fractal analysis (FA), may provide crucial and so far overlooked physiological information in healthy and pathological conditions. [32][33][34][35][36][37][38][39] We have shown that FD can be used to dissociate RSNs from EEG data 9 but, to our knowledge, no study has used FD in functional magnetic resonance imaging (fMRI)-derived RSNs. To fill this gap, we used a data-driven approach based on hierarchical clustering to differentiate RSNs based on their functional roles. 40 This approach has the advantage of containing very few assumptions about the data structure, allowing us to investigate whether FD as a feature of the data is able to map onto the functional roles of the RSNs. We were able to show that complexity features derived from BOLD permit the separation of RSNs into two main clusters, one mainly consisting of sensory networks, the second was more related to higher cognitive (HCN) functions.

Subjects
Fifty-two right-handed volunteers (age = 25.2 ± 3.5 years, 23 females) participated in a resting-state fMRI experiment, during which they were asked to lie still, keep their eyes open and think of nothing in particular. The protocol was approved by the Research Ethics Board of the University of Birmingham, and written informed consent was obtained from all participants.

fMRI data acquisition and preprocessing
The experiment was conducted at the Birmingham University Imaging Centre using a 3T Philips Achieva MRI scanner. An eight-channel phasedarray head coil was used to acquire T1-weighted anatomical image (1 mm isotropic voxels) and wholebrain T2*-weighted, functional EPI data ( . Pre-processing of the data included slice time correction, spatial realignment to correct for movement artefacts and motion by distortions interactions, and normalization to the MNI standard space and spatial smoothing with a Gaussian kernel of 6 mm full-width at half-maximum. To ensure consistency of results across the group, only the first 150 volumes of data were analyzed for each participant.

fMRI data analysis
After preprocessing, resting state fMRI data of all participants were analyzed using spatial independent component analysis (ICA) as implemented in the Group ICA of fMRI Toolbox (GIFT -http://trendscenter.org/software/gift/). GIFT decomposed the data into functional networks that exhibited a unique time course profile within each scan. Two data reductions steps were carried out using principal component (PC) analysis, subject-specific and group-level steps. 41 First, subject-specific data were reduced to the 30 PCs explaining maximum variance and then these subject-reduced data were temporally concatenated. For the subject-specific data we have selected 30 PCs following the rule that in this step more PCs need to be estimated respect to the number of independent components (ICs). 42,43 In particular, we have used the number of ICs (#ICs) estimated by the minimum description length (MDL), see below, and the following formulation: #PCs = #ICs+1/2(#ICs). Further, at the group level, data were reduced to 20 group ICs with the expectation-maximization algorithm, included in GIFT. 44 The number of ICs was estimated using the MDL criterion. 45,46 In our specific case, 20 ICs were indicated to be estimated using MDL. Subject-specific spatial maps and time courses were obtained using the back-reconstruction approach (GICA). 45 These time courses were then used in the subsequent FD analysis (see below). From the 20 ICs, we identified the relevant RSNs by applying a previously described procedure. 44 We first manually confirmed if the peak activation coordinates were located primarily in the gray matter, showing low spatial overlap with vascular, ventricular, and edge regions corresponding to artefacts. 44 This process resulted in eleven ICs representing meaningful RSNs (see Fig. 1 and Table 1) that together comprised six large-scale functional  networks, based on the spatial correlation between ICs and the template provided by GIFT Toolbox. 44 The RSNs were identified as: ventral default mode network (vDMN -IC 15) and both anterior and posterior regions of the dorsal DMN (adDMN -IC6; pdDMN -IC19); primary visual (PV -IC16) and higher visual network (HVN -IC14); ventral sensorimotor network (vSMN -IC9) and dorsal SMN (dSMN -IC10); right and left frontoparietal networks (rFPN -IC8 and lFPN -IC18); auditory network (AN -IC12) and anterior Salience Network (aSN -IC13).

Fractal analysis
A fractal is a shape that retains its structural detail despite scaling and this is the reason why complex objects can be described with the help of FD. 29,47,48 FD is a highly sensitive measure for the detection of hidden information contained in physiological time series. 29,49 In the time domain, FD represents the amplitude of the signal under consideration on a plane as a function of time. As a consequence of this planar geometry, its dimension is limited between 1 and 2, since a simple curve has dimension equal to 1, while a plane has dimension equal to 2. FD will be equal to 2 for time series with high complexity and no memory of the signal (i.e. the overall signal cannot be predicted using an individual portion of the signal), such as white random noise, while it will be equal to 1 for time series with low or no complexity in the signal, such as a ramp line, which is characterised by full memory (i.e. the overall signal using can be predicted using a portion of the signal itself). Moreover, it has been demonstrated that higher values of FD correspond to the presence of higher frequencies in the signal's Fourier spectrum and vice versa. [50][51][52] This relationship might not be linear: in our previous EEG study, 9 we showed that it was rather quadratic.
There are many methods used to calculate FD, but the widely accepted ones are Katz's and Higuchi's methods, respectively, 24,53 with the latter considered the most accurate to estimate FD. 54

Higuchi's fractal dimension
Higuchi's FD 24 is a quantitative measure of signal dynamics. Linear methods commonly used for signal analysis, such as the well-known FFT and wavelet transformation (WT), are good choices if the analyzed signals are stationary. However, neurophysiological processes are generally nonstationary and nonlinear by nature. Knowing that FD is an accurate numerical measure no matter what the properties (stationary, nonstationary, deterministic or stochastic) of the analyzed signal, it is reasonable to accept this advantage over widely used linear methods. 29 Higuchi's FD is a nonlinear measure of waveform complexity, here applied in the time domain. Discretised functions or signals can be analyzed as segment of data X(1), X(2), . . . , X(N ), where N is the total number of samples. From the starting time sequence, a new self-similar time series X k m can be calculated as The length, L m (k), of each curve X k m is calculated as k . In practice, the original curve or signal can be divided into smaller parts with or without overlap, called "windows". Then, the method for computing FD should be applied to each window where N should be seen as the length of the window. In that case, FD values are calculated for each window, with or without overlap. Individual FD values can be averaged across all windows for the entire curve (or data timeseries), and the mean FD value can be used as a measure of curve complexity.

An array of mean values L(k) was obtained and the FD was estimated as
Here, using the single-subject IC timecourses for each RSN, we calculated FD in nonoverlapping time windows of 100 s (corresponding to 50 of our fMRI volumes) as a good compromise between windows length, length of the data and computational time. The choice of the free parameter k has a crucial role in FD estimation. For each window we estimated twenty-four values of FD for k = 2, . . . , 25. The value 25 was equal to half of the samples within our 100 s window (i.e. 50 volumes and the maximum that can be chosen, k max is equal to half of the window length). There were three windows within our 150 volume scans, therefore we estimated three measures of FD at each value of k (e.g. FD 2 , FD 3 , FD 4 . . . FD 24 ). These three measures were averaged to give one mean value of FD for each k, for each subject. 29 The process was then repeated for every subject and every RSN.
Additional analysis demonstrated that the FD measurements were not dependent on the choice of window length or overlapping windows 29 (see supplementary figures in Appendix A -Figs. S1-S3). For the subsequent FD analysis, we set k = 12, which corresponds to the value just above the median. 9,29,38

Hierarchical clustering
Hierarchical clustering was used to subdivide the 11 ICs based on the similarity of the FD values derived from their BOLD signal. The relationships between ICs are represented by a hierarchical tree, in which branch length reflects the degree of similarity between the RSNs, as assessed by a pairwise similarity function, between the ICs. Notably, the separation of ICs into groups might also depend on the metric used to measure the distance between networks and on the linkage criterion, which specifies the dissimilarity of networks. Here, the Euclidean distance of FD values estimated among the ICs was used as similarity function for clustering the ICs similar to our previous work. 9 Given an m-by-n data matrix X, which is treated as m (1-by-n) row vectors x 1 , x 2 , . . . , x m , the various distances between the vectors x s and x t are defined as follows: T and the inner squared distance calculated by the minimum variance algorithm was used to measure the distance between clusters. 9 The cluster analysis in this study was applied as shown in the Statistics and Machine Learning Toolbox (https://it.mathworks.com/help/stats/hierarchicalclustering-12.html).

Classification of RSNs based on BOLD FD: Single-subject level
Two clusters were obtained by hierarchical clustering analysis: C1 (AN, PV, HVN, vSMN, dSMN,  aSN, and vDMN and C2 (lFPN, rFPN adDMN, and  pdDMN). We then investigated FD at the single subject level to understand the robustness of the classification and the consistency with which FD could identify the type of network, by defining a similarity index for the comparison between the C1 and C2 clusters, i.e. whether the FD feature was higher or lower for C1 compared to C2 at the single subject level. The similarity index was calculated as (C1−C2)/(C1+C2). This index is defined between −1 and +1, and was equal to zero only when the values of C1 and C2 are identical.

Statistical analysis
To test for FD differences between each pair of RSNs, two-sample permutation t-tests (10,000 permutations; p < 0.05) were performed. In case of the point-by-point two-sample permutation t-tests, the false discovery rate (FDR) was used to correct for multiple comparisons. 55

fMRI-Derived RSNs: Group level characterization by FD
FD was calculated on the BOLD signal timecourse of each RSN. Hierarchical cluster analysis ( Fig. 2 -Top, see Fig. S4 for a comparison with fractional amplitude of low-frequency fluctuations (fALFF) analysis) using the BOLD FD with k = 12 produced a dendrogram with two clusters: C1 (blue : AN,  PVN, HVN, dSMN, vSMN, vDMN, and aSN) and C2 (red: adDMN, pdDMN, lFPN, and rFPN). C1 mainly consisted of sensory RSNs (with the exception of vDMN and aSN), whereas C2 was composed of HCN RSNs (DMN and FPN). The group average FD for each spatial map is also shown (Fig. 2 -Middle). Mean and standard error (SE) for each RSN (blue and red) and for the mean of C1 and C2 clusters (black) are also shown, indicating a significant difference between C1 and C2, with higher FD for C1 with respect to C2 (p < 0.01). Finally, we analyzed FD values of the different RSNs (Fig. 2 -Bottom). This demonstrated that there were differences among RSNs (e.g. vSMN higher FD than any other  FPN (lFPN and rFPN)) (**indicate p < 0.01). For the Permutation Test all the T-stat values that did not reach the significance level (p = 0.05) were set to zero. RSN), but it was not enough to discriminate each single RSN and generally FD was consistent within each cluster. Since the choice of the free parameter k has a crucial role in FD estimation, we performed pointwise (for k = 2, . . . , 25) two-sample permutation t-tests (10,000 permutations, p < 0.05, FDR corrected) between C1 and C2. We found that the values from C1 were higher than those from C2 for all values of k (see Fig. 3) and that the clustering results were consistent at different values of k. This validated our choice of using k = 12 to report all main results. In particular, for k = 12 the permutation t-test and T -values were: p = 0.0002, t = 3.886.

Classification of RSNs based on BOLD FD: Single-subject level
The similarity index for C1 and C2 (Fig. 4) showed that for 50 subjects out of 52 the FD derived from C1 was significantly higher than that from C2, which led to 96.15% (50/52) accuracy at distinguishing C1 from C2. This finding suggests that it is possible to distinguish the C1 RSNs from C2 RSNs based on  their BOLD activity using FD. More specifically, an RSN characterised by a higher FD value is likely to belong to the C1 cluster with respect to the C2 cluster. This emphasises that FD can be used as a robust and parsimonious marker to define RSNs in pathological conditions.

Discussion
In this study, we have investigated how Higuchi's FD as a measure of multi-scale signal complexity can be used to characterise and differentiate the resting BOLD fMRI signal of RSNs in 52 healthy participants. As FD requires very few assumptions to be made about the underlying structure of the data (e.g. linearity, stationarity), it is an attractive tool for analyzing neuroimaging and neurophysiological data. FD is also a natural way to quantify complexity, and fits naturally with the current conceptual framework that models the brain as a complex system. 56,57 Using FD, we were able to identify two clusters of RSNs, one mainly consisting of perceptual networks (C1, including auditory, sensorimotor and VN) and the other associated with HCN functions (C2, including dorsal default mode network and frontoparietal networks). These clusters were defined in a completely data-driven manner using hierarchical clustering of the FD of RSNs BOLD signal, suggesting that quantification of BOLD signal complexity is able to identify meaningful physiological and functional variability between RSNs. The differentiation of these networks is consistent with previous EEG 9 and fMRI 40 studies, which found a similar dissociation between perceptual (PN) and HCN networks.
In our previous study 9 we followed a similar analysis approach, but we used 256-channel EEG data instead of fMRI data. Comparing spectral analysis, Shannon entropy and FD, we demonstrated that FD was able to provide a clearer dissociation between RSNs than the other signal analysis approaches. A very similar overall picture of lower FD in HCN compared to PN was found (see Figs. 2 and 3). The fact that FD can provide similar differentiation of RSNs and with the same description of complexity when based on either EEG or fMRI data is reassuring, and points to its sensitivity to underlying physiological variation. Ding and colleagues, 40 working with fMRI data, differentiated PN from HCN on the basis of network topology, again supporting the idea that these networks have different functional properties. Wang and colleagues 58 were also able to cluster RSNs based on quantification of entropy from fMRI data, although they did not provide an explicit comparison between HCN and PN. More work is needed to understand how complexity measures are related to functional connectivity and network topology 59 but the available data supports the distinction between HCN and PN on the basis of signal properties alone.
Furthermore, in our previous research we found a positive correlation between FD and gamma power, and a strong negative correlation between FD and theta power. 9 Given the putative relationship between gamma band activity and BOLD signal, 60,61 our findings suggest consistencies between EEG and fMRI results. In both cases, while using different neuroimaging modalities, PN showed higher FD than HCN. Whether this reflects the resting state data limiting the functional repertoire of HCN, or is a fundamental property of the networks, remains to be investigated in further studies. For example, Cottone and colleagues 52 demonstrated an increase in EEG signal FD with task performance in both motor and sensory cortices, suggesting that a more active state is reflected in increased complexity. Using multiscale entropy, McDonough and Nashiro, 59 suggested an interaction between complexity and time scale such that, for example, the complexity of DMN activity was low at short time scales and high at long time scales. Our data (Fig. 3) would not support this behavior with FD, which may represent a difference between the two measures of complexity and the fact that FD summarizes complexity over all scales. The clustering of the aSN with PN may also be related to the resting state. The aSN dynamically controls changes of activity in other networks, by mediating sensory and cognitive processing required for executive control. 62 In particular, it has been found to play a role in coordinating the activity between task-related networks and executive control networks, and might be particularly involved in the rapid and efficient engagement needed for motor control. 63 Dynamic analysis of FD, similar to ongoing work on dynamic functional connectivity, 64,65 as well as comparison of active vs resting state, 66 will be needed to understand more about the utility and interpretation of FD. This future work will benefit from longer fMRI scan durations than we used, in order to ensure that more prolonged patterns of self-similarity can be uncovered (e.g. using the one hour of fMRI data available as part of the Human Connectome Project). 67 Even longer time scales are only likely to be accessible using ambulatory EEG, but as we have demonstrated here and in our previous work 9 FD is a metric which allows direct comparison across modalities. In addition, given the complexity of the relationship between the BOLD signal and electrophysiology, future work will be needed to understand the neuronal and metabolic contributions to the fractal behaviour that we have used. 68 While the majority of the networks included in C1 were PN, and those in C2 were HCN, there were apparently contradictory specific cases (i.e. where an RSN included in C1 was more obviously a HCN). In particular, the ventral and dorsal DMN subsystems were assigned to different clusters, with the vDMN being more closely associated with the PN. The dissociation between DMN sub-systems based on functional connectivity has been widely demonstrated, [69][70][71][72][73] and it is clear that they have different functional roles. 74 For example, dorsal and ventral posterior cingulate cortices (PCC), which are a primary component of the vDMN/dDMN we investigated, are differentially impacted by task difficulty in a working memory task. 73 Dorsal PCC has increased integration with the other DMN regions as task difficulty increases, and a stronger anti-correlation with cognitive control networks. Ventral PCC shows the opposite response to increased task difficulty. However, it has not been suggested that either of these sub-divisions serves basic PN functions, and hence it is not at this stage clear why from a FD/complexity point of view the vDMN would cluster with PN. Ding and colleagues 40 did not examine sub-divisions of the DMN, so whether a similar differentiation would be observed from the point of view of network properties is unknown. The brain is a complex system close to criticality, which leads to the conclusion that characterization and quantification of complexity is likely to be an important method for understanding it. Higuchi's FD is an efficient and accurate method for that purpose, and we were able to demonstrate that it can dissociate RSNs in a way that is generally consistent with their function. FD represents a relatively straightforward index of complexity across multiple scales, which can summarize the high dimensionality data routinely acquired as part of neuroimaging experiments. We demonstrate consistency between classification of RSNs using fMRI data, and a previous EEG study, 9 supporting the idea that it is sensitive to underlying, physiologically meaningful variability.