Graph analysis of resting-state ASL perfusion MRI data: Nonlinear correlations among CBF and network metrics

NeuroImage 87 (2014) 265–275

Contents lists available at ScienceDirect

NeuroImage journal homepage: www.elsevier.com/locate/ynimg

Graph analysis of resting-state ASL perfusion MRI data: Nonlinear correlations among CBF and network metrics Xiaoyun Liang a,⁎, Alan Connelly a,b, Fernando Calamante a,b a b

Florey Institute of Neuroscience and Mental Health, Heidelberg, Victoria, Australia Department of Medicine, Austin Health and Northern Health, University of Melbourne, Melbourne, Victoria, Australia

a r t i c l e

i n f o

Article history: Accepted 5 November 2013 Available online 16 November 2013 Keywords: Small-world networks Arterial spin labeling Functional connectivity Network metrics Connectome

a b s t r a c t Human connectome mapping is important to understand both normal brain function and disease-related dysfunction. Although blood-oxygen-level-dependent (BOLD) fMRI has been the most commonly used method for human connectome mapping, arterial spin labeling (ASL) is an fMRI technique to measure cerebral blood flow (CBF) directly and noninvasively, and thus provides a more direct quantitative correlate of neural activity. In this study, investigations on properties of CBF networks using ASL perfusion data have been conducted on 10 healthy subjects. As with BOLD fMRI studies, the extracted networks exhibited small-world network properties. In addition, highly connected brain regions are shown to overlap mostly with hub regions detected from BOLD fMRI studies. Taken together, this demonstrates the capability of ASL fMRI for mapping the brain connectome. Furthermore, a sigmoid model was then employed to fit the extracted network metrics vs. CBF measurements. Interestingly, the relationships between 4 specific network metrics and region-wise CBF demonstrate that consistently nonlinear patterns exist across all subjects. In contrast to the positive nonlinear pattern of other network metrics (degree, vulnerability, and eigenvector centrality), the characteristic path length shows a negative nonlinear pattern, reflecting the mechanism underlying the small-world properties. To our knowledge, this is the first study to unravel the intrinsic relationships between specific network metrics and CBF estimates. This should have diagnostic and therapeutic implications for those studies focusing on patients who suffer from abnormal functional connectivity. © 2013 Elsevier Inc. All rights reserved.

Introduction The properties of small-world networks have been shown to have intermediate topologies between regular and random networks (Watts and Strogatz, 1998), characterized by small characteristic path lengths, as in the random network, but highly clustered, as in the regular one. In their seminal paper, Watts and Strogatz demonstrated that models of dynamical systems with small-world properties show enhanced signal-propagation speed, computational power, and synchronizability. Among many properties of small-world networks, some are considered to be key, including degrees, characteristic path length, clustering coefficient, centrality, and vulnerability. These small-world network properties have been shown in a broad range of complex networks, including biological, social, economic, communication, and transport networks (Latora and Marchiori, 2001). Since then, small-world networks have attracted intense interest from researchers from different disciplines (Bullmore and Sporns, 2009; Milo et al., 2002; Newman et al., 2002).

⁎ Corresponding author at: Florey Institute of Neuroscience and Mental Health, Melbourne Brain Centre, 245 Burgundy St., Heidelberg, Victoria, 3084, Australia. Fax: + 61 3 9035 7307. E-mail address: [email protected] (X. Liang). 1053-8119/$ – see front matter © 2013 Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.neuroimage.2013.11.013

In particular, it has become widely accepted that the human brain is one of the biological small-world networks (Bassett et al., 2006; Sporns et al., 2004; Stam et al., 2007). All brain regions cooperatively work to process local information in a small subset of localized neighborhoods, but transfer and integrate the processed information globally (Latora and Marchiori, 2001). Furthermore, abundant evidence shows that some brain regions play central roles, as shown by their higher degrees (connections) than other brain regions, as well as shorter characteristic path lengths; these regions are therefore often termed “brain hubs”, in an analogy to hubs in computer networks (Hagmann et al., 2008; van den Heuvel and Sporns, 2011). Further understanding of such mechanism has been facilitated by the development of small-world network analyses. Thereafter, a good number of associated studies have been investigated, such as mapping structural networks using diffusion MRI (van den Heuvel and Sporns, 2011) or cortical thickness measurements (He et al., 2008), mapping functional connectivity using blood-oxygen-level-dependent (BOLD) fMRI measurements (Achard et al., 2006; Liang et al., 2012c; Salvador et al., 2005) or electrophysiological techniques (EEG, MEG or MEA) (Stam et al., 2006). Mapping structural networks can reveal the fundamental interregional structural connections, whereas mapping functional networks may add more neurophysiological dynamics. This is especially important

266

X. Liang et al. / NeuroImage 87 (2014) 265–275

for some clinical studies focusing on dysfunctional connections, in which there may not be easily detectable structural connectivity changes. In fact, differences between structural and functional connectivity have shown that while the repertoire of structural motifs is small, the number and diversity of functional motifs are maximized to optimize the brain networks (Sporns and Kotter, 2004). Among the existing techniques for mapping functional networks, resting-state BOLD fMRI is the most commonly used, due to its simple MRI acquisition sequence and relative high spatial resolution. However, given that BOLD signal changes reflect a combination of effects from blood oxygenation, cerebral blood volume, cerebral blood flow (CBF), and metabolic rate of oxygen (Buxton et al., 2004), their interpretation can be complex, as BOLD contrast provides an indirect measure of neural activity; furthermore, it may include a large contribution from draining veins (Boxerman et al., 1995). A further limitation of EPI-based BOLD sequences (the more commonly used in practice) is that signal-dropouts in regions with high susceptibility-gradients corrupt the measured BOLD signal. Arterial spin labeling (ASL) is a non-invasive MRI technique to measure CBF (Calamante et al., 1999; Detre et al., 1992), and thus provides a more direct quantitative correlate of neural activity in comparison to BOLD (Luh et al., 2000). This typically comes at the expense of a reduced spatial and temporal resolution. However, previous results have demonstrated that ASL-based resting-state functional connectivity can successfully extract equivalent networks to BOLD fMRI data (Liang et al., 2012b). Furthermore, since ASL is now often used with spin-echo based methods, it can achieve better sensitivity in brain regions with high-susceptibility gradients (Liang et al., 2012b). Importantly, mapping resting-state networks using ASL can not only detect the spatial extent of the networks but also quantify CBF, which can be used to infer metabolic energy consumptions and thus provide a more biologically specific correlate of neural activity compared to BOLD fMRI (Liu and Brown, 2007). Measurements of metabolism and blood flow during resting-state in an early PET study found that certain brain regions consume more metabolic energy during resting-state than task-states; these regions were then referred to as the “default mode network” (DMN) (Raichle et al., 2001). Thereafter, brain regions within DMN, along with some other regions (such as temporal cortex, frontal cortex, and occipital cortex) have been intensively investigated and reliably classified as hub regions (Achard et al., 2006; Liang et al., 2012c). Interestingly, a recent study showed that significantly decreased CBF within some brain hubs was detected in subjects during a psychedelic state, which caused ‘unconstrained cognition’ (Carhart-Harris et al., 2012). It was then suggested that the effects could be elicited by decreased activity and connectivity between hub regions and other brain regions. This is in contrast with the highly connected brain hubs from normal subjects during resting-state, which reflects constrained cognition to maintain the self-referencing state. In fact, to maintain such a resting state, the energy consumption is almost identical to that of a task-related state (Fox and Raichle, 2007). The large costs of spontaneous activity in terms of metabolic energy consumptions and highly organized inter-region connections may suggest that spontaneous activity is crucial to maintain the normal self-referencing state. No consensus on the specific role of the spontaneous activity has been reached yet, although two possibilities have been considered. One is that spontaneous activity keeps track of previous memory, showing correlations between regions that have been modulated together in a taskdependent manner. The other is that spontaneous activity organizes and coordinates neuronal activity such that information can be transferred and integrated between regions in concert more efficiently. Whatever the role is, it is important to determine how the metabolic energy consumption is modulated to maintain a constrained cognition, characterized by a highly efficient and cooperative complex network. Given the above-mentioned interesting findings of resting-state networks and the potential advantages of ASL over BOLD fMRI, CBF

