1. Introduction
Among the most intricate networks in nature is the human brain, which transports signals between specific brain regions and responds to external stimuli. The study of brain connectivity, therefore, depends to a great extent, on the comprehension of brain functions and pathology (
Sporns & Zwi, 2004).
Functional Magnetic Resonance Imaging (fMRI) is an imaging technique applied to study human brain function and neurological diseases (
Smith, Matthews, & Jezzard, 2001). Blood Oxygenation Level Dependent (BOLD) contrast, on the basis of varied magnetic properties of oxygenated (diamagnetic) and deoxygenated (paramagnetic) blood, is applied by well-liked techniques in fMRI (
Fransson, 2005). fMRI benefits over the other functional imaging modalities, namely Electro Encephalography (EEG), Magneto Encephalography (MEG) and Positron Emission Tomography (PET) comprise its non-invasiveness, better spatial resolution compared to other modalities and short-time image attainment (
Culham & Kanwisher, 2001).
In literature, the term “network” has various definitions. In graph theory and complex networks, ‘‘network’’ bluntly implies a set of nodes and pair-wise edges, by which the nodes are connected. This sense is referred to as Graph theoretical analysis of brain networks. ‘‘Network’’, in neuroimaging, may designate a group of voxels or Regions of Interest (ROIs), that at resting state or in specific cognitive tasks, act identically (
Sun, Tong, & Yang, 2012).
From the standpoint of complex networks, there has been a growing interest, over the recent years, in studying the wide-ranging brain activity interaction structure with the use of graph theory and centered on fMRI in the areas of addiction (
Sutherland, McHugh, Pariyadath, & Stein, 2012), schizophrenia (
Cabral, Kringelbach, & Deco, 2012), brain injury (
Nakamura, Hillary, & Biswal, 2009), neuralgia (
Zhang et al., 2014), epilepsy (
Ponten, Bartolomei, & Stam, 2007), Alzheimer
(Supekar, Menon, Rubin, Musen, & Greicius, 2008) and the like (
Bullmore & Sporns, 2009).
A brain network, in graph theoretical analyses of fMRI data, is taken into account as an undirected graph, G=(V, E), where a node/vertex (V) in the graph delineates a brain region (i.e. ROI) and an edge/link (E) between two nodes is indicative of brain regions being functionally connected (
Deuker et al., 2009).
To create a graph, the concept and also the definition of the edge is a challenging stage. The edge definition methods and the importance of each is well evaluated and tested in Smith study (
Smith et al., 2011).
However, to assess the interaction strength between two brain regions for the edge definition, the LC (Pearson) coefficient of corresponding time series is most commonly utilized in functional connectivity studies, substantially in addiction research studies
(Xia & He, 2011). The drawback of such a choice may be that the linear correlation does not take into account the nonlinear dependences possibly occurring in the data. While linear measures, including the Pearson correlation coefficient or coherence are frequently applied, increased attention is being paid to potential benefits of nonlinear measures
(Donges, Zou, Marwan, & Kurths, 2009;
Kreuz et al., 2007). Mutual Information (MI) can establish a nonlinear relationship between fMRI time series and provides an effective pure nonlinear noise-robust correlation measure (
Bassett & Bullmore, 2009).
Research studies contemplating nonlinearity as an inherent feature of the brain dynamics are getting continuously more interested in nonlinear approaches toward the analysis of the brain signals, especially those measures in accord with the analysis of chaotic non-linear dynamical systems to analyze the resting state fMRI data, signifying that the presupposition of linearity might be oversimplifying (
Hlinka, Paluš, Vejmelka, Mantini, & Corbetta, 2011). On the other hand, the brain, with fractal structure complexity, is best modeled as a complex system (
Papo, Zanin, & Buldú, 2014). For the complexity of the brain signal, we can also assess the resting BOLD fMRI time series (
Warsi, 2012). The fMRI time series within any given voxel, rather than correlating the brain areas connectivity with the use of measures of linearity and nonlinearity, can reveal the resting state brain functional network complexity. By complexity we mean the spatial distribution of fMRI signals phase that should not be confused with brain complexity.
The most frequently applied method for analyzing the physiological signals complexity is fractal analysis (A
hmadi, Ahmadlou, Rezazade, Azad-Marzabadi, & Sajedi, 2013). Fractal Dimension (FD) analysis is likely to include new information on the functional connectivity of the brain (
Sporns, 2006). This information is not attained by applying any traditional linear, as well as nonlinear measures. However, it should be noted that selection of the appropriate method to determine the graph edges in functional connectivity studies and the impact of each
are defined and the brain network is constructed. Finally, statistical analysis is performed. Each phase is distinctly stipulated in the following sections.
2.1. Subjects
A total of 18 normal controls (NCs: 18, aged 23–46 years, mean±SD=31.67±7.98 years, right-handed) were enlisted from the local community. Additionally, we enrolled 17 age-matched methamphetamine-dependent individuals (MDIs: 17 Males, aged 22–39 years, mean±SD=30.52±4.57 years, right-handed). Demographic and clinical information regarding MDIs is demonstrated in Table 1. Mann-Whitney U test (
Bergfeldt, Jonsson, Bergfeldt, & Julin, 2015) showed no significant difference between MDIs and NCs (P=0.85).
All subjects were capable to realize and follow the protocol of the study during imaging. Neither MDIs nor NCs had any history of pain, neurologic or psychiatric disorders, head injuries, schizophrenia or an affective disorder, in line with their medical history. The present paper was performed in agreement with the Helsinki Declaration also, confirmed by the Research Ethics Review Board of Tehran University of Medical Sciences in Tehran, Iran. Prior to the MRI scanning, informed written consent was obtained from each subject.
2.2. Data acquisition
On a 3 Tesla Siemens Tim Trio scanner, MRI data were acquired in Medical Imaging Department at Imam Khomeini Hospital, Tehran, Iran. To acquire data for Resting-State Magnetic Resonance Imaging (RS-fMRI), a T2*-weighted gradient-echo Echo-Planar Imaging (EPI) sequence was used with the following parameters: Time of Repetition (TR)=3000 ms, Time of Echo (TE)=30 ms, flip angle=90°, matrix=64×64, field of view (FOV)=192 mm2, thickness/gap=4.5 mm, 22 axial slices covering the whole brain and 240 volumes attained in almost 8 minutes. Besides, via a T1-weighted 3D turbo-gradient-echo sequence (TR=1800 ms, TE=30 ms, flip angle=90°, matrix=256×256, FOV=230×230 mm2, thickness=1.0 mm, and 160 sagittal slices), brain structural images with 3D high resolution were obtained. The participants were supposed to have their eyes closed, keep calm, not to systematically think about anything, and not to fall asleep. None of the lights in the scanner room were on during the RS-fMRI scanning.
2.3. Data preprocessing
All the fMRI data were preprocessed using MATLAB software and DPARSF_V2.0 Toolbox (
Chao-Gan & Yu-Feng, 2010). For each individual, in the first place, the first 10 volume images were excluded from the RS-fMRI data to make the subjects adapt to the environment and to stabilize the scanner, letting 260 volumes remain for supplementary analysis. Then, to correct the acquisition time delay between slices within the similar TR, slice timing was administered. After that, for the correction of the inter-TR head motions, realignment to the first volume was done. Next, we spatially normalized the fMRI data to a standard Montreal Neurological Institute (MNI) (
Tzourio-Mazoyer et al., 2002) template and implemented resampling to a voxel size of 36×63 mm3. According to the studies conducted so far, no spatial smoothing was applied
(Achard & Bullmore, 2007;
Braun et al., 2012;
Wang et al., 2009). Eventually, to reduce the low-frequency drift and high-frequency physiological noise, band-pass filtering was performed for each voxel at the frequency of 0.01–0.08 Hz. The RS-fMRI data for each subject were checked for head motion. In accord with the criteria that the translation and rotation of head motion in any direction were not more than 1.5 mm or 1.5, none of the subjects were disqualified.
2.4. Graph construction
Topological properties of the brain were investigated through the method of binary graph G=(V, E), where a brain region (i.e. ROI) and an edge (E) between two nodes in the graph are indicated by a node/vertex (V). All the regarded graphs in this work are undirected. To analyze complex network, we applied a simple generalization of the graph called weighted graph
(Bollobás, 1998).
In this paper, based on the previously conducted studies, two common methods (i.e. LC (
Hlinka et al., 2011)and MI (
Hartman, Hlinka, Paluš, Mantini, & Corbetta, 2011)), as well as the proposed method (i.e. BCFD) were applied for edge definition. For graph construction, node and edge need to be defined, each of which is explained in detail as follows:
2.4.1. Node definition
In order to have the brain parcellated into 90 Regions Of Interest (ROIs) (45 in each hemisphere), an Automated Anatomical Labeling (AAL) atlas (
Tzourio-Mazoyer et al., 2002) was utilized to construct the brain functional networks for each participant.
2.4.2. Edge definition
At this stage, the connectivity, dependence and interactions between the brain regions, i.e. the graph nodes, have to be quantified. Edge definition phase to create graph functional connectivity is challenging to some extent, for which various methods, including linear and nonlinear are presented so far (
Ahmadlou & Adeli, 2011;
Ahmadlou, Ahmadi, Rezazade, & Azad-Marzabadi, 2013). It should be noted that the method considered as the best method for edge definition and as the edge definition gold standard is still a questionable issue (
Papo et al., 2014). According to the previous studies, we used LC and MI methods as the representatives of the linear and non-linear methods, respectively, together with the proposed method “BCFD” for quantifying the interactions between 90 regions selected as nodes based on AAL Atlas. Then, BCFD method was evaluated in comparison with LC and MI. To quantify the nodes connectivity, each of these three methods considers a different criterion.
LC method quantifies the interactions between the time series of each node with linearity supposition, while MI lacks the assumption of the linearity between the brain regions and proposes an index for the evaluation of the shared information between the two time series. But this method is sensitive to the length of time series being analyzed. Since, in functional connectivity studies, the length of time series is usually short, MI use can be somewhat challenging. BCFD was applied in the current study, as it is based on the concept of self-similarity and complexity, i.e. distribution of the points in the phase space, calculating the connectivity and interactions of the brain regions. In addition, it does not assume interactions linearity or sensitive to the length of the time series. The concepts of “self-similarity” and “complexity” have not been evaluated using the graph theory in the functional connectivity studies. We computed the time series for each ROI, by averaging the signals of all voxels within that region. By calculating the Pearson correlation coefficient, MI and BCFD in the residual time courses between all ROI-pairs, a 90×90 adjacency matrix was attained for each subject.
2.4.2.1. Pearson correlation coefficient
A value for the linear association between the time series to be quantified is Pearson correlation, by which the dependence structure Gaussianity presupposition is given. Regarding {X}={x1,x2,xND} and {Y}={y1,y2,yND} as the voxel’s time series, the number of components in each set would be. You will get the Pearson correlation coefficient r as the following equation (
Gómez, Vaquero, López-Mendoza, González-Rosa, & Vázquez-Marrufo, 2004) in which is the mean of X and is the mean of Y.
(1)
2.4.2.2. Mutual Information (MI)
An effectual noise-robust correlation value and a nonrandom association between time series can be developed and constructed by mutual information (
Khan et al., 2007). To predict MI, we can predict the average number of bits of one of the two time series by expressing the value of the other. Therefore, quantifying the variable X, the average number of the variable Y (i.e. mutual information of the variables X and Y, represented by I(X, Y)) can be predicted by measuring variable X. As a result, we will have MI as
(Bonita et al., 2014;
Ward & Mazaheri, 2008):
(2)
, and correspondingly designate the number of bins in the variable X’s histogram and the number of components in the variable Y’s histogram, in which does not vitally equal , in general. Besides, the probability that a component in the ith bin of the X axis segregation and the jth bin of the Y axis segregation would respectively constitute an (x, y) pair is the joint probability distribution, i.e. PXY (i,j).
2.4.2.3. Box-Counting Fractal Dimension (BCFD)
In order to predict fractal dimension, we employ self-similarity being conceptualized as the fundamental principle. The equation below is given for the FD of a bounded set A:
(3)
In this regard, the minimum number of A’s discrete copies in the scale r is designated as Nr. We can compute FD just for fractals that are deterministic. Via a box-counting method, a D estimate (i.e. the box-counting DB) can be computed. In the following, we explain the DBC method (
Li, Du, & Sun, 2009) . In this case, by plotting each time series of fMRI voxel’s against the other, we are able to create an image of the size M×M (M equals the length of fMRI time series). This image is regarded as a 3D spatial surface, in which the position of pixel on the plane of the image and the gray level of pixel are denoted by (x,y) and (z), the third coordinate. The plane, in this method, is partitioned by the s×s blocks. In each block’s scale (i.e. r=s), s represents an integer and 1
(4)
(i, j)=l−k+1
Additionally, having taken all blocks’ contribution into consideration, is given for the varied measures by:
(5)
Ultimately, the minimum squares linear fit of log (Nr) versus log (l/r) would result in FD being prognosticated.
2.5. Surrogate data
By generating the surrogate data, we aim at evaluating the performance of each method in several pairs of time series, before their applying in the real data. That is because these methods are used as the similarity measure of two time series from two fMRI voxels. We applied the outlined ideas, comparing the total MI, LC and BCFD between the signals in surrogate datasets. These surrogates are generated using the logistic equation (Arora & Santhanam, 2014) and randomization process.
This approach provides us with the opportunity to both test and quantify the deviation from linearity, providing a principled guide in judging the appropriateness of LC, MI and BCFD as measures of FC.
(6)
x(n+1)=A.x(n).(1-x(n))
We exploited logistic equation and increased the A value randomly between 3