studies that focus on the small-world properties of resting-state networks would allow us to investigate the relationships between network metrics (such as degrees, and characteristic path lengths) and CBF more directly and interpretably. However, to our knowledge, such studies focusing on small-world networks of human brain using ASL have not been reported. In this study, we investigate the mapping of resting-state CBF networks using ASL fMRI, which allows us to characterize the CBF network metrics of human brain during resting state based on direct blood flow measurements, as well as to investigate the relationship between CBF and network metrics. This should provide useful information to further our understanding of metabolic energy consumption, with which the brain can maintain normal cognition with the lowest cost. More specifically, we hypothesized that regional CBF is positively correlated with the degree, centrality and vulnerability, but negatively correlated with characteristic path length. Materials and methods Subjects Ten normal subjects were recruited to participate in this study and data were acquired on a 3T Siemens Tim Trio scanner with a 12channel head coil (Erlangen, Germany). All subjects provided written, informed consent, and all protocols were approved by the local Institutional Review Board. Data acquisition According to previous findings (Gonzalez-At et al., 2000), and our previous experience on resting-state functional connectivity using ASL (Liang et al., 2012b), a short post-labeling delay (PLD) ASL sequence is preferred to a long PLD sequence due to the relatively higher SNR of the former. However, short PLD is known to have problems for CBF quantification (Wang et al., 2003). Therefore, ASL datasets with both short PLD and long PLD were collected for the purposes of connectivity and CBF quantification estimation, respectively. Short PLD 3D GRASE pCASL sequence with k-space sharing to improve brain coverage (Liang et al., 2012b), TR/TE = 3750/56 ms, resolution = 4 × 4 × 6 mm3, 20 axial slices, matrix size = 64 × 51, PLD = 600 ms. Labeling position was determined according to (Aslan et al., 2010), with labeling duration = 1284 ms. Background suppression (BS) was achieved using inversion-times of 1913 ms and 523 ms (Garcia et al., 2005). For each subject, 120 images (60 label/control pairs) were acquired. Long PLD Labeling duration = 1584 ms, PLD = 1540 ms, BS inversion-times: 1800 and 520 ms (Fernandez-Seara et al., 2008), otherwise identical parameters to those of the short PLD sequence. Anatomical scan and ASL reference Structural images were acquired for each subject using a Magnetization-Prepared Rapid Gradient Echo (MPRAGE) sequence (TR/TE/TI = 1900/2.55/900 ms, flip-angle = 9°, 1 mm isotropic resolution, field-of-view = 160 × 214 mm2, 160 partitions). For the purpose of CBF quantification as well as accurate registration, a reference scan with higher intensity (M0) was performed with a long TR (8000 ms) with the same parameters as those used for ASL, except that no labeling or BS was included. Image analysis Preprocessing Image processing was performed by using SPM8. For each subject, images (with long and short PLD) were preprocessed as follows: (i) all

X. Liang et al. / NeuroImage 87 (2014) 265–275

ASL images were motion-corrected, spatially smoothed using an isotropic Gaussian kernel with 6 mm FWHM, followed by co-registration with the individual reference M0 image to the high-resolution anatomical image; (ii) data were then normalized to MNI standard space (resolution: 3 × 3 × 3, matrix size: 61 × 73 × 61); (iii) perfusion images were obtained by pairwise subtraction; (iv) perfusion images were detrended and band-filtered (0.01Hz b f b 0.1 Hz) (Fox and Raichle, 2007); (v) the posterior probability maps of gray matter (GM) and white matter (WM) were segmented from the individual anatomical image and coregistered to the MNI standard space (resolution: 3 × 3 × 3) using the MNI template as a reference. Independent component analysis (ICA) The aim of performing ICA in this study is to extract intravascular ‘artifacts’ since these are present in ASL data with short PLD (Gonzalez-At et al., 2000), and they could therefore affect the calculated network metrics. This was conducted at an individual-level as implemented by FSL MELODIC software. Based on our previous experience (Liang et al., 2012b), about 2 or 3 components were easily identified as intravascular-related components, which had high intensity on the spatial components distributed around arteries. These components were then removed by using the MELODIC denoising process, therefore minimizing the intravascular effects on our short PLD ASL data. Adjacency matrix construction Network nodes were defined by using the Anatomical Automatic Labeling (AAL) atlas (Tzourio-Mazoyer et al., 2002). From the commonly used 90 regions (i.e. after excluding the cerebellum Achard et al., 2006; Liang et al., 2012c), we considered 84 regions, after further excluding bilateral supplementary motor area, paracentral lobule and olfactory region (Table 1). Our exclusion criterion is based on the empirical observation of wrap-around artifacts (Yang et al., 2010) on the superior-most slices of the brain in our 3D GRASE acquisition. These artifacts are

Table 1 Eighty-four cortical and subcortical regions defined in our study from Automated Anatomical Labeling template. Index

Region

(1,2) (3,4)

Index

Precentral gyrus Superior frontal gyrus, dorsolateral (5,6) Superior frontal gyrus, orbital part (7,8) Middle frontal gyrus (9,10) Middle frontal gyrus, orbital part (11,12) Inferior frontal gyrus, opercular part (13,14) Inferior frontal gyrus, triangular part (15,16) Inferior frontal gyrus, orbital part

(51,52) (33,54) (55,56) (57,58) (59,60) (61,62)

(17,18) (23,24) (23,20) (27,27) (29,30)

Rolandic operculum Superior frontal gyrus, medial Superior frontal gyrus, medial orbital Gyrus recrus Insula

(63,64) (65,66) (67,08) (71,72) (73,74)

(31,32) Antezioi cmgulatcdandparacinzuLite gyri (33,34) Median cmgulate acdparacnrolate r.n (35,36) Postenor cingulated gyrus (37,38) Hippocampus (39,40) Parahippocampal gyrus

(75,76)

(41,42) Amygdala (43,44) Calcarine fissure and surrounding cortex (45,46) Cuneus

(85,86) (87,88)

Region

(47,48) Lingual gyrus (49,50) Superior occipital gyrus

(77,78) (79,80) (81,82) (83,84)

(89,90)

Middle occipital gyrus Inferior occipital gyrus Fusiform gyrus Postcentral gyrus Superior parietal gyrus Inferior parietal, but supramarginal and angular gyri Supramarginal gyrus Angular gyrus Precuteus Caudate nucleus Lenticular nucleus, putuucn Lenticular nucleus. pelLcum Thalamus Heschi gyrus Superior temporal gyrus Temporal pole superior temporal gyrus Middle temporal gyrus Temporal pole middle temporal gyrus Inferior temporal gyrus

267

inherent to the 3D imaging sequence, which corrupts the ASL signal changes on the superior slices, and affect four of the above-mentioned brain regions. With respect to the olfactory region, since it is located close to major arteries, it is often subject to severe intravascular artifacts. In fact, this region was observed in our data with abnormally high CBF values even with the long PLD sequence. This region was therefore excluded from our analysis. Region-wise mean time-courses were extracted across all voxels within each of the 84 brain regions for each subject. The region-wise mean time-courses were employed to estimate functional connections among regions (based on their Pearson's correlation coefficients). This resulted in one (84 × 84) correlation matrix for each subject representing interregional functional connectivity. To generate group-level connectivity results, network metrics were calculated at the individual-level (using the same threshold), and subsequently averaged to yield mean network metrics at group-level. Since the interregional correlations (84*83/ 2 = 3486) were multiple non-independent tests, a stringent threshold was employed by applying Bonferroni correction (p = 0.001/ 3486 = 2.87 × 10−7) according to the literature (Liao et al., 2010), which yielded r = 0.6 as the correlation threshold. However, to assess the effect of choice of correlation threshold, we also repeated the analysis for a range of threshold values (0.4–0.8). In this study, only the binary correlation matrix was considered; any voxel is assigned 1 if its correlation coefficient is larger than the threshold, and 0 otherwise (Note: in graph theory, this is usually called ‘adjacency matrix’ (Stam and Reijneveld, 2007)). Small-world network analysis Small-world network analysis was conducted on the individuallevel adjacency matrix. Network metrics were analyzed based on graph theory, with a graph described by a number of nodes and undirected connections between each pair of nodes if applicable (Stam and Reijneveld, 2007). Key metrics for characterizing small-world networks were calculated for each node, such as degrees, clustering coefficient, characteristic path lengths, vulnerability and centrality (see below). The small-worldness σ was assessed for the networks detected using ASL by comparing the characteristic path lengths (L) and clustering coefficient (C) with the same metrics measured in 1000 matched random networks (the same number of nodes, edges, and the same degree distribution as the real network) (Maslov and Sneppen, 2002; Sporns and Kotter, 2004; Sporns and Zwi, 2004). To express the small-worldness of a network, a small-world index is calculated σ ¼ γ=λ

ð1Þ

where γ = Cnet/Cran and λ = Lnet/Lran, and the subscripts ‘net’ refers to the network under consideration and ‘ran’ to random networks. Since it is expected that γ N 1 and λ ~ 1 for any network that possesses the small-world properties (Watts and Strogatz, 1998), σ is therefore larger than 1 for any small-world network. Considering an N × N adjacency matrix (N = 84), the aforementioned key metrics are defined as follows: Degree the degree Ki of a node i is defined as the number of connections to that node. Clustering coefficient The clustering coefficient Ci of a node i is defined as the number of connections of all its neighbors divided by all possible connections, which can be calculated as: Ci ¼

CEi K i ðK i −1Þ=2

ð2Þ

where CE i represents the number of existing connections (Watts and Strogatz, 1998). The clustering coefficient of a network, Cnet is the average of the C i in all nodes i in the network.

268

X. Liang et al. / NeuroImage 87 (2014) 265–275

Characteristic path lengths Characteristic path length of a node i can be calculated as: Lpi ¼

1 X L N−1 i≠ j∈G i; j

ð3Þ

where G is a graph and Li,j is the minimum number of edges that needs to be traversed to make a connection between nodes i and j. The characteristic path length of a network, Lnet, is the average of the characteristic path length Lpi in all nodes i in the network. Vulnerability The vulnerability Vi of a node i can be defined as the change in global efficiency of the network after the node i and all its edges are removed (Shu et al., 2009):   ′ V i ¼ Eglob −Eglob =Eglob

ð4Þ ′

where Eglob is the global efficiency of the network (and Eglob its value after removing node i): Eglob ¼

X 1 1 NðN−1Þ i≠ j∈G Li; j

ð5Þ

Eigenvector centrality It is a measure of the influence of a node in a network (Bonacich, 1972), hereafter, denoted as eigencentrality. For a given graph G, the eigencentrality of node i can be defined as: EC i ¼

1 X 1 X EC ¼ a EC λd n∈MðiÞ n λd n∈G i;n n

ð6Þ

where M(i) is a set of the neighbors of node i and n represents a node in M(i), and A = (ai,n) is the adjacency matrix, which can be written as eigenvector equation, Ax = λd x. Here λd is an eigenvalue of the adjacency matrix A (and, therefore, there can be different values; however, only the eigenvector corresponding to the largest eigenvalue results in the desired centrality measure (Newman et al., 2002)). Network efficiency Network efficiency can be classified in terms of global (see Eq. (5)) and local efficiency, i.e. characterizing its efficiency of global and local information processing and transfer. The global efficiency is given by the expression in Eq. (5); the corresponding expression for the local efficiency is: Eloc

1X ¼ E ðG Þ; N i∈G glob i

ð7Þ

Where Eglob (Gi) is the global efficiency of Gi, the subgraph of the neighbors of node i (Latora and Marchiori, 2001). Degree distribution fitting Three possible forms of distribution P(k) were employed to fit the degree distribution: the power law, p(k) ~ kα-1; the exponential, p(k) ~ e-αk and the exponentially truncated power-law, p(k) ~ kα-1ek/kc (Achard et al., 2006). To reduce the effects of noise on this relatively small dataset, cumulative distribution was employed for evaluating the goodness-of-fit with R2 values. CBF quantification corrections As a general step, absolute quantification of CBF was calculated according to the ref. (Wang et al., 2008). While the pCASL 3D GRASE sequence has been shown to have several advantages as an ASL sequence (Fernandez-Seara et al., 2007), there are a number of artifacts that can have a major influence on CBF quantification. Since the accuracy of CBF measurement can largely affect any subsequent analysis of

correlation of small-world metrics with CBF, these artifacts need to be minimized or corrected. In particular, CBF signal loss caused by imperfect slab profile in 3D imaging (Fernandez-Seara et al., 2007) and CBF affected by partial volume effects (Liang et al., 2013a) were corrected as follows: (1) Slab-profile effect (SPE) correction To correct for SPE, we used the following procedure. M0 images were employed to estimate the signal loss, based on the assumption that the mean intensity across the voxels within each slice is approximately constant along the z direction in an M0 image (excluding regions with very high intensity values, such as CSF). Slices located in the middle of the field-of-view (FOV) were assumed to have no signal loss, and they yield maximum mean intensity, which was divided by the corresponding value from the 20 slices, thus generating one correction coefficient for each slice. The yielded coefficients were then applied to improve the estimate of CBF with the long PLD sequence. To assess the proposed correction method, an extra acquisition with a shifted FOV was acquired in 2 subjects (FOV shifted downwards 30 mm, see Supplementary figure Fig. S1 (B)), such that the inferior slices in the original FOV acquisition (which were subject to severe SPE) are now located towards the middle of the FOV (and therefore with no or minimal SPE); see Fig. S1 for more details. By comparing the CBF measurements from these FOV acquisitions, the performance of the SPE correction method can be evaluated. (2) Partial volume correction A recently proposed improved partial volume effect (PVE) correction method using modified least trimmed squares (mLTS) was then employed to correct for PVE for each subject (Liang et al., 2013a). The mean CBF map after SPE correction along with the posterior probability maps of GM and WM in standard space were employed as inputs for this PVE correction method, thus generating a PVE corrected CBF map.

Region-wise mean CBF extraction For each subject, the mean CBF for each of the 84 brain regions defined above were extracted from long PLD CBF map (in MNI standard space, resolution: 3 × 3 × 3) after both SPE and PVE corrections. To assess the effect of the various correction steps, we also calculated the mean CBF values from long PLD CBF data with SPE correction only (i.e. without PVE correction), and without SPE correction. Statistical analysis To determine the relationship between CBF and ASL-based network metrics (degree, characteristic path length, vulnerability and centrality), a sigmoid function was employed to fit a predictive model to observed region-wise network metrics vs. CBF. Sigmoid functions have the following advantages: (1) they are bounded above and below; (2) by scaling the input to a sigmoid function, a range of slopes can be obtained; (3) by adding constants to the input, the ‘position’ of the curve can be moved; it is, therefore, suitable for fitting data lacking detailed description, such as ours (Gibbs and MacKay, 2000). (Note: a simpler linear model was investigated in an initial exploratory analysis of the data, but it was found to be unsuitable, based on the distribution of the residuals, which were found to be systematically primarily on one side of the fitted model for high and low CBF values – data not shown). The nonlinear fitting analyses were performed at group-level in Matlab using the routine NonLinearModel class. More specifically, CBF values were extracted for each subject and mean CBF was obtained by averaging CBF across the 10 subjects; network metrics were calculated for individual subject first, which were subsequently averaged to yield mean network metrics at group-level. Since three CBF values were extracted (see previous section), each analysis was conducted three times for each subject. To assess the significance of the nonlinear fitting

X. Liang et al. / NeuroImage 87 (2014) 265–275

results, a model-based resampling technique called residual bootstrap was employed (Wu, 1986) (see Appendix A for details). Mapping of network metrics and regional CBF For the purpose of illustrating the spatial relationship between network metrics and CBF, both the region-wise connections and the reciprocal characteristic path length of nodes, as well as region-wise mean CBF measurements, were mapped onto the cortical surfaces of the MNI standard template using the BrainNet Viewer software (www. nitrc.org/projects/bnv). The reciprocal of characteristic path length is employed, such that a shorter characteristic path length (e.g. as found in ‘hubs’) is represented by a larger circle.

269

properties across the whole range (i.e. sigma N 1). While both sigma and characteristic path length increase as r increases, the other 3 metrics (clustering coefficient, local and global efficiency) decrease as r increases (Figs. 2(C)–(E)).

Results Corrections of SPE for CBF Fig. 1(A) and (B) clearly demonstrates that the CBF decrease at the edge of the slab for the original unshifted FOV (blue line) was compensated with the method for SPE correction (green); this compensation was further validated by its agreement with the true slice profile, as calculated from the shifted FOD case (red). An example of effect of the SPE correction on the CBF maps can be found in supplementary materials (see Fig. S2). Small-world networks In this study, small-world properties (sigma, characteristic path length, clustering coefficient, local and global efficiency) of brain networks using ASL fMRI were calculated over a range of correlation threshold values (r = 0.4–0.8 with step size of 0.05) (see Fig. 2). Fig. 2(A) demonstrates that the extracted networks possess small-world

Fig. 1. Mean CBF measurements (averaged across each slice) for 2 typical subjects are shown in (A) and (B): before SPE correction (blue), after SPE correction (green) and with FOV shift downwards (30 mm relative to original position – see Fig. S1(B) in supplementary material) to reduce SPE in the inferior part of the original FOV (red) (see Fig. S1 (A)). By shifting the FOV, CBF measurements of slices 6th–15th (along inferior–superior direction) are considered as true measurements of mean CBF (i.e. without SPE). For simpler visual comparison, please note that the red curve corresponds to slices 6th–15th (from the shifted FOV case), which have been shifted along the x-axis in the graph by the appropriate amount (corresponding to the 30 mm FOV shift).

Fig. 2. Small-world network properties as a function of correlation threshold, r. (A) The small-world index sigma increases as r increases; (B) the characteristic path length increases as r increases; the clustering coefficient (C), local (D) and (E) global efficiency decrease as r increases.

270

X. Liang et al. / NeuroImage 87 (2014) 265–275

The relationship between each of the 4 network metrics and CBF at individual-level consistently demonstrated the nonlinear trends across the range of thresholds (data not shown – Note: group results are presented below in the “Network metrics vs. CBF” Results section). Fig. 3 shows the mean fitting values of R2 (calculated at individual level, with the error bars representing standard deviations across subjects). While the performance of fitting (in terms of mean R2) is generally improved with increasing thresholds, the trend reverses for very large thresholds, as shown by the increased error bars and/or decreased R2 (e.g. Fig. 3 (B)). Based on considering all of the data in Fig. 3, the threshold r = 0.6 was deemed optimal for the subsequent data analysis. For this threshold, the group-level small-worldness index is σ = 1.53 (see Fig. 2A), which is obtained by averaging the values of sigma across the group of 10 subjects. By fitting degree distribution with the three possible distributions described in the method section, R2 values were 0.65, 0.78 and 0.91 for the power law, the exponential, and the exponentially truncated power-law distribution, respectively. The exponentially truncated power-law is, therefore, the best-fit model, which is consistent with previous studies on both functional networks using BOLD fMRI (Achard et al., 2006; Hayasaka and Laurienti, 2010) and anatomical networks (He et al., 2008). Hub regions While brain regions with degrees larger than a normal threshold (mean + standard deviation) are usually considered as hub regions (van den Heuvel and Sporns, 2011), a less stringent threshold, mean + 1/2 standard deviation, was also chosen in order to not exclude regions that may be hub regions but which have a large variance due to relatively low SNR of ASL. Table 2 lists 39 brain regions as likely hub regions, with regions in bold representing hubs determined with the standard threshold (i.e. mean + standard deviation), which are consistent with brain connectivity studies using BOLD fMRI (Achard et al., 2006). As expected, all these possible hub regions are characterized by shorter characteristic path length, but higher vulnerability and eigencentrality. Network metrics vs. CBF Overall, the 4 network metrics vs. region-wise CBF measurements show nonlinear trends, of which, the 3 network metrics Ki, Vi and ECi follow positive nonlinear trends; Lpi, however, shows a negative trend. Fig. 4 shows the nonlinear fitting outcomes of 4 network metrics vs. CBF at group-level. The group-level R2 across the 10 subjects are 0.36, 0.31, 0.33 and 0.37 for nonlinear fitting outcomes of 4 network metrics, Ki, Lpi, Vi and ECi, vs. CBF with SPE correction, respectively. While a specific threshold was used for that analysis (r = 0.6, based on the results of Small-world network analysis Section above – see also Fig. 3), our data consistently showed nonlinear relationships between each of the 4 metrics and CBF across the range of thresholds (R = 0.4–0.8) (see Supplementary Fig. S3). As expected, more favorable fitting outcomes of Vi vs. region-wise CBF estimates are shown with than without SPE corrections (see Inline Supplementary Table S1); a pair-wise t-test showed that the fitting was significantly improved by using SPE correction (p = 5.3 × 10−14). In contrast, partial volume correction of CBF was not found to have a significant effect to improve the fitting (data not shown), which suggests that averaging the CBF across the ROIs is likely to minimize the benefits of partial volume correction for CBF. Therefore, only fittings of network metrics vs. CBF from data with SPE corrections are presented in the remainder of this work. Inline Supplementary Table S1 can be found online at http://dx.doi. org/10.1016/j.neuroimage.2013.11.013. Fig. 5 shows the results of the bootstrap analysis, where the R2 values are plotted for the nonlinear fittings (1000 resampled datasets using the residual bootstrap method) for each of the network metrics (Ki, Lpi, Vi

Fig. 3. Plots of the mean values of R2 (nonlinear fittings) of the 4 network metrics calculated from individual subject, with the error bars representing standard deviations across subjects. It can be appreciated that threshold r = 0.6 yielded relatively lower standard deviations and more robustness for all the 4 metrics than the rest of the thresholds.

and ECi) vs CBF. Furthermore, their 95% confidence intervals are also included, with green dashed-line representing upper confidence bound (Mean + 1.96*Standard deviation) and red dashed-line representing lower confidence bound (Mean-1.96*Standard deviation). These plots clearly show that R2 values are randomly distributed around their mean values, corresponding to 0.37, 0.33, 0.33 and 0.38 for Ki, Lpi, Vi and ECi vs CBF, respectively. Their confidence intervals are [0.21, 0.54], [0.15, 0.52], [0.17, 0.5] and [0.23, 0.54], respectively. As expected, most R2 values are within their associated confidence intervals. Furthermore, nonlinear relationship between each of the 4 network metrics was shown to be significant, with largest p values = 0.002, 0.03, 0.03 and 0.006 for Ki, Lpi, Vi and ECi vs CBF, respectively (see Supplementary Fig. S4 for p-values plots). These results obtained with the residual bootstrapping method, therefore, confirmed the significance of the

X. Liang et al. / NeuroImage 87 (2014) 265–275

271

Table 2 Hub regions of CBF networks ranked in order of decreasing degree (Kf). Brain regions in bold represent those with degrees larger than a normal threshold (mean + STD); the rest of the brain regions possess degrees larger than a less stringent threshold (mean + 1/2STD). Abbreviation

Regions

Ki

Lpi

Vi

ECi

CBF

MTG PCUN PCUN MTG CUN IFGorb SOG ITG SMG PoC'G CUN IFGorb STG ANG IFGniang SMG CAL ANG PoCG LING SFGmed SFG SFGmed IFGoperc MOG CAL PreCG FG MCC LING MFG FG IFGoperc IPL STG MOG ORBsup PreCG IFGtriang

Middle temporal gyrus (R) Precuneus (R) Precuneus (L) Middle temporal gyrus (L) Cuneus (R) Inferior frontal gyrus-orbital (L) Superior occipital gyrus (L) Inferior temporal gyrus (R) Suprainarginal gyrus (R) Postcentral gyrus (R) Cuneus (R) Inferior frontal gyrus-orbital (R) Superior temporal gyrus (R) Angular gyrus (R) Inferior frontal gyrus-Triangular (L) Supramarginal gyrus (L) Calcarine cortex (L) Angular gyrus (L) Postcentral gyrus (L) Lingual gyms (L) Superior frontal gyms-medial (R) Superior frontal gyrus (L) Superior frontal gyms-medial (L) Inferior frontal gyrus-opercular (L) Superior frontal gyrus (L) Calcarine cortex (R) Inferior temporal gyrus (L) Fusiform gyms (R) Middle cingulated cortex (R) Middle cingulated cortex (L) Middle Frontal gyrus (L) Fusiform gyms (R) Inferior frontal gyrus-opercular (R) Inferior parietal lobule (L) Superior temporal gyrus (L) Middle occipital gyrus (R) Superior frontal gyrus-orbital (R) Precentral gyms (R) Inferior frontal gyrus-triangular (R)

59 55 55 54 54 54 53 53 52 52 52 52 52 52 52 50 50 50 50 50 49 49 49 49 49 49 4S 4S 4S 4S 4S 47 47 47 47 47 47 47 47

1.37 1.36 1.38 1.36 1.36 1.38 1.39 1.40 1.40 1.41 1.40 1.40 1.40 1.42 1.43 1.43 1.42 1.45 1.45 1.43 1.44 1.45 1.43 1.44 1.44 1.43 1.4S 1.44 1.46 1.45 1.47 1.49 1.45 1.49 1.47 1.48 1.49 1.51 1.47

0.0048 0.0032 0.0034 0.0029 0.0030 0.0037 0.0028 0.0026 0.0024 0.0023 0.0024 0.0027 0.0024 0.0023 0.0022 0.0021 0.0021 0.00 IS 0.0019 0.0020 0.0020 0.0018 0.0022 0.001S 0.0019 0.0017 0.0014 0.0017 0.0016 0.0015 0.0039 0.0014 0.0016 0.0012 0.0014 0.0015 0.0013 0.0012 0.0012

0.1477 0.1383 0.1406 0.1362 0.1350 0.1263 0.1366 0.1350 0.1319 0.1356 0.1282 0.1159 0.1142 0.1299 0.1251 0.1237 0.1236 0.1240 0.1291 0.1203 0.1160 0.1207 0.1174 0.1153 0.1264 0.1185 0.1204 0.1152 0.1195 0.1157 0.1163 0.1141 0.1071 0.1187 0.1023 0.1212 0.1138 0.1184 0.1065

84 80 76 67 81 60 80 67 83 72 84 57 79 94 64 6S 74 82 64 63 63 55 5S 66 S4 70 63 51 61 64 61 4S 64 72 64 90 4S 64 57

nonlinear relationships between the 4 network metrics and CBF measurements. Fig. 6 visualizes the intrinsic relationships among region-wise connections, the characteristic path length of the nodes and their CBF values, by mapping all 3 measurements onto the standard MNI template at group-level. Each line linking 2 nodes represents a realistic connection at a given threshold (r = 0.75. Note: this threshold is higher than that used for deriving network metrics (i.e. 0.6, see Adjacency matrix construction Section) for the sole purpose of improved visualization of connections by controlling network sparsity). As expected, hub regions are characterized by denser connections, shorter path length and higher CBF values (see Fig. 6). This figure illustrates one of the advantages of ASL compared with BOLD data, given that a physiological parameter, namely CBF, can be assigned to each of the nodes, thus enhancing the available information in the connectome. Discussion In this study, we investigated the small-world properties of CBF networks using ASL fMRI data. ASL data with both short PLD (for mapping brain connectome) and long PLD (for CBF quantification) were acquired, which allowed us for the first time to measure network metrics of human brain during resting-state based on ASL data, as well as to investigate the intrinsic relationships among CBF network metrics with CBF measurements. Consistent with previous studies using BOLD fMRI, the small-world index of CBF network is larger than 1 (~ 1.53 in our study), which is

considered an important characteristic of small-world networks (Achard et al., 2006; Watts and Strogatz, 1998). The small-world properties of the CBF networks have also been demonstrated by the fact that the exponentially truncated power law is the best-fit model for degree distribution (Achard et al., 2006; He et al., 2008). Such small-world properties of CBF networks have been shown recently by using SPECT imaging (Melie-Garcia et al., 2012). However, that study adopted a very different approach from the present work: they investigated the concurrent changes between CBF measurements in brain structures across subjects and, importantly, no temporal dynamics of brain networks was considered due to the lack of temporal information in their SPECT data. The capability of ASL fMRI for mapping the brain connectome has been further corroborated by the evidence that the brain regions with top ranks in degree are mostly overlapped with those brain regions that are considered as hub regions using BOLD fMRI (Achard et al., 2006). However, the correspondence is not one-to-one, suggesting that differences in the underlying mechanisms of ASL and BOLD fMRI might give rise to some differences in the network characteristic identified. Given that there are some potential advantages of ASL relative to BOLD fMRI (Liu and Brown, 2007; Luh et al., 2000), the perfusionbased methodology has attracted increasing interest in neuroimaging studies. However, its applications in functional connectivity have not been widely used to date; this is likely due to a number of factors, including the lower SNR of ASL signal, reduced temporal resolution, and difficulties with whole brain coverage. Nevertheless, it has been recently demonstrated that a 3D GRASE pCASL sequence with k-space sharing is

272

X. Liang et al. / NeuroImage 87 (2014) 265–275

Fig. 4. Nonlinear fitting of 4 network metrics (Vi, Lpi, Ki and ECi) vs. CBF across a group of 10 subjects.

capable of mapping whole-brain functional connectivity (Liang et al., 2012b). Using the same technique in the current study, we have shown that the relationships among network metrics and CBF follow nonlinear patterns at the group-level. Notably, the inverse relationship between Lpi and CBF compared to the other 3 metrics (Ki, Vi and ECi) vs. CBF reflects the mechanism underlying the small-world properties of the brain. This finding is in line with previous findings that ‘hub’ regions tend to have higher values for degree, vulnerability, and centrality, but lower values for characteristic path length (Achard et al., 2006; Gong et al., 2009). Moreover, regardless of the lack of investigations on explicit relationships among network metrics and CBF measurements, our findings of brain regions with higher degrees demanding heavier CBF supply are also consistent with the previously reported higher metabolic energy consumption (and CBF) in hub regions (Raichle et al., 2001). This, therefore, has diagnostic and therapeutic implications when applied to those patients who suffer from abnormal functional connectivity. This study integrates the physiological measurements of both CBF network metrics and CBF measurements, and thus provides more direct biological information than can be provided by functional connectivity using the more indirect BOLD signal. In fact, this advantage can be readily appreciated from Fig. 6, where ASL can be used not only to determine the connectivity (and therefore the properties of the network), but also to assign an important physiological parameter to each node. Our findings are also consistent with recent studies that showed that the metabolic costs of a node are proportional to its degree and/or its centrality, by fitting data of centrality rank and the glycolytic index in 41 Brodmann areas of the cerebral cortex (Bullmore and Sporns, 2012; Vaishnavi, 2010; Vlassenko et al., 2010). While a linear relationship was investigated in a recent study, there is no evidence from the data plots (Bullmore and Sporns, 2012) that a linear model would be more suitable than a nonlinear model. In fact, the distribution of residuals in the data plots indicates that a nonlinear model might describe the relationship between centrality (and/or degree) and the metabolic costs better. Interestingly, a very recent study reported a linear relationship

between regional CBF and functional connectivity strength detected with BOLD fMRI (Liang et al., 2013b). Once again, however, the residuals in the plots suggest that nonlinear models might provide a better fit to the data. We, therefore, believe that these findings would be consistent with ours on the nonlinear relationships among the 4 network metrics and CBF measurements. In fact, the observed nonlinear patters could be expected based on biological grounds. It has been widely accepted that small-world networks are self-organizing (Rubinov et al., 2011), and self-organization is an important component for a successful ability to establish networking whenever needed. Moreover, the dynamic of a self-organizing system is typically characterized by a nonlinear system, because of circular or feed-back relations between components (Heylighen, 2001). Therefore, these previous findings provide strong support for the detection of nonlinear patterns of 4 specific network metrics vs. CBF, revealing an intrinsically nonlinear mechanism underlying the CBF network.

Methodological issues and future work There are a few methodological considerations that may be relevant for this study. One main factor that would affect mapping functional connectivity is the chosen threshold for the Pearson's correlation coefficients. In this study, the threshold (r = 0.6) was reasonably chosen with the consideration of Bonferroni correction so as to ensure that the false positive rate is b 0.001 (Note that this choice of threshold was also considered appropriate based on the results from Fig. 3). However, choice of the optimal threshold remains somewhat ‘arbitrary’ (Wang et al., 2009). Varying the threshold may affect the small-world network properties to a certain degree. Furthermore, the same threshold may not be optimal for each subject in the group. Therefore, applying the same threshold to all subjects may unfavorably affect the small-world properties at group level. Therefore, further work needs to be carried out focusing on developing more robust threshold free methods.

X. Liang et al. / NeuroImage 87 (2014) 265–275

273

intravascular artifacts could still be present, which would affect network metrics; this could potentially explain the presence of some ‘outlier’ data points in the results shown in Fig. 4. The development of improved head coils (e.g. with 32 channels or more) can greatly enhance SNR of ASL data, opening up the possibility of performing reliable ASL functional connectivity studies based only on long PLD ASL data (i.e. for the purposes of both CBF quantification and connectivity analysis). Given that ASL has intrinsically low SNR (Calamante et al., 1999), the reliability of extracting networks and thus that of measuring network metrics may be compromised. While our results show that it is less an issue for studies focusing on region-based functional connectivity, it may be problematic for voxel-wise functional connectivity studies using ASL data, especially at individual level. If that were the case, denoising methods that can enhance the SNR of the ASL image (e.g. (Wells et al., 2010)) may play an important role. Lastly, we note that while, as a conceptual study, a more common method (namely Pearson cross correlation) was employed to detect functional connectivity, partial correlation estimation has advantages over this approach. As a further step, advanced methods combined with sparse solutions will be developed and applied to improve the network analysis. Conclusion In this study, ASL fMRI has been employed to investigate the smallworld properties of CBF networks within the normal human brain, whereby intrinsic relationships among 4 common network metrics and region-wise CBF measurements have been sought. Our study shows that Ki, Vi and ECi show positive nonlinear relationships to regional CBF, but a negative trend for Lpi. To our knowledge, this is the first study that unravels the intrinsic relationships between specific network metrics and CBF estimates. This should provide useful information to further our understanding of the metabolic energy consumptions with which the brain can maintain normal cognition with the lowest cost. It should also have diagnostic and therapeutic implications for those studies focusing on patients who suffer from abnormal functional connectivity. Acknowledgments

Fig. 5. The 95% confidence interval for R2 values for the nonlinear fitting of Ki, Lpi, Vi and ECi vs CBF using the sigmoid model based on a bootstrap analysis (see text for more details). Here CB denotes confidence bound, with red and green dashed lines representing lower and upper confidence bounds, respectively.

Another factor to be considered in this study is that interregional functional connectivity is investigated with brain regions defined in the AAL template (Achard et al., 2006; He et al., 2008), which reduces the data from a voxel resolution to a regional resolution. This is usually for the purpose of reducing the computational burden based on the simplistic assumption that all voxels within each region possess the same connectivity properties. Realistically, this may not be the case and there may be intra-regional connectivity between sub-parts of brain regions (Hayasaka and Laurienti, 2010; van den Heuvel et al., 2008). This would affect the detection of brain networks by simply averaging the time courses across all voxels in each region. A voxel-wise approach that is able to localize the ‘hub’ at voxel-level will be explored in our future work. Despite the benefits of short PLD ASL for functional studies (Gonzalez-At et al., 2000), it inevitably introduces intravascular artifacts. To address the issue, we used ICA to minimize the problem. Components related with intravascular artifacts identified were subsequently removed from the ASL data. This makes it feasible to analyze small-world networks with short PLD ASL data by alleviating the effect from intravascular contamination. It should be noted, however, that residual

We are very grateful to Dr. Maria Fernandez-Seara (University of Navarra, Spain) and Dr. David Feinberg (Advanced MRI Technologies, USA) for providing the initial 3D GRASE pCASL sequence. We are also grateful to the National Health and Medical Research Council (NHMRC) of Australia, the Australian Research Council (ARC), and the Victorian Government's Operational Infrastructure Support Program for their support. Appendix A Statistical bootstrap is a non-parametric procedure to estimate the uncertainty of a given statistic (Efron, 1979). A variant of this approach has been previously applied to perfusion MRI to estimate the precision of CBF estimation calculate using dynamic susceptibility contrast MRI (Calamante and Connelly, 2009). In the current study, however, we apply the bootstrap method to assess the statistical significance of the nonlinear fitting of network metrics vs. CBF. In particular, a modelbased resampling approach, named residual bootstrap, fits the data to the model and then residuals are resampled, without the requirement of the repeated acquisitions (Wu, 1986). In brief, the approach and implementation used in our study is as follows: a sigmoid model was employed:

y¼

a þd expð−b  xÞ þ c

ðA1Þ

274

X. Liang et al. / NeuroImage 87 (2014) 265–275

Fig. 6. (A) axial view and (B) sagittal view of connections (Pearson correlation coefficients N 0.75, chosen for improved visualization – see text for further details) and reciprocal of characteristic path length of 84 nodes mapped onto the standard MNI template by using the BrainNet viewer software. Each node is represented by a sphere, with its size inversely proportional to the characteristic path length. The relationship between the 2 network metrics and CBF has been further shown by including region-wise CBF measurements (indicated by the sphere's color, corresponding to the colorbar range (in ml/100 g/min) shown on the right, with hot colors representing high values and cold colors for low values).

where a, b, c and d are the parameters to be fitted, and x and y the observed data. Generally, the relationship between dependent variable y and independent variable x can be described as follows: y ¼ f ðx; pÞ

ðA2Þ

where p is the m vector of adjustable model parameters to be estimated from the observed datasets. ^ is estimated by nonlinear fitting, the raw residual Once the signal y can be calculated as: ^ ^ ε ¼ y−y

ðA3Þ

In general, data fitting can be denoted as: ^ ¼ Hy y

ðA4Þ

where H is the hat matrix. The diagonal elements of the hat matrix are the leverages describing the influence each observed value has on the fitted value. For nonlinear models, H can be calculated as follows (Kuzmic et al., 2004):  −1 T T J ; H¼J J J

ðA5Þ

where J is the Jacobian matrix: J i; j ¼

^Þ ∂f ðxi ; p ∂p j

ðA6Þ

^ is the vector of best-fit model parameters. where p Once H is calculated, the raw residual ^ε can be corrected to yield the m modified residual ^ ε : ^ εi m ^ ffi εi ¼ pffiffiffiffiffiffiffiffiffiffiffiffi 1−hii

ðA7Þ

New bootstrapped residual values εnew were randomly selected m from ^ ε , whereby εnew were added to fitted signal to generate a new ‘bootstrap’ dataset: ^ þ εnew ynew ¼ y

ðA8Þ

By repeating this procedure 1000 times, 1000 independent datasets were obtained, which were subject to nonlinear fitting to assess the statistical significance of R2 values.

Appendix B. Supplementary data Supplementary data to this article can be found online at http://dx. doi.org/10.1016/j.neuroimage.2013.11.013. References Achard, S., Salvador, R., Whitcher, B., Suckling, J., Bullmore, E., 2006. A resilient, lowfrequency, small-world human brain functional network with highly connected association cortical hubs. J. Neurosci. 26, 63–72. Aslan, S., Xu, F., Wang, P.L., Uh, J., Yezhuvath, U.S., Van Osch, M., Lu, H., 2010. Estimation of labeling efficiency in pseudocontinuous arterial spin labeling. Magn. Reson. Med. 63, 765–771. Bassett, D.S., Meyer-Lindenberg, A., Achard, S., Duke, T., Bullmore, E., 2006. Adaptive reconfiguration of fractal small-world human brain functional networks. Proc. Natl. Acad. Sci. U. S. A. 103, 19518–19523. Bonacich, P., 1972. Factoring and weighting approaches to status scores and clique identification. J. Math. Sociol. 2, 113–120. Boxerman, J.L., Bandettini, P.A., Kwong, K.K., Baker, J.R., Davis, T.L., Rosen, B.R., Weisskoff, R.M., 1995. The intravascular contribution to fMRI signal change: Monte Carlo modeling and diffusion-weighted studies in vivo. Magn. Reson. Med. 34, 4–10. Bullmore, E., Sporns, O., 2009. Complex brain networks: graph theoretical analysis of structural and functional systems. Nat. Rev. Neurosci. 10, 186–198. Bullmore, E., Sporns, O., 2012. The economy of brain network organization. Nat. Rev. Neurosci. 13, 336–349. Buxton, R.B., Uludag, K., Dubowitz, D.J., Liu, T.T., 2004. Modeling the hemodynamic response to brain activation. NeuroImage 23 (Suppl. 1), S220–S233. Calamante, F., Connelly, A., 2009. Perfusion precision in bolus-tracking MRI: estimation using the wild-bootstrap method. Magn. Reson. Med. 61, 696–704. Calamante, F., Thomas, D.L., Pell, G.S., Wiersma, J., Turner, R., 1999. Measuring cerebral blood flow using magnetic resonance imaging techniques. J. Cereb. Blood Flow Metab. 19, 701–735. Carhart-Harris, R.L., Erritzoe, D., Williams, T., Stone, J.M., Reed, L.J., Colasanti, A., Tyacke, R.J., Leech, R., Malizia, A.L., Murphy, K., Hobden, P., Evans, J., Feilding, A., Wise, R.G., Nutt, D.J., 2012. Neural correlates of the psychedelic state as determined by fMRI studies with psilocybin. Proc. Natl. Acad. Sci. U. S. A. 109, 2138–2143. Detre, J.A., Leigh, J.S., Williams, D.S., Koretsky, A.P., 1992. Perfusion imaging. Magn. Reson. Med. 23, 37–45. Efron, B., 1979. Bootstrap methods: another look at the jackknife. Ann. Stat. 7, 1–26. Fernandez-Seara, M.A., Wang, J., Wang, Z., Korczykowski, M., Guenther, M., Feinberg, D.A., Detre, J.A., 2007. Imaging mesial temporal lobe activation during scene encoding: comparison of fMRI using BOLD and arterial spin labeling. Hum. Brain Mapp. 28, 1391–1400. Fernandez-Seara, M., Wang, J., Feinberg, D., Detre, J., 2008. Achieving late inflow delay in pseudo-CASL 3D GRASE using a hybridized labeling and background suppression scheme. International Society of Magnetic Resonance in Medicine 16th Scientific Meeting, Toronto, Canada, p. 3343. Fox, M.D., Raichle, M.E., 2007. Spontaneous fluctuations in brain activity observed with functional magnetic resonance imaging. Nat. Rev. Neurosci. 8, 700–711. Garcia, D.M., Duhamel, G., Alsop, D.C., 2005. Efficiency of inversion pulses for background suppressed arterial spin labeling. Magn. Reson. Med. 54, 366–372. Gibbs, M.N., MacKay, D.J.C., 2000. Variational Gaussian process classifiers. IEEE Trans. Neural Netw. 11, 1458–1464. Gong, G., He, Y., Concha, L., Lebel, C., Gross, D.W., Evans, A.C., Beaulieu, C., 2009. Mapping anatomical connectivity patterns of human cerebral cortex using in vivo diffusion tensor imaging tractography. Cereb. Cortex 19, 524–536.

X. Liang et al. / NeuroImage 87 (2014) 265–275 Gonzalez-At, J.B., Alsop, D.C., Detre, J.A., 2000. Cerebral perfusion and arterial transit time changes during task activation determined with continuous arterial spin labeling. Magn. Reson. Med. 43, 739–746. Hagmann, P., Cammoun, L., Gigandet, X., Meuli, R., Honey, C.J., Wedeen, V.J., Sporns, O., 2008. Mapping the structural core of human cerebral cortex. PLoS Biol. 6, e159. Hayasaka, S., Laurienti, P.J., 2010. Comparison of characteristics between region-and voxel-based network analyses in resting-state fMRI data. NeuroImage 50, 499–508. He, Y., Chen, Z., Evans, A., 2008. Structural insights into aberrant topological patterns of large-scale cortical networks in Alzheimer's disease. J. Neurosci. 28, 4756–4766. Heylighen, F., 2001. The science of self-organization and adaptivity. In: Kiel, L.D. (Ed.), The Encyclopedia of Life Support Systems ((EOLSS). Eolss Publishers, Oxford. Kuzmic, P., Hill, C., Janc, J.W., 2004. Practical robust fit of enzyme inhibition data. Methods Enzymol. 383, 366–381. Latora, V., Marchiori, M., 2001. Efficient behavior of small-world networks. Phys. Rev. Lett. 87, 198701. Liang, X., Tournier, J.D., Masterton, R., Connelly, A., Calamante, F., 2012b. A k-space sharing 3D GRASE pseudocontinuous ASL method for whole-brain resting-state functional connectivity. Int. J. Imaging Syst. Technol. 22, 37–43. Liang, X., Wang, J., Yan, C., Shu, N., Xu, K., Gong, G., He, Y., 2012c. Effects of different correlation metrics and preprocessing factors on small-world brain functional networks: a resting-state functional MRI study. PLoS One 7, e32766. Liang, X., Connelly, A., Calamante, F., 2013a. Improved partial volume correction for single inversion time arterial spin labeling data. Magn. Reson. Med. 69, 531–537. Liang, X., Zou, Q., He, Y., Yang, Y., 2013b. Coupling of functional connectivity and regional cerebral blood flow reveals a physiological basis for network hubs of the human brain. Proc. Natl. Acad. Sci. U. S. A. 110, 1929–1934. Liao, W., Zhang, Z., Pan, Z., Mantini, D., Ding, J., Duan, X., Luo, C., Lu, G., Chen, H., 2010. Altered functional connectivity and small-world in mesial temporal lobe epilepsy. PLoS One 5, e8525. Liu, T.T., Brown, G.G., 2007. Measurement of cerebral perfusion with arterial spin labeling: Part 1. Methods. J. Int. Neuropsychol. Soc. 13, 517–525. Luh, W.M., Wong, E.C., Bandettini, P.A., Ward, B.D., Hyde, J.S., 2000. Comparison of simultaneously measured perfusion and BOLD signal increases during brain activation with T(1)-based tissue identification. Magn. Reson. Med. 44, 137–143. Maslov, S., Sneppen, K., 2002. Specificity and stability in topology of protein networks. Science 296, 910–913. Melie-Garcia, L., Sanabria-Diaz, G., Sanchez-Catasus, C., 2012. Studying the topological organization of the cerebral blood flow fluctuations in resting state. NeuroImage 64C, 173–184. Milo, R., Shen-Orr, S., Itzkovitz, S., Kashtan, N., Chklovskii, D., Alon, U., 2002. Network motifs: simple building blocks of complex networks. Science 298, 824–827. Newman, M.E., Watts, D.J., Strogatz, S.H., 2002. Random graph models of social networks. Proc. Natl. Acad. Sci. U. S. A. 99 (Suppl. 1), 2566–2572. Raichle, M.E., Macleod, A.M., Snyder, A.Z., Powers, W.J., Gusnard, D.A., Shulman, G.L., 2001. A default mode of brain function. Proc. Natl. Acad. Sci. U. S. A. 98, 676–682. Rubinov, M., Sporns, O., Thivierge, J.P., Breakspear, M., 2011. Neurobiologically realistic determinants of self-organized criticality in networks of spiking neurons. PLoS Comput. Biol. 7, e1002038. Salvador, R., Suckling, J., Coleman, M.R., Pickard, J.D., Menon, D., Bullmore, E., 2005. Neurophysiological architecture of functional magnetic resonance images of human brain. Cereb. Cortex 15, 1332–1342.

275

Shu, N., Liu, Y., Li, J., Li, Y., Yu, C., Jiang, T., 2009. Altered anatomical network in early blindness revealed by diffusion tensor tractography. PLoS One 4, e7228. Sporns, O., Kotter, R., 2004. Motifs in brain networks. PLoS Biol. 2, e369. Sporns, O., Zwi, J.D., 2004. The small world of the cerebral cortex. Neuroinformatics 2, 145–162. Sporns, O., Chialvo, D.R., Kaiser, M., Hilgetag, C.C., 2004. Organization, development and function of complex brain networks. Trends Cogn. Sci. 8, 418–425. Stam, C.J., Reijneveld, J.C., 2007. Graph theoretical analysis of complex networks in the brain. Nonlinear Biomed. Phys. 1, 3. Stam, C.J., Jones, B.F., Manshanden, I., Van Cappellen Van Walsum, A.M., Montez, T., Verbunt, J.P., De Munck, J.C., Van Dijk, B.W., Berendse, H.W., Scheltens, P., 2006. Magnetoencephalographic evaluation of resting-state functional connectivity in Alzheimer's disease. NeuroImage 32, 1335–1344. Stam, C.J., Jones, B.F., Nolte, G., Breakspear, M., Scheltens, P., 2007. Small-world networks and functional connectivity in Alzheimer's disease. Cereb. Cortex 17, 92–99. Tzourio-Mazoyer, N., Landeau, B., Papathanassiou, D., Crivello, F., Etard, O., Delcroix, N., Mazoyer, B., Joliot, M., 2002. Automated anatomical labeling of activations in SPM using a macroscopic anatomical parcellation of the MNI MRI single-subject brain. NeuroImage 15, 273–289. Vaishnavi, C., 2010. Diagnostic approach to Clostridium difficile infection. Indian J. Gastroenterol. 29, 137–139. Van Den Heuvel, M.P., Sporns, O., 2011. Rich-club organization of the human connectome. J. Neurosci. 31, 15775–15786. Van Den Heuvel, M.P., Stam, C.J., Boersma, M., Hulshoff Pol, H.E., 2008. Small-world and scale-free organization of voxel-based resting-state functional connectivity in the human brain. NeuroImage 43, 528–539. Vlassenko, A.G., Vaishnavi, S.N., Couture, L., Sacco, D., Shannon, B.J., Mach, R.H., Morris, J.C., Raichle, M.E., Mintun, M.A., 2010. Spatial correlation between brain aerobic glycolysis and amyloid-beta (Abeta) deposition. Proc. Natl. Acad. Sci. U. S. A. 107, 17763–17767. Wang, J., Alsop, D.C., Song, H.K., Maldjian, J.A., Tang, K., Salvucci, A.E., Detre, J.A., 2003. Arterial transit time imaging with flow encoding arterial spin tagging (FEAST). Magn. Reson. Med. 50, 599–607. Wang, Z., Aguirre, G.K., Rao, H., Wang, J., Fernandez-Seara, M.A., Childress, A.R., Detre, J.A., 2008. Empirical optimization of ASL data analysis using an ASL data processing toolbox: ASLtbx. Magn. Reson. Imaging 26, 261–269. Wang, J., Wang, L., Zang, Y., Yang, H., Tang, H., Gong, Q., Chen, Z., Zhu, C., He, Y., 2009. Parcellation-dependent small-world brain functional networks: a resting-state fMRI study. Hum. Brain Mapp. 30, 1511–1523. Watts, D.J., Strogatz, S.H., 1998. Collective dynamics of ‘small-world’ networks. Nature 393, 440–442. Wells, J.A., Thomas, D.L., King, M.D., Connelly, A., Lythgoe, M.F., Calamante, F., 2010. Reduction of errors in ASL cerebral perfusion and arterial transit time maps using image de-noising. Magn. Reson. Med. 64, 715–724. Wu, C.F.J., 1986. Jackknife, bootstrap and other resampling methods in regressionanalysis — Discussion. Ann. Stat. 14, 1261–1295. Yang, R.K., Roth, C.G., Ward, R.J., Dejesus, J.O., Mitchell, D.G., 2010. Optimizing abdominal MR imaging: approaches to common problems. Radiographics 30, 185–199.