has been widely used in high-throughput biomedical data analysis with the goal to remove systematic measurement error generated in the complex data acquisition and preprocessing steps and to improve the validity and reproducibility of the following statistical analyses. It has been discussed that a preprocessing step of global signal regression could shift the distributions of the correlations and influence the statistical inferences (Fox et al., 2009; Murphy et al., 2009; Weissenbacher et al., 2009). There may be many other latent factors to affect the quantification of the connectivity metrics as well. Therefore, we feel that normalization toward connectivity metrics should be introduced.

### 5.1. Quantification of Brain Functional Connectivity Metrics

Different from the high-throughput "omics" data, the brain functional connectivity is not measured directly but rather calculated by some statistics/metrics based on a pair of time courses from fMRI data. It is unclear how the calculated statistics/metrics can appropriately reflect the true connectivity strength and are comparable across subjects, regardless what statistic is chosen (e.g., correlation coefficient or mutual information coefficient). It is possible to obtain extremely large absolute value correlations between two white noise vectors, which gives rise to the false positive discovery. From the statistical perspective, most connectivity statistics can be proved to follow a known distribution asymptoticly and accordingly the p-values are calculated with both type I and II errors. Comparing with the conventional normalization method such as quantile normalization, the empirical Bayes mixture model lends itself to incorporating the false positive concept into quantification of the functional connectivity expression and provides a (posterior) probability based scale. The data driven (rather than a deterministic linear/nonlinear transformation) quantification method could provide a more comparable scale for group level connectivity inferences. For example, a 0.1 difference in original correlations could be mapped to around 0.5 difference in the normalized correlations at the interaction between two components due to the increase of posterior probability of true positive. The amplified difference tend to improve the subtle difference detection, because it can better represent connectivity strength. The computational techniques for the mixture model estimation have been developed for local fdr estimation by Efron (2004) and Wu et al. (2006), which provides us a convenient tool to calculate the subject-specific normalization function. The only assumption of our method is that the majority (p0 > 0.9) of connectivity expressions are from the null distribution, which needs to further verified with more rs-fMRI studies. The assumption is generally valid, and all connectivity metric distributions of the data sets we tested follow such pattern. If the assumption is violated, Wu et al. (2006) provides promising numerical solution using nonparametric curve fitting methods. Moreover, another obvious advantage of the normalization method is that it maps the correlations to the range of [0, 1] by the empirical Bayes posterior probability normalization function, which avoids the information loss due to hard thresholding of correlations in complex network analysis using graph theoretical models (Rubinov and Sporns, 2011).

The appropriate brain connectivity metric normalization method improves the power to detect the truly differentially expressed features and yield less false positive findings. In the simulation study, we compare the test results based on different connectivity metrics with reference to ground truth, and it shows the empirical Bayes normalized correlation has the lowest type I and II errors and is more robust to systematic shifts. When applying our method to the data example, the analysis results based on normalized connectivity metrics detect hyper-connectivity between pairs of regions from distant functional subsystems for the ASD group with comparing to TC group. Such features are not detected by using the nonnormalized correlations. The findings align with the results by Supekar et al. (2013) which performs between region connectivity analysis for several autism studies from different sites. Supekar et al. (2013) also provides explanation of these findings from the perspectives of neuroscience and the link to clinical symptoms of ASD. The practical brain connectivity study using neuroimaging technology often involves multiple steps of numerical analysis which are subject to many unavoidable errors and noises, and we feel that the empirical Bayes normalization improves both power and reliability of statistical analysis.

### 5.2. Anticorrelations

The anticorrelations in rs-fMRI data have drawn attention of many neuroimaging researchers (Fox et al., 2009; Murphy et al., 2009; Weissenbacher et al., 2009; Chai et al., 2012). The discussion has not reached to the agreement whether the anticorrelations are "true positive" or "false positive." The proposed normalization method provides a pathway to automatically detect the "true positive" anticorrelation

component by classifying the "true positives" and "false positives" based on the empirical distribution of connectivity metric. **Figure 8** shows that correlated and anticorrelated components can be identified, if existing, we could assign either "+" or "−" sign to anticorrelated connectivity metric depending on different following analyses. Generally, "−" sign suits the regression analysis or statistical tests better, because anticorrelation could be considered as the opposite of correlation. When applying the graph theoretical model based network analysis using normalized connectivity, two separate analyses should be conducted for correlations and anticorrelations (with "+" sign) if both components are detected, with the normalized connectivity metric range of [0, 1] (in **Figure 8B**). Thus, the results include two parts of inferences: properties of correlated networks and anticorrelated networks. Although in our data example there is no anticorrelation component detected, that normalization method can be also applied to deal with anticorrelations in practical data analysis. Yet, out normalization method could be combined with pre-processing steps (e.g., global signal regression), as the normalized connectivity is probability and shift-invariant.

## 6. Conclusion

In summary, a new rs-fMRI connectivity metric normalization method has been developed and applied to functional

### References

brain connectivity analysis. The better connectivity normalization/quantification methods yield generally higher reproducibility. Although we utilize the Pearson correlation coefficient as connectivity metric and rs-fMRI for demonstration, we are optimistic that the developed method are ready to be applied to the task-induced fMRI connectivity study and other connectivity metrics because the empirical Bayes framework is flexible to fit various distributions of connectivity metrics.

### Acknowledgments

The authors would thank the researchers of the ABIDE project for sharing their clinical and rs-fMRI data at http://fcon\_1000. projects.nitrc.org/indi/abide/. Chen's research is supported in part by UMD Tier1A seed grant. Kang's research was partially supported by NIH grant 1R01MH105561. The authors also want thank Dr. Luiz Pessoa from University of Maryland, College Park and Dr. F. DuBois Bowman from Columbia University for constructive discussions.

# Supplementary Material

The Supplementary Material for this article can be found online at: http://journal.frontiersin.org/article/10.3389/fnins. 2015.00316

Mesulam, M. M. (1998). From sensation to cognition. Brain 121, 1013–1052.

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

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

# A Parcellation Based Nonparametric Algorithm for Independent Component Analysis with Application to fMRI Data

### Shanshan Li 1 \*, Shaojie Chen2 , Chen Yue2 and Brian Caffo2

*1 Department of Biostatistics, Indiana University Fairbanks School of Public Health, Indiana University, Indianapolis, IN, USA, 2 Department of Biostatistics, Johns Hopkins Bloomberg School of Public Health, Johns Hopkins University, Baltimore, MD, USA*

### Edited by:

*Bertrand Thirion, Institut National de Recherche en Informatique et Automatique, France*

### Reviewed by:

*Xi-Nian Zuo, Chinese Academy of Sciences, China Alexis Roche, Siemens Healthcare/Centre Hospitalier Universitaire Vaudois, Switzerland Ronald Phlypo, Grenoble Images Parole Signal Automatique-Lab, France*

\*Correspondence:

*Shanshan Li<EMAIL_ADDRESS>
### Specialty section:

*This article was submitted to Brain Imaging Methods, a section of the journal Frontiers in Neuroscience*

Received: *14 May 2015* Accepted: *11 January 2016* Published: *29 January 2016*

### Citation:

*Li S, Chen S, Yue C and Caffo B (2016) A Parcellation Based Nonparametric Algorithm for Independent Component Analysis with Application to fMRI Data. Front. Neurosci. 10:15. doi: 10.3389/fnins.2016.00015* Independent Component analysis (ICA) is a widely used technique for separating signals that have been mixed together. In this manuscript, we propose a novel ICA algorithm using density estimation and maximum likelihood, where the densities of the signals are estimated via p-spline based histogram smoothing and the mixing matrix is simultaneously estimated using an optimization algorithm. The algorithm is exceedingly simple, easy to implement and blind to the underlying distributions of the source signals. To relax the identically distributed assumption in the density function, a modified algorithm is proposed to allow for different density functions on different regions. The performance of the proposed algorithm is evaluated in different simulation settings. For illustration, the algorithm is applied to a research investigation with a large collection of resting state fMRI datasets. The results show that the algorithm successfully recovers the established brain networks.

Keywords: blind source separation, density estimation, functional MRI, p-spline bases, signal processing

# 1. INTRODUCTION

This manuscript puts forward two innovations. Firstly, we demonstrate a fast, likelihood motivated and straightforward method for applying independent components analysis (ICA). Secondly, we propose a parcellation based adjustment when the source signals distribute differently across regions. Our work is routed in the context of understanding human brain networks, and we use functional magnetic resonance imaging (fMRI) data for illustration in this manuscript.

We approach our study of fMRI by simultaneously analyzing all voxels. This is in contrast to regional or seed-based approaches (Buckner et al., 2005; Wang et al., 2006; Allen et al., 2007) that restrict attention to carefully chosen locations. Such approaches require strong assumptions on the choice of seeds or parcellation used to define region. Hence voxel-wise approaches are important complementary procedures. Given the volume of voxels under study (usually on the order of fifty thousand non-background ones), flexible yet parsimonious models approaches are required. However, even with parsimonious models, whole brain voxel-level techniques are more empirical and exploratory than their more hypothesis driven regional and seed-based counterparts. Thus, exploratory factor-analytic models are common approaches in voxel-level investigations.

Independent components analysis (ICA) is a factor-analytic approach that has been frequently utilized for the analysis of functional neuroimaging data, because of its success in discovering key benefits of ICA are its exploratory nature and its often considered reasonable underlying generative model. Specifically, it models collected signals, X, as linear weighted combinations of independent sources, S1, S2, ...Sp. Thus, we can write the noisefree ICA model as X = AS, where S = [S1, S2, ..., SQ] and A is a Q × Q full rank matrix, the so-called mixing matrix. The goal of ICA is to recover the underlying signals S1, S2, ..., SQ from their observed mixtures X1, X2, ..., XQ. Note that, in the context of fMRI, the independent components S1, S2, ..., SQ are often interpreted as brain networks and A is the mixing matrix characterizing the temporal pattern of the corresponding brain networks.

Various algorithms for ICA have been proposed in the literature. See Hyvärinen et al. (2001); Comon and Jutten (2010); Risk et al. (2014) for comprehensive reviews. One common procedure is to postulate a parametric family for the source distributions and then obtain the independent components (ICs) by optimizing a contrast function that measures the distribution property of the output (Samworth and Yuan, 2012). The contrast functions can be selected based on different measures, such as entropy, mutual independence, high-order decorrelations, divergence between the joint distribution of the output and some model, etc. (Cardoso, 1998). These include the popular FastICA algorithm (Hyvärinen and Oja, 2000), the JADE algorithm (Cardoso, 1999), the Pearson ICA algorithm (Karvanen et al., 2000), and a few other algorithms proposed by Comon (1994); Amari and Cardoso (1997); Li and Adali (2010). An alternative procedure is to assume smooth densities for the source distributions and use nonparametric or semiparametric approaches to estimate those density functions. The mixing matrix can then be derived using maximum likelihood method. For example, Bach and Jordan (2003) developed a nonparametric estimation approach based on canonical correlations in a reproducing kernel Hilbert space. Hastie and Tibshirani (2002) expressed the source distribution by an exponentially tilted Gaussian function and used the fixed-point algorithm for estimation of the mixing matrix. Boscolo et al. (2004) used kernel density estimation techniques to model the underlying densities and quasi-Newton method for optimization. Guo and Pagnoni (2008) used Gaussian mixture models for the source distribution and provided an expectation-maximization (EM) framework for estimation, assuming Gaussian noise in the model. Eloyan et al. (2013) estimated the source distribution by using mixture density estimates, and proposed a constrained EM algorithm for estimation.

The benefit of the likelihood-based ICA algorithm is that, as a byproduct of the ICA algorithm, one obtains the fully specified likelihood of the ICA model which can be used for further statistical inference. For example, based on the fully specified likelihood, one can conduct Bayesian analysis or perform likelihood based model selection. However, the existing likelihood-based ICA algorithms are mostly semi-parametric and are usually computationally intensive. In this manuscript, we aim to develop a likelihood-based algorithm that is exceedingly simple and truly blinded to the source distributions.

We propose to estimate the density function of the ICs via histogram smoothing, following a well-known approach in the penalized spline literature. At its core, likelihood-based ICA requires estimation of the mixing matrix and flexible density estimation for the ICs. Our approach, like many other likelihood-based approaches, iteratively estimates these components separately using block maximization. In contrast to other approaches, we use an exceedingly simple density estimation technique via histogram smoothing. Specifically, we assume the bin counts of the frequency histogram follow a Poisson distribution and express the mean counts as sum of B-spline bases via generalized linear model. To smooth the histogram, we follow Eilers and Marx (1996) to construct a penalized likelihood with a difference penalty on coefficients of adjacent B-splines. Apart from its simplicity, a benefit of this approach is speed. Density estimation and evaluation for tens of thousands of voxels is time consuming, and worse, is performed within an iterative algorithm. Using histogram smoothing, the voxel-level calculation reduces to estimating a histogram, a very fast process.

We briefly mention that, in our primary area of application, fMRI, we focus entirely on noise-free group spatial independent component analysis. By assuming noise-free model, noise in the data is absorbed into the estimated ICs and the mixing matrix. By using spatial ICA model, the fMRI data is decomposed into spatial maps multiplied by their respective time courses, where the maps are drawn from spatial distributions that are statistically independent (Calhoun et al., 2001a). The spatial independence assumption is well suited to the sparse nature of the spatial pattern for typical brain activation (McKeown and Sejnowski, 1998; Guo and Pagnoni, 2008). The time courses estimated from spatial ICA describe the temporal characteristics of functional networks, i.e., areas of temporal correlation in the BOLD signal. For multi-subject fMRI data, we assume common spatial maps for all subjects and subject-specific mixing matrices, therefore, we can concatenate all subjects' data in the temporal domain, and apply ICA to the aggregated data matrix. The group mixing matrix is the concatenated time course for all subjects. Individual mixing matrices can be backreconstructed by partitioning the group mixing matrix into submatrices corresponding to each subject.

The remainder of the paper is organized as follows. Section 2 describes the p-spline based ICA algorithm and considers relaxation of the i.i.d signal assumption. Section 3 shows the performance of the proposed algorithm in simulation study. Section 4 provides the application of the proposed algorithm to the 1000 Functional Connectome Project (https://www.nitrc.org/ projects/fcon\_1000/), while Section 5 gives a discussion.

# 2. METHODS

# 2.1. Description of ICA Methodology

Independent component analysis models collected signals as linear weighted combinations of independent sources. Notationally, let Xi be a T × V matrix for subject i = 1, ...,I. In the context of fMRI, T indicates scans while V indicates voxels. Assume the number of ICs is Q. The ICA model specifies Xi = AiS, where Ai is a T × Q mixing matrix and S is a Q × V matrix of ICs. By assuming common spatial maps across subjects, we can stack the individual matrices in the temporal domain. Let X = [X T 1 , X T 2 , ..., X T I ] T be the TI × V group data matrix, and A = [A T 1 , A T 2 , ..., A T I ] T be the TI × Q group mixing matrix. Spatial group ICA simply specifies the standard model

$$X = A\text{S.}\tag{1}$$

We use parentheses to index matrices so that X(t, v) is element (t, v) of X and define X(t, ·) as row t of X and X(·, v) as column P v. Then, model (1) could be rewritten as X(t, v) = Q q=1 A(t, q)S(q, v) and X = PQ q=1 A(·, q)S(q, ·).

We assume that E[X] = µx = 0 and hence E[S] = µs = 0. If this assumption were not made, the ICA model would imply X − µx = A(S − µS), which is exactly an ICA model with a centered data matrix and the ICs having mean 0. Hence, X is demeaned prior to analyses and µS is assumed to be zero. Similarly, since A(·, q)S(q, ·) = {A(·, q)/c} ∗ {cS(q, ·)}, ICs are only identified up to scalar multiplication. Thus, we assume that Var{S(q, v)} = 1 for q = 1, . . . , Q and v = 1, . . . ,V.

ICA gets its name by assuming that S(q, ·) á S(q ′ , ·) when q 6= q ′ , where á implies statistical independence. However, standard variations of ICA also assumes that {S(q, v)} V v=1 is an i.i.d collection, which we also adopt for now. The i.i.d assumption will be relaxed later in the next subsection. As a consequence of these assumptions, X(·, v)áX(·, v ′ ) when v 6= v ′ ; yet note that X(t, ·) is not (necessarily) independent of X(t ′ , ·).

Typically, Q < TI and Equation (1) is overdetermined. A two-stage dimension reduction is often performed to reduce the computational load and avoid overfitting (Calhoun et al., 2001a; Beckmann and Smith, 2005; Guo and Pagnoni, 2008; Eloyan et al., 2013; Risk et al., 2014). Specifically, in the first stage, an SVD is performed in the temporal domain within subject, where the first R eigenvectors are retained. The dimension for the group data matrix then becomes RI × V. In the second stage, an SVD is performed on the group data matrix obtained from the first stage and the first Q eigenvectors are retained to force a determined linear system for the group ICA model. This discards information in the data. However, one hopes that by selecting the first Q singular vectors, the most relevant features of the data will be retained. The choice of R and Q could be based on various criteria, including variance explained, informationtheoretic criteria, and practical considerations. This is not a major concern in this article.

## 2.2. ICA Through Fast Nonparametric Density Estimation

ICA estimates S by seeking an unmixing matrix, say Bˆ, such that BXˆ is a good approximation to the original sources S. Let B = A −1 be the estimand of interest. Notationally following Hyvärinen et al. (2001), if fq is the density for S(q, v) for v = 1, . . . ,V, and f = (f1, ..., fQ), then standard multivariate random variable transformation results imply that the joint density of X(·, v) is

$$\begin{aligned} \{\mathcal{S}(\mathcal{X}(\cdot,\nu))\} &= |\det(\mathcal{B})| \prod\_{q=1}^{Q} f\_q \{\mathcal{S}(q,\nu)\} \\ &= |\det(\mathcal{B})| \prod\_{q=1}^{Q} f\_q \{\mathcal{B}(q,\cdot)\mathcal{X}(\cdot,\nu)\}, \end{aligned}$$

therefore the joint log-likelihood including all contributions for v = 1, . . . ,V is

$$\mathcal{L}(B, f) = \sum\_{\nu=1}^{V} \sum\_{q=1}^{Q} \log[f\_q(B(q, \cdot)X(\cdot, \nu))] + V \log|\det(B)|.$$

It is generally not possible to solve the joint likelihood for the parameters in fq and B simultaneously. Instead, an iterative optimization is often performed. Specifically, given the current estimate of B at iteration k, say Bˆ (k) , one can get an estimate for S via Sˆ (k) = Bˆ (k)X. Given Sˆ (k) (q, ·), density estimation techniques can be used to obtain ˆ f (k) q , the estimate of ˆ fq at iteration k.

We suggest the use of histogram smoothing as the density estimation technique, where the bin counts of the frequency histogram are assumed to follow a Poisson distribution and a penalized likelihood is constructed to produce a smooth density estimate. The details of histogram smoothing can be found in Eilers and Marx (1996), and we provide a sketch below. (Readers not familiar with statistical smoothing may skip the rest of this paragraph). Notationally, let c (k) (q, 0) < c (k) (q, 1) < . . . < c (k) (q, J) be equidistant histogram cutpoints, where c (k) (q, 0) = −ǫ+min Sˆ (k) (q, ·) and c (k) (q, J) = ǫ+max Sˆ (k) (q, ·). The number ǫ is added to avoid numerical boundary effects. Let n (k) P (q, j) = V v=1 I{c (k) (j − 1) < Sˆ (k) (q, ·) ≤ c (k) (j)}, for j = 1, . . . , J, be the count of values between cutpoints j − 1 and j for row q of Sˆ (k) . Define the midpoints of intervals [c (k) (q, j − 1),c (k) (q, j)] by m(k) (q, j) for j = 1, . . . , J. We obtain a density estimate via the log-linear model n (k) (q, j) ∼ Poisson{λ (k) (q, j)}, where log{λ (k) (q, ·)} = PL l=1 D (k) {m(k) (q, ·), l}β (k) (q, l). Here the log function is presumed to act component-wise on vectors, D (k) is a B-spline basis design matrix, L is the number of knots for B-splines, and β (k) (q, ·) is a vector of coefficients. To avoid overfitting the B-spline model, and to avoid sensitivity to the degrees of freedom, we choose a large value for the degrees of freedom and put a squared penalty on the coefficients. Let µ (k) (q, j) denote the expectation of n (k) (q, j), then the penalized log likelihood takes the form (Eilers and Marx, 1996)

$$\begin{aligned} \mathcal{L} &= \sum\_{j=1}^{f} n^{(k)}(q,j) \ln \mu^{(k)}(q,j) - \sum\_{j=1}^{f} \mu^{(k)}(q,j), \\ &- \delta \sum\_{l=3}^{L} \frac{\{\Delta^2 \beta^{(k)}(q,l)\}^2}{2}, \end{aligned}$$

where δ is a parameter controlling the smoothness of the fit, 1 denotes the difference operator, 12β(·, l) = β(·, l) − 2β(·, l − 1) + β(·, l − 2). The resulting model is then a generalized linear mixed model on the counts. The B-spline basis is evaluated at the midpoint of the cutpoint interval. However, via interpolation, the smoother gives an estimate for all values, thus yielding a continuous function, say ˆ f (k) q (s), which is the density estimate.

Using generalized linear mixed models to penalize smoothing has become standard practice and is well described in Ruppert et al. (2003). Histogram smoothing as a density estimate appears to be less commonly used. However, we note that this pspline based density smoother has very attractive properties (Eilers and Marx, 1996). First, it results in a proper density. Secondly, it elegantly handles boundary issues, unlike other density estimators (such as kernel density estimator). Thirdly, the estimated density conserves the first few empirical moments (means and variances) of the histogram, depending on the order of the B-splines. More details regarding these properties can be found in Eilers and Marx (1996). Note that, conservation of moments is an important property that guarantees the identifiability of the ICA model. We choose a cubic B-spline which then conserves the first two moments of the histogram.

Furthermore, due to the convenient differentiation properties of B-spline bases and the simple exponential (Poisson) model, the first and second derivatives of ˆ f (k) q are immediately available, where d ˆ f (k) q = exp{ ˆ f (k) q }β (k) (q, ·)dD(k) . Thus, derivatives of L(B) are available in closed form, making gradient- and Hessian-based optimization algorithms easy to implement. This is useful for the stage of the algorithm for obtaining the next iterate of B. Accordingly, we use a Newton-Raphson method to update the mixing matrix. Specifically, let L ′ and L ′′ denote the first and second derivatives of the log likelihood. At the kth iteration, we update B by

$$B^{(k+1)} = B^{(k)} - \mathcal{L}^{\prime\prime} (B^{(k)})^{-1} \mathcal{L}^{\prime} (B^{(k)}).\tag{2}$$

The starting values of B should satisfy the condition that the underlying ICs are the same for all subjects. Following Eloyan et al. (2013), we decompose the full matrix X using the population value decomposition X = U6V T (Crainiceanu et al., 2011), and the starting values of the Bi are chosen as the ith block of the rows of U6. Thus, given a starting value for B, histogram smoothing is used to obtain fq, then given the update for fq, the natural gradient algorithm is used to obtain B and these steps are iterated until convergence. Let P denote B (k) (B (k+1)) −1 . We use the Amari metric between B (k+1) and B (k) as our convergence criterion (Amari, 1998), where the metric is defined as

$$\begin{aligned} d\langle B^{(k)}, B^{(k+1)} \rangle &= \frac{1}{2Q} \sum\_{i=1}^{Q} \left( \sum\_{j=1}^{Q} \frac{|P\_{ij}|}{\max\_{j} |P\_{ij}|} - 1 \right) \\ &+ \frac{1}{2Q} \sum\_{j=1}^{Q} \left( \sum\_{i=1}^{Q} \frac{|P\_{ij}|}{\max\_{i} |P\_{ij}|} - 1 \right). \end{aligned}$$

The Amari metric is useful, as it is invariant to permutation of the ordering of the ICs, a necessary condition for a convergence metric to be useful.

### 2.3. ICA Based on Parcellation

Most ICA algorithms (such as the commonly used fastICA, JADE, etc.) assume that {S(q, v)} V v=1 is an i.i.d collection for all q = 1, ..., Q. Intrinsically, this is to assume that the values of the ICs are independent draws from a density. The i.i.d assumption is made for simplicity, but it may not hold for fMRI data. Calhoun et al. (2001b) considered possible violations of the independence assumption for task-based fMRI data. They found that the ICA algorithm was successful when the correlation in the signal was small, but it may fail when the signals are highly dependent. However, for most task-based fMRI and resting-state fMRI data, the correlation between voxels is negligible. Therefore, we do not pursue the approach to deal with violation of the independence assumption here. Instead, we consider relaxation of the identically distributed assumption.

Specifically, we propose to account for the difference in the activity across the brain by allowing different density distribution in different regions. To this end, we adopt the functional parcellation of the brain activity map proposed by Yeo et al. (2011). The parcellation includes 17 functional networks in the cerebral cortex, that is, I = 18 ROIs for the whole brain. We assume the signals are i.i.d within region but could be differently distributed across region. Under this assumption, the density function fq can be written as the sum of the region-specific density function, that is,

$$f\_q(s) = \sum\_{i=1}^{I} I(s \in R\_i) f\_{iq}(s),$$

where Ri denotes the ith ROI, fiq is the density function on Ri . Thus, fiq takes positive values on the ith region and zero elsewhere. The density estimate of fiq can be obtained using the same procedure as proposed in Section 2.2, confined to the ith region. The estimate for fq can be constructed by taking the sum of ˆ fiq. The rest of the ICA algorithm follows the proposed procedure in Section 2.2.

The proposed ICA algorithm can be summarized as follows:

	- a. Let S = BX.
	- b. For each IC q, calculate the density function fiq(s) on the ith ROI, i = 1, 2, ...,I, using the p-spline based density estimation algorithm.
	- c. Get fq(s) = PI i=1 I(s ∈ Ri)fiq(s).
	- d. Update the mixing matrix B using the Newton-Raphson method, see Equation (2).

Note that, in the special case that f1q = f2q = ... = fIq, the above algorithm reduces to the algorithm proposed in Section 2.2 assuming i.i.d signals across the entire brain.

### 3. SIMULATION

We conduct simulation studies to evaluate the performance of the proposed ICA algorithm. We consider four settings where

### TABLE 1 | The average computation time (in seconds) per simulation replication using different algorithms in the simulation study.

data are generated using different distributions. We compare the results of the proposed algorithm with fastICA (Hyvärinen et al., 2001), JADE (Cardoso, 1999), Pearson ICA (Karvanen et al., 2000), ProDenICA (Hastie and Tibshirani, 2002), and HDICA (Eloyan et al., 2013). We implement the algorithms fastICA, ProDenICA, JADE, PearsonICA using the R packages "fastICA" (Marchini et al., 2013), "ProDenICA" (Hastie and Tibshirani, 2010), "JADE" (Nordhausen et al., 2014), and "PearsonICA"

(Karvanen, 2006). The proposed p-spline based ICA algorithm and the HDICA (Eloyan et al., 2013) are also implemented in R.

The computation environment is a multi-core Linux cluster with more than 680 cores running in the average of 2.5 GHz speed and 4.4 TB of memory. On average, the contrast-function based algorithms (fastICA, PearsonICA, JADE) perform much faster than the likelihood-based algorithms (p-spline ICA, ProDenICA, HDICA). (See **Table 1** for a summary of the computation time using different algorithms.) However, since those are essentially two different sets of algorithms, we restrict the comparison of the computational intensity within the category of likelihood-based algorithm.

In the first set of simulation studies, we assume there are Q = 3 independent components, and they are generated by S(1, ·) ∼ Weibull(1, 1), S(2, ·) ∼ Gamma(1, 1), and S(3, ·) ∼ Gamma(2, 2), respectively. Standard Gaussian noises are added to the generated ICs. The mixing matrix is given by

$$A = \begin{pmatrix} 2 & 1 & 2 \\ 3 & 3 & 1 \\ 1 & 2 & 2 \end{pmatrix}.$$

**Figure 1** summarizes the simulation results based on 200 replications. The spatial correlation is the absolute correlation between the estimated spatial map and the true spatial map without noise. The Amari error is computed to evaluate the accuracy of the estimated mixing matrix (Amari, 1998). It is seen from the boxplots of the spatial correlation and the Amari errors that the proposed ICA algorithm performs equally well as fastICA, JADE, and PearsonICA, and all

these algorithms perform substantially better than ProDenICA algorithm. ProDenICA fails probably due to the extreme values introduced by the noise (See more discussion in Risk et al., 2014). This shows that the ProDenICA is sensitive to extreme values, while our algorithm is robust to extreme values. The average computation time per replication is 6.21 s using the

FIGURE 4 | The underlying signals for the fourth simulation setting: ICs 1, 2, and 3 (left to right).

p-spline ICA, 3.12 s using ProDenICA, and 308.34 s using HDICA.

In the second setting, we assume the number of source signals Q = 2, and we generate the signals based on parcellation. Specifically, we partition the real line into 10 intervals, with cutoffs at the 10th, 20th, 30th, 40th, 50th, 60th, 70th, 80th, and 90th percentiles of the normal distribution. For the first IC, the density function is uniformly distributed within each interval, but the overall shape is approximately normal. For the second IC, the density function follows Laplace distribution within each interval, and the overall shape is approximately normal. The mixing matrix is given by

$$A = \begin{pmatrix} 2 & 1\\ 3 & 2 \end{pmatrix}.$$

The boxplots of the spatial correlation and the Amari errors based on 200 replications are summarized in **Figure 2**. Under the second scenario, the underlying signals have region specific densities, and the overall density functions for both components are approximately normal. All the competing algorithms considered in the simulation studies show substantial bias. These algorithms fail to recover the true signals because they heavily depend on the non-gaussianity assumption. On the contrary, the proposed algorithm accounts for the effect of parcellation and recovers the true signals with relatively high accuracy. The proposed algorithm substantially outperforms all the competing algorithms under the second setting. The average computation time per replication is 5.89 s using p-spline ICA, 1.67 s using ProDenICA, and 76.59 s using HDICA.

In the third setting, we generate multi-subject data with number of subject I = 3. The source signals are the same as those in the second setting, and the mixing matrices for the three subjects are given by

$$A\_1 = \begin{pmatrix} 1 & 0.25 \\ 0.5 & -0.5 \end{pmatrix}, \ A\_2 = \begin{pmatrix} 1 & 2 \\ 0.5 & -0.5 \end{pmatrix}, \ A\_3 = \begin{pmatrix} 1 & 0.5 \\ 0.5 & 2 \end{pmatrix}.$$

The simulation results are summarized in **Figure 3**, where in each simulation replication, the Amari error is calculated as the average of the Amari errors for all three subjects. The results show that, for multi-subject data, the proposed algorithm successfully recovers both the common spatial signals and the individual mixing matrices. In addition, for similar reasons as in the second setting, the proposed algorithm substantially outperforms all the competing algorithms.

In the fourth setting, we generate the ICs and mixing matrices by mimicking signals from real fMRI data. Specifically, we run fastICA on 10 subjects from the NITRC 1000 Connectome dataset to get twenty ICs (networks). Three of the twenty networks are chosen as the true signals, and they are shown in **Figure 4**. The time courses are also signals from real data, obtained in a similar way as in Calhoun et al. (2009). They are shown in **Figure 5**. We first apply a two-stage dimension reduction using the method as described in Section 2.1. Then we apply the proposed atlas-based ICA algorithm using the brain parcellation proposed by Yeo et al. (2011). The correlation matrix

between the true signals and the estimated signals using the proposed algorithm is

$$
\begin{pmatrix}
0.028 & 0.999 & -0.006 \\
\end{pmatrix}.
$$

The results indicate that our proposed p-spline based ICA algorithm is successful in recovering signals from real fMRI data.

### 4. APPLICATION

We apply our proposed algorithm to the 1000 Functional Connectomes Project dataset, which consists of thousands of resting state scans combined across multiple sites with the goal of facilitating discovery and analysis of brain networks (Biswal et al., 2010). It is one of the largest freely available fMRI datasets. The fMRI scans were collected when the subjects stayed in the scanner for 2.2–20 min in resting state. Scanning parameters used to acquire the data from each site are detailed elsewhere (for complete information see https://www.nitrc.org/projects/fcon\_ 1000/).

As the quality and scanning parameters vary across sites, we focus on data from the largest site, Cambridge, which contains I = 50 subjects. For the subjects used in this analysis, the number of time points is T = 119. We use the MNI template to remove the background noise and to retain voxels that are in the actual brain. For each subject, we have a T × V dimensional matrix Xi . The group data matrix X is obtained by concatenating I subjects' data in the temporal domain.

Following Biswal et al. (2010), we assume there are Q = 20 independent components in this application. An SVD is

performed to reduce the dimension of the aggregated data matrix to Q×V. The ICA algorithms are then applied to the reduced data matrix and the Python toolbox Nilearn (Abraham et al., 2014) is used for visualization of the estimation results. Specifically, the estimated ICs using the proposed p-spline based ICA algorithm are shown in **Figure 6**. Several main brain networks including the default mode network (DMN) and the control network are successfully identified by the proposed algorithm. As a comparison, the results from fastICA are shown in **Figure 7**. The ICs estimated by fastICA and the p-spline ICA are matched by correlation. Of the 20 pairs, the highest correlation is 0.99, the lowest correlation is 0.52, and the median correlation is 0.93. Specifically, the following is a list of these correlations for the major brain networks: visual network (0.99), auditory network (0.98), DMN (0.96), and control network (0.92).

As suggested by an anonymous reviewer, we investigate the impact of the dimension of the reduced space on the final results. Specifically, we select different values of R and Q (the number of eigenvectors in the subject-level and group-level dimension reduction), and rerun the ICA algorithm on the dimension

reduced dataset. We set R = 15, 20, 30 and Q = 15, 20, 30, respectively. Similarly as in Li et al. (2007), we find that the IC estimates are well separated when Q = 15, 20. When Q = 30, the estimation of the major networks shows degradation and a few of the other estimated components seem to be noise. Specifically, the correlations for the major brain networks using R = 20, Q = 15, and R = 20, Q = 20 are as follows: visual network (0.96), auditory network (0.73), DMN (0.86) and control network (0.84). In addition, the correlations for the major brain networks using R = 20, Q = 30, and R = 20, Q = 20 are as follows: visual network (0.78), auditory network (0.61), DMN (0.88) and control network (0.69). In summary, we find that, although the estimation results depend on the number of components, the major networks appear to be robust against the choices of number of components.

### 5. DISCUSSION

Independent component analysis is a factor-analytic approach that is commonly used in analyzing fMRI data. In this manuscript, we present a novel and simple ICA algorithm that is fast, likelihood based and straightforward to program. The algorithm is nonparametric, data-driven, and is blind to the particular distribution of the underlying signals. As a byproduct of the algorithm, we obtain the likelihood function of the ICA model which can be used for further statistical inference. It should be noted that, the likelihood function in our algorithm is a profile likelihood, since we are mainly interested in the mixing matrix estimates and the parameters over the spline basis are nuisance parameters. Indeed, one could also study the coefficients on the spline basis in a full likelihood, but this is not the goal of this manuscript, hence the variance of the estimator of the mixing matrix depends on the variance of the nuisance parameters.

The proposed algorithm is extended to allow for region specific IC density functions, on the rationale that most signals of interest are reasonably confined to a subset of the entire anatomical brain space (Guo and Pagnoni, 2008). When the source signals distribute identically across brain, the estimation accuracy of the parcellation-based estimator is similar to that of the full-brain estimator, because it becomes equivalent to the full-brain estimator. However, when the source signals distribute differently across brain, the full-brain estimator may result in substantial bias while the parcellation-based estimator can successfully recover the source signals. It should be noted that the parcellation based adjustment can be applied to other ICA algorithms as well. Indeed, for any gradient-based ICA, one can do the adjustment by taking a weighted sum over the updates of each of the parcellations, where the weights account for the number of samples in the parcellations. This flexibility ensures the generalizability of the proposed parcellation based adjustment.

Simulation studies show that our proposed algorithm works well in both the simple and complex situations, and it

### REFERENCES

substantially outperforms the existing ICA algorithms when the identically distributed assumption of the source signals is violated. In applying the proposed algorithm to the fMRI data, we choose to account for the difference in brain activities across regions by using the brain parcellation proposed by Yeo et al. (2011). Our data application results show that the proposed algorithm successfully identifies the main brain networks in the 1000 Functional Connectomes Project dataset.

There are a few directions for future research. Firstly, the testretest reliability of the intrinsic brain networks is an important issue and has been studied extensively in recent years. For example, Zuo et al. (2010) found that a few functionally relevant components (such as the default mode, auditory-motor and executive control) show the highest reliability across all components. It would be interesting to compare different ICA algorithms in identifying and characterizing those functionally relevant components. Secondly, there are a variety of existing brain parcellation schemes, including those proposed by Tzourio-Mazoyer et al. (2002); Fischl et al. (2004); Beckmann et al. (2009); Yeo et al. (2011). It would be interesting to study the optimal choice of parcellation under different scientific scenarios. Thirdly, as pointed out by an anonymous reviewer, pre-whitening, although a standard pre-processing procedure, may result in loss of information and bias in estimation (Cardoso, 1994). It would be interesting to investigate alternative preprocessing procedures to avoid the bias introduced by prewhitening.

### ACKNOWLEDGMENT

This work was partly supported by a grant from the Simons Foundation (#354180, Shanshan Li).

during the stroop color-naming task. Proc. Natl. Acad. Sci. U.S.A. 95, 803–810. doi: 10.1073/pnas.95.3.803

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

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

# Identification of Voxels Confounded by Venous Signals Using Resting-State fMRI Functional Connectivity Graph Community Identification

Klaudius Kalcher 1, 2 \* † , Roland N. Boubela1, 2 †, Wolfgang Huf 1, 2 , Christian Našel 3 and Ewald Moser 1, 2, 4

1 Center for Medical Physics and Biomedical Engineering, Medical University of Vienna, Vienna, Austria, 2 MR Centre of Excellence, Medical University of Vienna, Vienna, Austria, 3 Department of Radiology, Tulln Hospital, Karl Landsteiner University of Health Sciences, Tulln, Austria, 4 Brain Behaviour Laboratory, Department of Psychiatry, University of Pennsylvania Medical Center, Philadelphia, PA, USA

Edited by: Brian Caffo, Johns Hopkins University, USA

### Reviewed by:

Hidenao Fukuyama, Kyoto University, Japan Xin Di, New Jersey Institute of Technology, USA

\*Correspondence: Klaudius Kalcher<EMAIL_ADDRESS>
† These authors have contributed equally to this work.

### Specialty section:

This article was submitted to Brain Imaging Methods, a section of the journal Frontiers in Neuroscience

Received: 15 August 2015 Accepted: 25 November 2015 Published: 16 December 2015

### Citation:

Kalcher K, Boubela RN, Huf W, Našel C and Moser E (2015) Identification of Voxels Confounded by Venous Signals Using Resting-State fMRI Functional Connectivity Graph Community Identification. Front. Neurosci. 9:472. doi: 10.3389/fnins.2015.00472 Identifying venous voxels in fMRI datasets is important to increase the specificity of fMRI analyses to microvasculature in the vicinity of the neural processes triggering the BOLD response. This is, however, difficult to achieve in particular in typical studies where magnitude images of BOLD EPI are the only data available. In this study, voxelwise functional connectivity graphs were computed on minimally preprocessed low TR (333 ms) multiband resting-state fMRI data, using both high positive and negative correlations to define edges between nodes (voxels). A high correlation threshold for binarization ensures that most edges in the resulting sparse graph reflect the high coherence of signals in medium to large veins. Graph clustering based on the optimization of modularity was then employed to identify clusters of coherent voxels in this graph, and all clusters of 50 or more voxels were then interpreted as corresponding to medium to large veins. Indeed, a comparison with SWI reveals that 75.6±5.9% of voxels within these large clusters overlap with veins visible in the SWI image or lie outside the brain parenchyma. Some of the remaining differences between the two modalities can be explained by imperfect alignment or geometric distortions between the two images. Overall, the graph clustering based method for identifying venous voxels has a high specificity as well as the additional advantages of being computed in the same voxel grid as the fMRI dataset itself and not needing any additional data beyond what is usually acquired (and exported) in standard fMRI experiments.

Keywords: fMRI, BOLD, graph analysis, graph clustering, physiological signals, brain, veins

# 1. INTRODUCTION

Any interpretation of fMRI results as indirect measures of neuronal activation rests on the assumption that magnetization changes caused by changes in blood oxygenation are due to brain activity in the immediate vicinity. Whether and to what extent this assumption holds, however, has been the matter of much debate from the first days of fMRI onwards. While the discussion of this "brain or vein" question has generated a wealth of research on identifying veins and signals originating from them, the majority of fMRI studies still ignore the issue and take no measures to assess or reduce the influence of signals from major vessels (Menon, 2012).

Quantification of the influence of draining veins on fMRI results has provided evidence that it is most pronounced at low magnetic field strengths, and the relative influence of microvasculature to the MR signal increases as the field strength increases (Duong et al., 2003): while at 1.5 T the signal originates virtually entirely in the macrovasculature (Lai et al., 1993), the ratio shifts in favor of the microvasculature as a major signal source at 3T, 4T, and 7T. Still, even at higher field strengths, protocols using gradient echo EPI are highly sensitive to signal changes originating from larger veins, to the point that no significant improvement can be observed between 3T and 7T results (Geißler et al., 2013).

Several efforts have been made to reduce the influence of venous signals on fMRI measurements, and a number of effective ways for increasing the specificity of signals for the microvasculature have emerged from them. The use of spinecho instead of gradient echo sequences can drastically reduce extravascular signal contributions (Duong et al., 2003), but only at a steep cost in terms of signal-to-noise and contrast-to-noise ratio (Norris, 2012). Specific corrections to eliminate venous signals based on phase images have also been developed (Menon, 2002; Rowe and Logan, 2005; Curtis et al., 2014), but are rarely included in the sequences used by typical fMRI studies, and yield the risk of introducing errors through over-correction (Nencka and Rowe, 2007).

The development of these methods for identifying veins is to some extent related to that of methods for eliminating physiological influences, leading to converging developments. One of the most common approaches to address physiological signal contamination is RETROICOR (Glover et al., 2000), which uses externally measured respiration and cardiac signals for a regression-based correction of fMRI data. A later method termed CompCor (Behzadi et al., 2007) eliminates the need for externally measured signals by identifying potential regressors from either ventricular and white matter signals, or from signals in voxels with higher-than-average time course standard deviation—a feature typically seen in voxels containing larger veins. The most recent development, Highcor, merges this line of research with the work on phase-based venous suppression mentioned above by using the correlation between phase and magnitude image time series to identify venous voxels that can be used to extract regressors for physiological noise reduction (Curtis and Menon, 2014).

The reasons for identifying venous voxels are more complex than the elimination of global physiological noise, however, as any BOLD effects necessarily carry the potential for signal changes further downstream along draining veins. It is thus desirable to not only reduce global physiological noise as in these regression approaches, but also to reliably identify voxels with potential venous signal contributions to help interpreting signal changes seen in and around them. The most direct way of localizing macrovascular effects might be the creation of venous maps using susceptibility weighted imaging (SWI; Reichenbach et al., 2000; Haacke et al., 2004). Still, extravascular signal influences in fMRI measurements might extend beyond the delineation of veins in the SWI image, and imperfect coregistration of SWI to fMRI images can further limit the precision of this method of localization. An immediate measure proposed to minimize macrovascular influences is the elimination of voxels with time series coefficients of variation much larger than the local average of their surrounding voxels [a concept related to the second of the two approaches used by Behzadi et al. (2007) in CompCor, see above], which empirically corresponds to regions next to large vessels, as done in the minimal preprocessing pipeline (Glasser et al., 2013) of the Human Connectome Project (Van Essen et al., 2013). This method has the advantages of requiring neither separate protocols nor specialized measurement techniques, and estimates affected voxels directly from standard fMRI time series; however, there exists no evaluation of the correspondence of the time series standard deviation with venous effects.

The question of which voxels are influenced by large vessels is not a purely theoretical exercise. While in the beginnings of fMRI, the brain-vs.-vein debate was settled by a general acknowledgment that with the imaging resolutions then available, the practical relevance of knowing whether a signal originates in the microvasculature in the cortex or in the draining vein at its surface was limited. There are multiple reasons why this argument should no longer be used as an excuse for avoiding the question. First, many veins run in sulci and might cause a signal change whose causal origin is on the gyrus on one side of the sulcus be misattributed to the gyrus on the other side. While earlier imaging techniques might not have allowed to make such distinctions regardless of whether the measured signal originated in the parenchyma or in a draining vein, many current analysis techniques like surface projections for surfacebased analyses rely on correct attribution of signals at this level. The second, somewhat related, reason is that improvements in spatial resolution at higher field strengths and with the use of more sophisticated acceleration techniques (Simultaneous Image Refocused EPI, Multiband EPI) have led to the possibility of imaging at sub-millimeter resolutions (Feinberg et al., 2010), but this improvement in nominal spatial image resolution can only lead to interpretable gains if the physiological basis for the measured effect matches this granularity. Indeed, Turner (2002) suggested that due to dilution effects, draining vein effects might not be seen at more than a few millimeters distance from the gray matter region drained where the effect originated. But while a spatial gap between neuronal origin and the immediate source of the measured BOLD signal of 4 mm might be considered of limited relevance when imaging at a spatial resolution of 3 mm, the existence of unavoidable spatial discrepancies of this magnitude would render advances into higher resolutions entirely pointless. Finally, it is possible to show that at least in some cases, draining vein effects might occur at much larger distances from their neuronal origins, as is the case for signal changes in the basal vein of Rosenthal (BVR) next to the amygdala (Boubela et al., 2015).

In the analysis of the BVR signals, it also became apparent that one distinguishing feature of venous voxels wastheir resting-state functional connectivity pattern (Boubela et al., 2015), exhibiting very strong positive and negative correlations to other voxels in the macrovasculature. Thus, in an approach similar to functional parcellation methods of the brain (Eickhoff et al., 2015) this resting-state connectivity structure between voxels could be used to distinguish voxels in the macrovasculature from others. One such approach consists in using graph learning tools on the connectivity graph (where vertices correspond to voxels or sets of voxels, and edges link vertices with correlated time courses together). Previously, connectivity analyses have rarely been performed on a voxel-wise level, among others for computational reasons: if, for example, 150,000 voxels lie within the brain mask, the complete voxel-by-voxel correlation matrix would consist of 2.25 · 1010 entries, taking up about 167 GB (in practice, the number of voxels is typically reduced by restricting analysis to gray matter voxels or by resampling the data to a coarser resolution). Not all of these entries actually need to be stored to perform analyses, in particular when analyzing a relatively sparse connectivity graph, and efficient tools for tackling similar problems on large datasets have been developed in other fields.

In this work, graph based cluster analysis was performed to show how these tools can be applied to solve a practical problem of fMRI data analysis. Voxel-by-voxel correlations are computed for all in-brain voxels to create a voxelwise connectivity graph. Resampling to a coarser resolution as well as limiting analysis to a subset of voxels (e.g., gray matter voxels) are avoided as they would hamper the specific research question: resampling would lead to a loss of specificity in that affected voxels would be averaged with adjacent voxels to form the larger voxels of the coarser grid, and a gray matter mask might exclude parts of the venous structure, effectively hindering the identification of venous voxels based on their connectivity to voxels within other veins. Based on previous observations (Boubela et al., 2015), the largest clusters with the strongest (positive and negative) correlations among their voxel's time-series emerging from a clustering of this graph could be expected to reflect medium to large veins.

### 2. MATERIALS AND METHODS

### 2.1. Subjects

Fifteen healthy subjects (8 females/7 males, mean age 35.3, SD 13.3) were recruited at Medical University of Vienna. Exclusion criteria were prior psychiatric or neurologic illnesses, as well as the usual exclusion criteria for MR studies. All subjects gave written informed consent prior to the scan and the study was approved by the local institutional review board (Ethikkommission der Medizinischen Universität Wien).

### 2.2. Data Acquisition Protocols

All MRI scans were performed on a 3 Tesla TIM Trio using the standard 32-channel head coil and whole-body gradients (Siemens Medical Solutions, Erlangen, Germany). First, a highresolution anatomical image was acquired using MPRAGE with 1 × 1 × 1.1 mm3 resolution, and 160 sagittal slices (TE = 4.21 ms, TR = 2300 ms, flip angle 90◦ , inversion time 900 ms). Second, BOLD fluctuations at rest were measured with a short-TR multiband EPI-sequence (Feinberg et al., 2010) using 1.7 × 1.7 × 2 mm3 resolution, 2 mm slice gap (matrix size 128 × 128, 32 axial slices, TE = 31 ms, TR = 333 ms, flip angle 30◦ , multiband factor 8, bandwith = 1776 Hz/Pixel) collecting 1200 volumes. Finally, susceptibility weighted images (SWI) were acquired at 0.6 × 0.6 × 2.0 mm resolution (matrix size 384 × 384, 52 slices per slab, 1 slab, TE = 29 ms, TR = 42 ms, flip angle 15◦ ) to visualize medium to large venous vessels.

### 2.3. Preprocessing

To keep closely to the original images, only minimal preprocessing was applied to functional data, including only skull stripping using FSL BET, motion removal using FSL MCFLIRT, and band-pass filtering. For the latter, the pass-band used was 0.01–0.2 Hz, to avoid as far as possible influence from high-frequency respiratory or cardiac fluctuations. SWI images were segmented using FSL FAST for vein delineation, coregistered to the EPI weighted images, with the vein masks generated from segmentation also being transformed into EPI space using the resulting transformation parameters, using trilinear interpolation to ensure that all voxels in EPI space with some overlap with veins from the SWI mask have non-zero values. This vein map in EPI space was then binarized to generate a vein mask for the EPI images.

### 2.4. Graph Generation

Pairwise Pearson correlation coefficients were computed between all voxels within the brain mask, using GPUs for the calculation of the correlation coefficients (Boubela et al., under revision) and splitting the dataset into tiles to allow for the computations to fit within GPU memory (6 GB). For each subject, this correlation matrix was thresholded to generate the adjacency matrix of a graph using a correlation threshold such that S < 4, with

$$\mathcal{S} = \frac{\log E}{\log K}$$

(where E is the number of edges and K the average node degree). This thresholding criterion results in a rather sparse graph of only the strongest correlations, which are more likely to reflect adjacent voxels along the same vein or otherwise highly congruent voxel signals (as opposed to the more subtle long-distance connections of brain networks of neuronal origin; see below in the discussion for more details on the effect of the sparsity criterion). For each subject, the largest correlation threshold fulfilling S < 4 was computed iteratively by decreasing the threshold in steps of 0.01, starting at 1. The threshold was applied to the absolute value of the correlation coefficients to take into account both positive and negative correlations exceeding a certain correlation strength. The resulting connectivity graph had all in-brain voxels as vertices, and each correlation between two voxels that was above the threshold resulted in edges between the two corresponding vertices, with the correlation coefficient used as edge weight. Graphs were represented using the package igraph (version 1.0.1) in R (version 3.1.1). Self-loops and

multiple edges were eliminated using the igraph function simplify.

## 2.5. Graph Cluster Identification

Community identification on the graph was performed using the method based on modularity optimization by Newman (2006) as implemented in the igraph function cluster\_fast\_greedy. The optimization of graph modularity means that the resulting clusters are defined by their voxels having maximum connectedness among each other and minimum connectedness to voxels outside their own cluster. All voxels from all clusters that individually contained 50 or more voxels were then pooled into a single mask, which thus contained all voxels with time-courses strongly correlated (either positively or negatively) with those of a large number of different voxels. It should be noted that Newman's method is intended to detect communities in connected networks and that its application on sparse networks as used here might result in some cases in entire connected components being categorized as clusters. Nonetheless, for our purposes, this is still sufficient to detect groups of voxels with highest relative connectedness to each other considering the general sparsity of the graph.

# 2.6. Validation

To show that these voxels correspond mostly to vasculature, the overlap with the vein mask from SWI was computed (in EPI space); the proportion of voxels from the graph clustering map that overlaps with the SWI vein mask can be seen as a measure for the specificity of the graph clustering method, though it should be kept in mind that it is not the true specificity because the segmented SWI is not the ground truth for the identification of venous voxels: coregistration imperfections can lead to spatial deviations in the localization of these voxels, and not all low signal intensities in SWI originate from veins since other factors like iron levels (higher in the basal ganglia than in the rest of the brain) or proximity to air cavities or bone affect susceptibility. The latter observation also implies that it is impossible to make any meaningful quantification of the sensitivity of the graph clustering method by using SWI, as it means that an accurate map of venous voxels should not indiscriminately include all voxels with low signal intensity in SWI.

# 3. RESULTS

Overall, of the voxels within the brain masks (between 142,800 and 172,300 for the different subjects, mean 157,600, SD 10,540), 17,730 ± 5069 voxels (or 11.2 ± 3.0%) were identified by the graph clustering algorithm as being part of large highly coherent networks (see **Table 1**).

Single-subject images of the graph clustering masks overlaid over SWI are shown in **Figure 1**. The spatial distribution of the voxels identified by the method seems to exhibit a consistent pattern. Most of the voxels within the brain follow the path of veins visible in the SWI underlay. Another set of voxels delinates areas of low signal quality in orbitofrontal regions subject to susceptibility artifacts or at the edge of the brain, in either case such voxels could be discarded for fMRI analyses interested in neuronal effects.

This observation can also be quantified by comparing the mask gained from graph clustering with a binarized mask gained from segmenting the SWI image (see **Figure 2**), and the average proportion of voxels of the mask identified via graph clustering overlapping with veins in SWI is 0.67±0.05. A further significant proportion of voxels not directly overlapping veins lies on the edge of the brain mask, as identified by eroding the brain mask with a 5 × 5 × 5-voxel kernel (see **Figure 4**) using the R package mmand (Clayden, 2014), outside of what can be recognized as the brain itself (BET seems to be rather conservative in skullstripping), possibly reflecting signals from superficial veins, raising the overlap proportion to 0.77 ± 0.06 (see also **Table 1**).

In some brain areas (e.g., in the medial prefrontal region in **Figure 1**), the locations of the vein recognizable in SWI on one hand and the voxels of the graph clustering mask on the other can be observed not to overlap perfectly. Still, the similarity of the shape between the two features, only shifted by 1–2 voxels, strongly suggests that they are caused by the same underlying structure (i.e., the same vein). Such discrepancies are not necessarily worrying. Since the graph clustering mask is generated from the EPI voxel timecourses themselves as opposed to the SWI images acquired in a separate measurement, they can be seen as yielding potentially valuable complementary information on the effect of a vein on the EPI measurement.

Comparing the time series standard deviations of voxels within and outside the clustering brain mask reveals that voxels

The minimum, first quartile, median, mean, third quartile, and maximum across all subjects is given for the correlation threshold obtained for the binarization of the correlation matrix to a graph adjacency matrix, the absolute numbers of voxels within the brain mask and within the graph clustering mask, the percentage of voxels of the brain mask that were within the clustering mask, as well as for the percentage of voxels of the graph clustering map that lie within the mask generated from SWI segmentation (which can be interpreted as a measure of specificity), without and with the additional inclusion of brain edge voxels.

Figure 1). The underlay is the SWI, coregistered, and resampled to the EPI space. Voxels within the mask gained from the segmented SWI are yellow, voxels from the graph clustering mask that are within the SWI mask are red, voxels from the graph clustering mask that do not overlap with the SWI mask are blue. Note how the SWI segmentation mask tends to be rather unspecific to veins in regions with susceptibility-related low signal intensities of different origin, most notably in the basal ganglia due to their high iron levels.

within the mask indeed have on average significantly higher standard deviations (p < 2 · 10−16 in all individual subjects), in congruence with the underlying assumption of the method used by Glasser et al. (2013). However, the overlap between the signal standard deviations in macro- and microvasculature is very pronounced in all subjects, suggesting that a threshold based on the signal standard deviation alone might not be sufficient to discriminate between the two types of voxels (see **Figure 3**).

# 4. DISCUSSION

The results presented here show that graph based brain network analysis on a voxelwise basis can yield important insights on the origins of the underlying signals. In particular, network clustering yields a set of voxels defined by strong connections among each other and weak connectivity to voxels outside of this set that contains mostly voxels in medium to large veins (as identified by SWI) as well as extra-cerebral voxels. The very structure of these clusters, voxels spread across the whole brain with a large number of connections with correlation strengths above 0.7–0.95 (depending on the subject), strongly suggests that these voxels contain little information related to local neuronal activity, and rather reflect blood oxygenation changes on a larger (physiologic) scale. These voxels thus violate the underlying assumption of most fMRI analyses that BOLD signal changes in a voxel can be interpreted as an indirect measure of local neuronal activations, and should thus be excluded in this type of analysis.

FIGURE 3 | Comparison between the distributions of time series standard deviation in voxels outside (black) and within (red) the clustering mask over all subjects. Note that while time series standard deviations in voxels identified as veins tends to be higher than in those outside veins, the boudary is not clear-cut, highlighting the limits of venous voxel identification based on time series standard deviations.

It is worth noting that some correlation patterns typically arising in fMRI datasets are not reflected in the results computed here and do not seem to confound the correlation-based vein identification. The first of those is that in task fMRI, large vessels often correlate with task signals in active brain regions as they are draining the blood from those regions. One might thus expect that conversely, there should be a correlation between the venous signals investigated here and the signals from voxels within the parenchyma that are drained by these veins, making it difficult to distinguish venous voxels from parenchyma voxels. This does not seem to occur here, as indicated by the overlap of most voxels in the clustering mask with venous voxels identified by SWI. One reason for this might be that the correlation strength between the signals within the vein and each of the individual regions drained by the vein is substantially below the relatively high correlation thresholds that emerge from the high S threshold used in network creation. For medium to large veins, each gray matter voxel ultimately drained by them contributes only a small part to the signal in voxels within the vein, thus leading to lower correlation strengths between parenchyma voxels and veins of this scale in accordance with theoretical models of downstream dilution of effects in veins. This might explain why the clustering mask includes only larger vessels, and fails to identify some of the smaller vessels appearing in the SWI image. The reason for the absence of the correlations between signals from large veins and large areas of activation that typically occur in task fMRI might be that in task fMRI, there is an artificially high coherence of a particular (set of) brain region(s) with the signal in the vein due to the task-induced structure in these activations. In resting-state data as used in this study, however, the patterns along which all regions draining into a particular vein contribute to its signal are less coherent among each other, with different regions potentially contributing differently to the venous signal, and thus having lower individual correlations with it.

The second type of correlation pattern that might be expected to be visible in a functional connectivity graph are restingstate networks previously described in the literature, such as the default-mode network (DMN) or the left and right frontoparietal networks. The reason for them not appearing in the cluster mask is that in unblurred datasets as those used in this study to compute the correlation graphs, the correlations between parenchyma voxels in these networks are much lower than those between venous voxels, and the high correlation threshold used to construct the graphs ensures that only the latter are reflected in it. This subtlety of voxel time course correlation patterns is easily lost when using only blurred datasets, but can be visualized

effectively using the DMN as an example. **Figure 5** shows the functional connectivity of two voxels in the posterior cingulate cortex (PCC), a main constituent of the DMN, one of them a voxel clearly in a vein as identified by SWI (its functional connectivity map being shown in the top part of the figure), the other being an adjacent voxel outside the vein (its functional connectivity map is shown in the bottom part of the figure). For the venous voxel, the correlation coefficients in adjacent venous voxels exceed 0.63, which was the threshold for graph construction in that particular subject (based on the network sparcity criterion of S = 4; this was the lowest connectivity treshold for all subjects, see **Table 1**), as well as a number of other voxels in typical DMN regions in the medial prefrontal cortex as well as bilaterally in the parietal cortices. On closer inspection, one notices that all of those voxels can be related to veins identifiable on the SWI image (see zoomed inserts). Correlations with other voxels in DMN regions, including voxels further from large visible veins, also exist, but with far lower correlation strengths. When using a seed outside of veins, as exemplified by the connectivity map in the lower part of the figure of an adjacent voxel in the parenchyma, correlations above the threshold can be found neither in other DMN regions nor even in adjacent voxels.  This is consistent with results typically obtained from blurred resting-state datasets: signal time courses from veins draining one part of a resting-state network can be seen as reflecting the signal time courses of the voxels in the gray matter that these veins drain in a way similar to how the time course of a voxel in a blurred dataset reflects the time courses of the voxels in its neighborhood, and time courses of voxels in the draining veins of different parts of a restingstate network are more strongly correlated with each other in the same way as correlation strengths between voxels in different parts of a resting-state network are increased in spatially blurred datasets. Correlation strengths between parenchyma voxels of resting-state networks are lower, and thus, if the correlation threshold used to generate a binarized graph from the voxelwise correlation matrix is high enough, these parenchyma voxel correlations are not reflected in the graph analyses performed after binarization, and only the connections between venous voxels remain in the graph and, ultimately, define the graph modules.

Of course, this work is not the first attempt at identifying artifactual signals in fMRI datasets. Previous attempts include both the use of complementary measurements (including SWI) as well as the localization of its effects based on EPI measurements alone. The advantage of using EPI-based identification over complementary measurements is the more immediate relationship to the application of the results to the fMRI analyses in question. In the results presented here, this can be seen in the slight differences in localization of venous voxels between the clustering and SWI-based masks in the context of general correspondence between the two (see for example **Figure 1**). While the general similarity between the features seen in the SWI and clustering masks corroborates the theory of venous origins of the signals seen in the graph cluster voxels, the difference in location highlights that imperfect correspondence between two different modalities, emerging from imperfect coregistration, geometric distortion or other origins, necessarily limit the use of inference from complementary measurements for the identification of affected voxels in the EPI image. In this sense, vein identification methods using EPI and SWI complement each other as different, independent modalities, and both are important to provide a link between the results from EPI images to signal sources not directly visible in BOLD EPI such as most veins. While SWI yields the more anatomically accurate maps of the venous architecture in a subjects' brain, methods based on post-processing techniques applied on the BOLD EPI data themselves add to this a more direct view on the immediate effect of these veins on fMRI measurements.

Among methods based on EPI measurements, two categories can be distinguished, the first being methods using externally measured signals and correlating them with time courses from the EPI measurement, and the second being methods based on the analysis of EPI time courses by themselves. The first category typically uses high-frequency physiological nuisance signals (usually heart rate and respiration monitoring), which, however, has a different spatial distribution than the venous signals investigated here (Windischberger et al., 2002): highfrequency physiological noise tends to be concentrated near arteries and the CSF, which is subject to the same pulsations, while the low-frequency physiological signals investigated here tend to be localized in or near the venous macrovasculature. Physiological low-frequency signals are acquired and analyzed only in very few studies, but studies using them have shown them to be quite useful in identifying blood flow related phenomena in fMRI datasets (Tong and Frederick, 2012; Tong et al., 2014). The use of peripheral measurements, however, has one practical and one more fundamental limitation. The practical limitation is the necessity of additional hardware and measurement overhead for their acquisition leading to such measurements not always being available for all fMRI datasets, and the potential for additional error sources in their acquisition. While this issue can be overcome in any given study if the necessary steps are taken prospectively, it cannot be employed when analyzing datasets acquired without measuring these peripheral physiological signals, as is often the case in investigations using data shared by other researchers (Biswal et al., 2010; Kalcher et al., 2012). A more fundamental issue, though, is the time delay involved between the recording of the physiological signals at the external measurement location (e.g., the fingertip or toe for pulse oxymetry) and the brain, or, to be more precise, different locations in the brain. With standard EPI sequences as currently used, with a TR of between 2 and 3 s, the effect of this issue is rather limited, but with the current development toward short-TR multiband EPI sequences with higher temporal resolutions, the difficulties arising from the delay between the peripheral acquisition of physiological signals and their effect in the brain become more pronounced.

Finally, the use of measures directly derived from the EPI time series has been mostly confined to computationally less complex methods (e.g., the voxelwise time course standard deviation, a variant of which has been used in the Human Connectome Project), in part due to the lack of tools for tackling the computational challenges posed by more sophisticated methods like the voxel-by-voxel graph clustering approach presented here. Readily available tools from the domain of big data analysis can be applied to overcome computational obstacles and open the way to more comprehensive analysis tools. The comparison of time course standard deviations within and outside the graph clustering mask (see **Figure 3**) confirms the rationale behind the Human Connectome Project's preprocessing step of eliminating voxels with higher than normal standard deviations, but at the same time suggests that a one-dimensional measure not taking into account the connection structure between voxels might not yield a clear-cut discrimination threshold, as values of this score for normal brain tissue voxels with relatively high signal standard deviation and venous voxels with relatively low signal standard deviations overlap substantially.

In contrast, voxelwise graph analysis can be a useful tool to identify voxels in the macrovasculature by their highly correlated low-frequency signals. This latter point should be highlighted, as the band-pass filter applied (0.01–0.2 Hz) eliminates the possibility that the correlated signals in those voxels can merely be attributed to large-scale physiological noise (e.g., of respiratory or cardiac origin), which would have a higher frequency signal spectrum. Instead, they might exhibit more problematic signal fluctuations in the low-frequency domain, easily misattributed to local low-frequency fluctuations. In addition, the presence of such fluctuations might also be indicative of a risk of seeing downstream activations due to venous drainage of activations at more distant voxels, as it occurs during some emotional-visual tasks in the BVR (Boubela et al., 2015). The identification of voxels at risk is thus a powerful tool to increase specificity in the interpretation of fMRI BOLD activations.

### REFERENCES

### ACKNOWLEDGMENTS

This study was financially supported by the Austrian Science Fund (P22813, P23533).

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

Copyright © 2015 Kalcher, Boubela, Huf, Našel and Moser. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# The Voxel-Wise Functional Connectome Can Be Efficiently Derived from Co-activations in a Sparse Spatio-Temporal Point-Process

Enzo Tagliazucchi 1, 2, 3 \*, Michael Siniatchkin1 , Helmut Laufs 2, 4 and Dante R. Chialvo5, 6

1 Institute for Medical Psychology, Christian-Albrechts University, Kiel, Germany, 2 Department of Neurology and Brain Imaging Center, Goethe University Frankfurt am Main, Germany, 3 Department of Sleep and Cognition, Netherlands Institute for Neuroscience, Amsterdam, Netherlands, 4 Department of Neurology, University Hospital Schleswig-Holstein, Christian-Albrechts-University Kiel, Kiel, Germany, 5 Consejo Nacional de Investigaciones Científicas y Tecnológicas (CONICET), Buenos Aires, Argentina, 6 Center for Multidisciplinary Complex Systems Studies and Brain Sciences (CEMSC3), Escuela de Ciencia y Tecnología, Universidad Nacional de San Martín, Buenos Aires, Argentina

### Edited by:

Brian Caffo, Johns Hopkins University, USA

### Reviewed by:

Joshua T. Vogelstein, Johns Hopkins University, USA Jian Kang, Emory University, USA

\*Correspondence: Enzo Tagliazucchi<EMAIL_ADDRESS>
### Specialty section:

This article was submitted to Brain Imaging Methods, a section of the journal Frontiers in Neuroscience

Received: 24 June 2015 Accepted: 04 August 2016 Published: 23 August 2016

### Citation:

Tagliazucchi E, Siniatchkin M, Laufs H and Chialvo DR (2016) The Voxel-Wise Functional Connectome Can Be Efficiently Derived from Co-activations in a Sparse Spatio-Temporal Point-Process. Front. Neurosci. 10:381. doi: 10.3389/fnins.2016.00381 Large efforts are currently under way to systematically map functional connectivity between all pairs of millimeter-scale brain regions based on large neuroimaging databases. The exploratory unraveling of this "functional connectome" based on functional Magnetic Resonance Imaging (fMRI) can benefit from a better understanding of the contributors to resting state functional connectivity. In this work, we introduce a sparse representation of fMRI data in the form of a discrete point-process encoding high-amplitude events in the blood oxygenation level-dependent (BOLD) signal and we show it contains sufficient information for the estimation of functional connectivity between all pairs of voxels. We validate this method by replicating results obtained with standard whole-brain voxel-wise linear correlation matrices in two datasets. In the first one (n = 71), we study the changes in node strength (a measure of network centrality) during deep sleep. The second is a large database (n = 1147) of subjects in which we look at the age-related reorganization of the voxel-wise network of functional connections. In both cases it is shown that the proposed method compares well with standard techniques, despite requiring only data on the order of 1% of the original BOLD signal time series. Furthermore, we establish that the point-process approach does not reduce (and in one case increases) classification accuracy compared to standard linear correlations. Our results show how large fMRI datasets can be drastically simplified to include only the timings of large-amplitude events, while still allowing the recovery of all pair-wise interactions between voxels. The practical importance of this dimensionality reduction is manifest in the increasing number of collaborative efforts aiming to study large cohorts of healthy subjects as well as patients suffering from brain disease. Our method also suggests that the electrophysiological signals underlying the dynamics of fMRI time series consist of all-or-none temporally localized events, analogous to the avalanches of neural activity observed in recordings of local field potentials (LFP), an observation of potentially high neurobiological relevance.

Keywords: functional connectome, functional connectivity, dimensionality reduction, point process, resting state fMRI

# INTRODUCTION

The human brain comprises an interconnected network of cortical and sub-cortical regions globally linked by anatomical long-range tracts of connections. The mapping of the corresponding functional connections at a particular spatial scale (termed connectome in contemporary neuroscience; Sporns et al., 2004; Sporns, 2011) is an important ingredient in the process of understanding how the human brain can perform diverse cognitive functions. Furthermore, many neurological and psychiatric diseases can be understood in terms of deviations from a healthy connectome (Fox and Greicius, 2010; Kelly et al., 2012).

Advances in neuroimaging methods, such as Diffusion Tensor Imaging (DTI) and Diffusion Spectrum Imaging (DSI) allow the in vivo mapping of the human structural connectome at a largescale (Hagmann et al., 2008). Blood oxygenation level-dependent (BOLD) functional Magnetic Resonance Imaging (fMRI) allows for a functional counterpart of the anatomical connectome, a notion first introduced about a decade ago (Sporns et al., 2004; Eguiluz et al., 2005; Salvador et al., 2005) by computing the statistical covariance between all pairs of BOLD signals. This functional connectome contains information on how all pairs of regions (at a certain spatial scale) relate dynamically and collectively with each other.

These two approaches are being applied by international coordinated efforts to systematically map connectomes in very large populations of subjects and at the highest temporal and spatial resolution currently available (see for instance Biswal et al., 2010; Smith et al., 2013; Van Essen et al., 2013). These efforts will eventually lead to the availability of large-scale databases useful to account for potential inter-subject variability caused by different demographical variables, as well as to reduce the harmful effect of noise and artifacts through massive averaging.

These collaborative efforts need to be paralleled by methodological developments facilitating efficient extraction of relevant information from the data. Common strategies are based on averaging BOLD signals over brain parcellations comprising extended regions, thus reducing the dimensionality of the problem as well as the number of required computations. However, there are several problems inherent to this approach. First, all detail of the functional connectome inside each region of the parcellation is lost. Second, partitions are usually arbitrary and therefore might sub-divide a functionally coherent region into many regions. Different studies have addressed how the properties of parcellation-based networks can change depending on region selection (Wang et al., 2009; Zalesky et al., 2010). Third, efforts to increase the spatial resolution of fMRI sequences are pointless if data will be down-sampled after acquisition by averaging BOLD signals inside a small number of regions in a parcellation.

The objective of this paper is to show how a very sparse representation of brain activity, namely a discrete spatiotemporal point-process, is able to estimate the whole brain voxelwise functional connectome. This point-process is derived from the times at which the BOLD signals reach some maximum level of activity, either by detecting crossings of an arbitrary threshold, or by the identification of local peaks, i.e., the pointprocess comprises large amplitude events in the data. At its core, our proposed method is based on identifying a basis of discrete contributions to resting state fMRI signals, in analogy to other neural recording modalities (such as spikes in intraand extra- cellular recordings). Following this analogy, once the relevant events are identified, much of the signal (i.e., the stereotypical response associated with a discrete event) can be disregarded without reducing their information content, facilitating data storage, manipulation and interpretation (this analogy is limited, however, since neural recordings provide more sampling points than fMRI recordings and hence a larger number of discrete events). The main merit of this method is to reduce the continuous representation of BOLD signals into a series of timings associated with events of interest, thus (1) drastically reducing the dimensionality of the data, (2) abstracting the relevant information from sources of noise.

It has been shown previously that this method suffices to reproduce large-scale patterns of coordinated activity (Tagliazucchi et al., 2011, 2012a) termed Resting State Networks (RSN; Beckmann et al., 2005) and is essentially identical to the de-convolution of the signals as a series of discrete impulse functions (Petridou et al., 2013). Furthermore, de-convolution into a point-process can lessen the impact of hemodynamic lags for the estimation of causality between BOLD signals (Wu et al., 2013). Here we contribute a systematic evaluation of the capacity of this method to reproduce all bivariate relationships between signals (i.e., whole-brain correlation matrices). This validation is obtained, for the first time, from a large database of subjects n = 1147) scanned with different parameters at different locations, thus supporting its universal validity.

We also investigated whether abstracting the signal into a point-process could yield benefits from the perspective of reducing confounds and noise in the data. For this we adopted a practical, classification-based approach, investigating how accurately connectivity matrices derived from the point-process and from linear correlations could distinguish two groups of subjects (younger and older subjects from the n = 1147 database). We hypothesized that keeping high-amplitude events in the data could disregard low-amplitude noise and result in a better classification accuracy than the one obtained using full BOLD signals.

# MATERIALS AND METHODS

We will first describe all steps of the proposed method and then introduce different datasets used for validation as well as to show possible applications. The general procedure followed to estimate correlation networks via the point-process analysis is graphically outlined in **Figure 1**.

# Voxel-Wise Correlation Matrix

Consider an fMRI measure consisting of N voxels and T volumes, represented as Fn(t), with 0 ≤ n ≤ N and 0 ≤ t ≤ T. Thus, Fn(t) represents the BOLD signal at voxel n and time t. The common

definition of voxel-wise correlation matrix (Eguiluz et al., 2005) is as follows,

$$\mathcal{R}\_{\text{ij}} = \frac{<\{\mathcal{F}\_{\text{i}} - <\mathcal{F}\_{\text{i}} >\} (\mathcal{F}\_{\text{j}} - <\mathcal{F}\_{\text{j}} >)}{\sigma \left(\mathcal{F}\_{\text{i}}\right) \sigma \left(\mathcal{F}\_{\text{j}}\right)} \tag{1}$$

where < Fi > and σ(Fi) represent the mean value and the standard deviation of the BOLD signal at the voxel i, respectively. Note that according to this definition, for the computation of Rij, Equation (1) must be evaluated N(N − 1)/2 times (although not serially in efficient implementations). Often these calculations are used to define functional connectivity networks which in turn allow for further analysis of the resulting graphs.

### Constructing the Point-Process

The approach here proposed starts with converting the BOLD signal at every voxel into its z-score, F˜ i = Fi − <Fi> σ(Fi) . This is done under the assumption that, according to our formalism, the absolute amplitude of the BOLD signal carries less information than its temporal evolution (for the biological underpinnings of this assumption please see the Discussion section). To define the point-process, the a priori arbitrary threshold γ is selected and the spatio-temporal process PPi(t) is defined as follows:

$$\text{PP}\_{\mathbf{i}}(\mathbf{t}) = \begin{cases} 1 \text{ if } \mathbb{F}\_{\mathbf{i}} \left( \mathbf{t} \right) < \chi \text{ and } \mathbb{F}\_{\mathbf{i}} \left( \mathbf{t} + \mathbf{l} \right) > \chi \\\\ 0 \text{ otherwise} \end{cases} \tag{2}$$

This point-process was introduced in a previous publication (Tagliazucchi et al., 2012a) where we showed that it suffices to replicate the topographical features of the major canonical RSN, even though for most values of t and i, PPi (t) will be zero (indeed, taking γ = 1, for a signal of T = 240 on average the pointprocess is non-zero for 15 ± 3 time points, or approximately 6% of the data—see Tagliazucchi et al., 2012a). Note that once the point-process is constructed much of the data can be discarded. From a signal comprising 240 values, only a series of (on average) 15 numbers needs to be retained, namely, the timings of the events in the point-process. Clearly, this results in a considerable compression of the fMRI data.

Alternatively, PPi(t) can be defined by the (high amplitude) local peaks of the BOLD signal. For this, BOLD signals are also converted to z-scores and all sufficiently large peaks (for instance, those above an arbitrary threshold) are the points represented in PPi(t). The formal definition is as follows,

$$\text{PP}\_{\mathbf{i}}(\mathbf{t}) = \begin{cases} 1 \text{ if } \mathbb{\tilde{F}}\_{\mathbf{i}}(\mathbf{t}) > \mathbb{\tilde{F}}\_{\mathbf{i}}(\mathbf{t} - \mathbf{1}) \text{ and } \mathbb{\tilde{F}}\_{\mathbf{i}}(\mathbf{t}) > \mathbb{\tilde{F}}\_{\mathbf{i}}(\mathbf{t} + \mathbf{1}) \text{ and } \\\\ \mathbb{\tilde{F}}\_{\mathbf{i}}(\mathbf{t}) > \mathbb{\tilde{Y}} \\\\ 0 \text{ otherwise} \end{cases} \tag{3}$$

Although formally both methods are justified, it will be shown later that either definition of the point-process leads to similar results.

### Estimating Correlations from the Point-Process

After converting Fi into PPi (t) we introduce the following framework to generalize the methods introduced in Tagliazucchi et al. (2012a), from the estimation of seed based correlations to the efficient computation of all pairs of correlations between voxels. We first define the co-activation matrices Aij(t) as follows:

$$\mathbf{A}\_{\mathbf{i}\rangle}(\mathbf{t}) = \text{PP}\_{\mathbf{i}}(\mathbf{t})\,\text{PP}\_{\mathbf{j}}(\mathbf{t})\tag{4}$$

Note that according to this definition, Aij(t) only has two possible values: Aij (t) = 1 if at time t the point-process is non-zero both at voxels i and j, and Aij (t) = 0 otherwise.

The co-activation matrices defined in Equation (4) can be used to estimate the functional connectivity between all pairs of voxels in the brain by performing a simple matrix addition. Two highly synchronized signals will cross the threshold together most of the time, thus a measure of coupling between the signals can be obtained by counting the number of times the signals crossed the threshold together. This is formalized simply by,

$$\mathbf{C\_{i\circ}} = \sum\_{\mathbf{t}=1}^{\mathrm{T}} \mathbf{A\_{i\circ}}\left(\mathbf{t}\right) = \sum\_{\mathbf{t}=1}^{\mathrm{T}} \mathrm{PP\_{i}}\left(\mathbf{t}\right) \mathrm{PP\_{j}}\left(\mathbf{t}\right) \tag{5}$$

In matrix notation, this can be succinctly summarized as C = PP PPT , considering PP as a matrix with voxels as rows and time as columns and containing the point-process. The matrix Cij contains in its i, j entry the number of shared co-activations between BOLD signals at voxels i and j. Note that since all Aij are symmetrical matrices, then Cij is also symmetrical. Note also that the matrices Aij(t) contain valuable information about instantaneous co-activations between voxels and as such their analysis might be important to understand the temporal evolution of large-scale synchronization between brain regions (Tagliazucchi et al., 2012b; Hutchison et al., 2013).

The main issue with this matrix as a measure of functional connectivity is that it is not normalized, therefore there is no way to directly decide (for instance) if a perfect synchronization between signals has been reached. An appropriate normalization for this matrix would be as follows,

$$\tilde{\mathbf{C}}\_{\text{ij}} = \frac{\mathbf{C}\_{\text{ij}}}{\max\left(\sum\_{t=0}^{T} \mathbf{P} \mathbf{P}\_{\text{i}}, \sum\_{t=0}^{T} \mathbf{P} \mathbf{p}\_{\text{j}}\right)} = \frac{\mathbf{C}\_{\text{ij}}}{\max\left(\mathbf{C}\_{\text{ii}}, \mathbf{C}\_{\text{jj}}\right)} \tag{6}$$

This definition of C˜ ij is reasonable since Cij achieves its highest possible value if all threshold crossings are also shared between both voxels. However, one voxel could have all its threshold crossings in common with the other, whereas the opposite might not be true (since the other voxel could have a larger number of crossings in total (this can be the case only if Cii 6= Cjj), thus normalizing using the maximum between the number of crossings at both voxels is required. Also, C˜ ij is symmetrical with this normalization.

The normalization presented in Equation (6) requires the maximum value between the numbers of threshold crossings at all pairs of voxels. If normalization is needed, then a more efficient approximate solution is to divide by the number of threshold crossings without taking the maximum value, for instance, across rows or columns of the matrix, and then symmetrizing (if needed) the result by averaging with the transpose:

$$\mathbf{\dot{C}\_{\rm ij}} = \frac{1}{2} \left[ \frac{\mathbf{C}\_{\rm ij}}{\mathbf{C}\_{\rm ii}} + \frac{\mathbf{C}\_{\rm ji}}{\mathbf{C}\_{\rm ii}} \right] \tag{7}$$

Note that Cij Cii deviates from a symmetrical matrix only in the case of different numbers of threshold crossing between voxels (Cii 6= Cjj). Note also that normalization might not be necessary if comparing fixed-length recordings between two populations, under the assumption that the rate of events in the point-process is not different between groups.

For the computation of C˜ ij all steps can be performed efficiently in vectorized form in any language with matrix manipulation capabilities (for instance, MATLAB or Python with NumPy), in particular, after constructing the point-process in Equation (2), the operations involved consist of a single matrix multiplication (Equations 4 and 5), multiplication by scalars and matrix symmetrization (Equation 7). In this work, all computations were performed using a 8 core CPU running at 2400 MHz with a total of 128 GB built-in memory.

After introducing the core methods, we now discuss the methodology for the validation of our results.

### Measures Derived Whole Brain Voxel-Wise Correlations Used for Method Validation

The number of connections derived in a voxel-wise analysis complicates easy visualization of networks and their changes across conditions. Thus, in the many applications of functional connectomes found in the literature, rarely whole-brain voxelwise networks are directly visualized. Instead, lower-dimensional metrics are to be derived, which are easy to visualize as 3D maps overlaid on brain anatomy. One possible choice is to assess measures of network centrality, this is, how important nodes are in the network, thus collapsing all connections attached to a node into a single number. A straightforward definition in a weighted network is the strength (Barthelemy et al., 2005), defined as:

$$\mathbf{S}\_{\mathbf{i}} = \sum\_{\mathbf{j}=1}^{N} \mathbf{R}\_{\mathbf{i}\mathbf{j}} \tag{8}$$

In the present case, using the point-process to estimate correlations, Rij is replaced by C˜ ij. Nodes with the highest strength values are termed hubs and their reorganization has been repeatedly linked to different brain pathologies (Crossley et al., 2014), such as coma (Achard et al., 2012) or Alzheimer's disease (Buckner et al., 2009).

Note that the evaluation of Equation (8) requires the whole brain correlation network. In the case of a voxel-wise network, centrality of nodes (i.e., voxels) can be easily visualized as a 3D map overlaid on an anatomical image.

Another measure employed for validation of our method is the interhemispheric or homotopic connectivity. This is defined as the correlation between the BOLD signal of every voxel and the contralateral voxel. Interhemispheric connectivity is in particular useful to quantify re-organization of functional connectomes for which left-right asymmetries are expected (as in the case of aging, see Dolcos et al., 2002).

# Datasets

To demonstrate the validity of our proposal two different datasets from previously published studies will be used. The first dataset comprises BOLD fMRI recordings from the 1000 Functional Connectomes database, and the second dataset comprises recordings from a recently published study in which combined EEG, EMG, BOLD-fMRI, and physiological data were obtained from 71 subjects.

The Connectome dataset was downloaded from the 1000 Functional Connectome Project online database (http://fcon\_1000.projects.nitrc.org). Demographics, scanning parameters, and experimental conditions are described in the database website as well as in Tagliazucchi and Laufs (2014). Only epochs of wakefulness were employed in the present analysis. For more information on sleep vs. wakefulness classification in this dataset (see Tagliazucchi et al., 2012c; Tagliazucchi and Laufs, 2014). Since individual data presents variable length in this data set, normalization (Equation 7) was always required.

Data from a previously published study (Tagliazucchi and Laufs, 2014) was used for the sleep dataset. A total of 71 subjects were selected from a larger dataset on the basis of successful multimodal polysomnographic data recording and quality (written informed consent, approval by the local ethics committee). All subjects were scanned during the evening and instructed to close their eyes and lie still and relaxed. A group of 55 subjects was formed out of the original dataset of 71 subjects by excluding subjects who did not fall asleep. Hypnograms obtained via expert sleep staging based on AASM rules (American Academy of Sleep Medicine, 2007) were scanned for contiguous epochs of wakefulness, N1, N2, and N3 sleep lasting 250 volumes (∼ 2 min), resulting in 84 epochs of wakefulness, 16 epochs of N1 sleep, 19 epochs of N2 sleep, and 20 epochs of N3 sleep. Sleep epochs are present (by construction) fixed length in this data set (250 volumes), therefore normalization (Equation 7) was not required under the assumption that sleep does not modify the rate of points in the data.

EEG was recorded via a cap (modified BrainCapMR, Easycap, Herrsching, Germany) during fMRI acquisition (1505 volumes of T2∗ -weighted echo planar images, TR/TE = 2080/30 ms, matrix 64 × 64, voxel size 3 × 3 × 2 mm3 , distance factor 50%; FOV 192 mm2 ) at 3 T (Siemens Trio, Erlangen, Germany) with an optimized polysomnographic setting [chin and tibial EMG, ECG, EOG recorded bipolarly (sampling rate 5 kHz, low pass filter 1 kHz), 30 EEG channels recorded with FCz as the reference (sampling rate 5 kHz, low pass filter 250 Hz), and pulse oxymetry, respiration recorded via sensors from the Trio (sampling rate 50 Hz)] and MR scanner compatible devices (BrainAmp MR+, BrainAmp ExG; Brain Products, Gilching, Germany).

MRI and pulse artifact correction were performed based on the average artifact subtraction (AAS) method (Allen et al., 1998) as implemented in Vision Analyzer2 (Brain Products, Germany) followed by objective (CBC parameters, Vision Analyzer) ICAbased rejection of residual artifact-laden components after AAS resulting in EEG with a sampling rate of 250 Hz. Good quality EEG was obtained, which allowed sleep staging by an expert according to the AASM criteria (American Academy of Sleep Medicine, 2007).

# fMRI Preprocessing

Using Statistical Parametric Mapping (SPM8) EPI data were realigned, normalized (MNI space) and spatially smoothed (Gaussian kernel, 8 mm3 full width at half maximum). The data were band-pass filtered in the range 0.01–0.1 Hz using a sixth order Butterworth filter. The same procedure was applied to the sleep dataset and to the 1000 Functional Connectomes dataset.

# Multivariate Classification

We compared the accuracy of a Random Forest classifier with 100 estimators (implemented in scikit-learn, http://scikit-learn.org/ stable/) to distinguish younger (<20 years) and older (>40 years) subjects from the 1000 Functional Connectomes dataset. This was based both on strength and interhemispheric connectivity maps obtained via normalized co-activation matrices (derived from the point-process) and standard linear correlation matrices. We applied a 5-fold cross validation procedure combined with feature selection (F-test to retain the top 10, 25, 50, 75% features), as well as with all features. Accuracy was reported as the area under the receiver operator characteristic (ROC) curve (AUC).

# RESULTS

### Correlations between C˜ ij and Rij

We obtained the point-process for both datasets following the procedure illustrated in **Figure 1** and in the methods section. In the case of the 1000 Functional Connectomes dataset we repeated calculations both for voxel-wise networks and for networks based on time series extracted from the AAL template. Using this data, we first evaluated the similitude in the estimation of the connectivity matrix by both methods (point-process analysis with normalization and linear correlations) as a function of the threshold γ used to define the point-process (see Equation 2). Results are shown in **Figure 2** (left) for the average correlation between connectivity networks estimated by both methods as a function of γ. Correlations peaked at 0.6 and were highest for ≈ 0.7. The histogram of all 1147 correlations obtained using γ = 1 (**Figure 2**, center) revealed a sharp peak around the mean value. The plot of the entries of the estimated correlation (values of C˜ ij) and the linear correlation (entries of Rij) is shown in **Figure 2** (right). A monotonously increasing relationship was present between both quantities, even though the functional dependency between them was not linear. For low linear correlation values, the point-process co-activation increased slowly and did so more quickly for larger linear correlation values.

We compared the performance of computing voxel-wise functional connectivity matrices using the proposed pointprocess based method vs. standard linear correlations. In **Figure 2B**, left, the percentage of the time required using linear correlations (corrcoef.m MATLAB function, average time 131.48 s on a reference system) was plotted as a function of the threshold. At every threshold value a total of 100 iterations were performed for a single subject and results were then averaged. For

FIGURE 2 | (A) Left: Correlation between Rij and C˜ ij as a function of the threshold [γ in Equation (2); mean ± SEM]. Connectivity networks were derived from 116 time series extracted from the AAL template in all subjects from the 1000 Functional Connectomes dataset (n = 1147). Center: Histogram of all correlation values at γ = 1, P = probability. Right: Average (mean ± SEM) plot of the linear correlation coefficient between brain regions (entries of Rij) and the estimate from the point-process analysis (entries of C˜ ij). The inset shows the plot for each one of the 1147 subjects. (B) Left: Performance of the point-process based estimation of functional connectivity as a function of the threshold γ (mean ± SD). Elapsed computation times were obtained for a single subject across 100 repetitions and compared with the performance using linear correlations. Right: Percentage of the original number of data points retained after converting the data to a sparse point-process with γ = 1, plotted as a function of the threshold (for all subjects in the 1000 Functional Connectomes dataset). (C) Left: Cumulative computation time required to compute whole-brain voxel-wise connectivity matrices from 1000 subjects extracted from the Functional Connectomes dataset. An un-normalized point-process with γ = 1 was used. Right: Cumulative space required to store 1000 subjects from the Functional Connectomes dataset, both for the full data and for a sparse representation based on a point-process with γ = 1.

thresholds larger than approximately 1 standard deviation, the point-process based method slightly outperformed the standard computation, with performance becoming increasingly better as the threshold was increased and less points were included in the analysis. However, more evidence needs to be gathered to confirm that the method outperforms the standard linear correlation approach, considering that the routines have not been properly optimized. In **Figure 2B** (right) we plot the percentage of data points retained after conversion to the point-process. Even for the smallest threshold values, only about 6% of the data was retained. Thus, this very sparse representation of fMRI data contained sufficient information to capture all the differences during deep sleep and in the 1000 Functional Connectomes dataset (see below), requiring but a small fraction of the original time series. Specifically, the required information consists of the (discrete) timings of the events in the point-process (i.e., at which volumes the "points" appear).

To gauge the usefulness of our approach in a real setting, we computed the cumulative time and space required to process (i.e., obtain whole-brain voxel-wise connectivity matrices) and store 1000 subjects extracted from the Functional Connectomes dataset. Results are shown in **Figure 2C**. An un-normalized point-process with a threshold of γ = 1 resulted in a reduction of computation time (reference system) from a total of ≈30 h to ≈19 h. However, we note again that more careful experiments need to be performed to compare the time performance of both methods.

We also investigated the sparseness (defined as the percentage of zeros) in the point-process time series and in the associated normalized connectivity matrices (derived via point-process co-activations). The results are shown in **Figure 3A**. Not only the time series are very sparse (≈95% zeros for a threshold of 1 S.D.) but also the connectivity matrices (≈50% zeros for the same threshold). This results in dramatically smaller file sizes when both the time series and the connectivity matrices are stored (**Figure 3B**).

# Strength Maps in Wakefulness vs. Deep Sleep

To compare results obtained by both methods, we applied them to derive the strength maps (Equation 8) from the estimated whole brain voxel-wise correlations in the sleep dataset and to reveal changes between wakefulness and deep sleep. A total of 20 2-min epochs of deep sleep and 84 epochs of wakefulness could be extracted. After deriving the correlation networks, Equation (8) was applied to obtain the voxel-wise spatial distribution of strengths. Results for the contrast wakefulness > deep sleep are shown in **Figure 4A**, both for normalized and un-normalized co-activation matrices, as well as for the point-process derived from BOLD signal peaks instead of threshold crossings. Spatial patterns of decreased strength in deep sleep (comprising frontal, cingulate, primary visual, motor, and auditory cortices) were

(B) Cumulative file size (in bytes) of fMRI time series (left) and pair-wise connectivity matrices (right) derived using linear correlations (from the full data) and co-activations (from the point-process with threshold equal to 1).

captured equally well by both methods, as well as by the peakbased point-process. In particular, since fixed epoch lengths were used (250 volumes) results were reproduced with and without normalization of connectivity matrices as derived from the pointprocess. This similitude can also be seen in **Figure 4B**, in which a joint 3D rendering of both maps shows their spatial agreement. The main plots in **Figure 4C** show node strength values at all voxels computed using the point-process method (entries of C˜ ij) vs. those computed using linear correlations (entries of Rij). The functional dependency was clearly monotonously increasing on average, both for wakefulness and sleep, although two individual epochs of sleep displayed an opposite trend.

### Strength Maps in Young vs. Older Subjects

We then studied changes in node strength in the 1000 Functional Connectomes dataset, in particular, we compared a group of subjects younger than 20 years with an older group of subjects older than 40 years. Results can be found in **Figure 5A**. For both methods an increase of (normalized) functional connectivity strength in the older group was observed, comprising a network of regions that included the right parietal cortex, inferior frontal cortex, insula, and the precentral and postcentral gyrus.

Driven by the asymmetry observed in the strength differences between age groups, and by the proposal that the right hemisphere shows accelerated functional decline with aging (Dolcos et al., 2002), we applied linear correlations and the pointprocess analysis to quantify interhemispheric or homotopic connectivity between groups and compared the respective values. Results are shown in **Figure 5B**. Increased interhemispheric connectivity was observed for the older group of subjects by both methods, comprising areas in the parietal and temporal cortex, as well as in the precentral gyrus.

Finally, an additional calculation was performed to allow for further evaluation of our method. We regressed subject age vs. strength values in two regions of interest extracted from the analysis of young vs. older subjects (right Inferior Parietal Cortex—IPC, right and insular cortex). Strength values were obtained both from connectivity matrices obtained with linear correlations and with the point-process. Results are shown in **Figure 5C**. The plots show a moderate increase in strength with age, which suddenly increased for more mature subjects (age > 40 years approximately). Spearman's rank correlation coefficients were higher for the strength values computed using the pointprocess.

of age (in years) vs. strength values (derived from linear correlations and the normalized point-process) extracted from two regions of interest (right Inferior Parietal

Cortex—IPC, and right insular cortex; mean ± SEM). An almost monotonous (but clearly non-linear) relationship between age and network centrality is observed.

# Classification of Young vs. Older Subjects

We implemented the classifier described in the methods to investigate how accurately subjects could be classified by age using strength and interhemispheric connectivity maps, computed with both linear correlations and normalized point-process co-activations. Results are presented in **Figure 6**. We observed similar classification accuracy for the computation based on inter-hemispheric connectivity, and higher classification accuracy for point-process co-activations vs. linear correlations for the computation based on strength maps.

# DISCUSSION

We are witnessing in recent times how neuroscience, and in particular neuroimaging, is moving at a fast pace toward the accumulation and analysis of very large volumes of data. A

number of international collaborations is aiming to break new ground in the scale and speed of data collection, including the 1000 Functional Connectomes Project (Biswal et al., 2010), the NIH BRAIN Initiative (Insel et al., 2013), as well as the Human Connectome Project (Van Essen et al., 2013). These studies span hundreds of subjects scanned at high temporal resolution, resulting in very large datasets. Exploratory analyses of this data may thus benefit from biologically principled dimensionality reduction.

While it is obvious that having large volumes of data reduces the negative effect of noise, artifacts and the relative importance of the mathematical models employed to analyze it [a position eloquently defended by Halevy et al. (2009) in their seminal article "The Unreasonable Effectiveness of Data"], it is also true that the handling of redundant data might may be inefficient, both from a computational perspective and in terms of distinguishing the real contributors to the signal from sources of noise. In this line of thought, we have shown that the introduction of a sparse representation of fMRI datasets can reproduce findings obtained from full time series while keeping on the order of 1% of the original data. With respect to vulnerability to noise, sudden head movements can induce spurious points in the process, however, these can be identified from the realignment parameters and erased, following the strategy of scan censoring (Siegel et al., 2014) but eliminating single points (instead of continuous segments of data) from the analyses (see Tagliazucchi et al., 2014 for an application). A consequence of defining the point-process based on high amplitude excursions of the signal is that the impact of physiological noise sources affecting low amplitude fluctuations (Cordes et al., 2002) will be lessened.

# Sleep Validation Dataset: Loss of Connectivity in the Thalamus, Frontal, Midline, and Auditory Cortices

We validated our method by first computing correlation between connectivity matrices as obtained by both methods over > 1000 subjects in the Functional Connectomes dataset, as well as by comparing voxel-wise network strength (a measure of centrality computed from the voxel-wise network of functional connections) between wakefulness and deep sleep and between two age groups extracted from the 1000 Functional Connectomes dataset. In this latter dataset we also obtained the distribution of voxel-wise inter-hemispheric connectivity. The maps of altered network strength in deep sleep and the age-dependent effect observed in the 1000 Functional Connectomes dataset are of biological relevance themselves, as we are not aware of prior reports of these results. Deep sleep resulted in a loss of connectivity across all voxels located in frontal and cingulate cortices, as well as in the primary auditory cortex (Heschl's gyrus) and the thalamus. These are plausible correlates of reduced awareness (frontal and cingulate cortex) and loss of sensory engagement with the environment (primary auditory cortex and thalamus) resulting in increased arousal thresholds (Tagliazucchi et al., 2013).

# Age Groups Validation Dataset: Increased Connectivity with Age in Inferior Parietal and (Pre-)Frontal Cortices

With respect to the two different age groups extracted from the 1000 Functional Connectomes database, regions central to working memory processes (inferior parietal and frontal cortices, prefrontal cortex) showed "over-connectivity" in the older group of subjects. The meaning of this result is less clear, especially in the light of reports showing an inverse relationship between seed-based functional connectivity and age (Sambataro et al., 2010). However, voxel-based strength maps do not require any a priori anatomical hypotheses (i.e., seed selection) and thus might be capable of capturing more global changes in connectivity as opposed to the aforementioned approach. Interestingly, changes in the node strength values were mostly located in the right hemisphere. It has been noted by Dolcos et al. (2002) that the right hemisphere shows a more marked decline with aging, a fact supported so far by evidence from working memory neuroimaging experiments. The changes observed by the authors were hypothesized to be of compensatory origin, which is compatible with the outcome of our analyses (increased overall connectivity in the right hemisphere of older subjects). Prompted by this observation, we also found differences in interhemispheric connectivity located in a set of regions overlapping with those involved with changes in node strength.

# Why Few Points Are Sufficient to Reproduce Functional Connectomes

It is worthwhile discussing the reasons underlying the effectiveness of our approach, since it might be surprising that a small fraction of the data suffices to capture all bivariate relationships between BOLD signals (functional connectome) without sacrificing (and even enhancing) classification accuracy.

From a signal processing perspective the answer is relatively straightforward: keeping large amplitude events can increase the signal-to-noise ratio, since it discards low-amplitude activity containing a larger noise component. This non-linear filtering selectively amplifies the importance of those time points at which the signal amplitude becomes relatively large and therefore the signal-to-noise ratio increases. Physiological artifacts have been shown to affect BOLD signals at low frequencies and low amplitudes (Cordes et al., 2002) and signals measured in white matter and cerebrospinal fluid (which do not reflect activity of neural origin and are commonly employed as proxies for physiological confound time series) present smaller amplitude fluctuations compared to those in gray matter (see for instance Tagliazucchi et al., 2013). This situation can result in selective down-weighting of physiological noise when only large-amplitude excursions of the signals are considered.

From a biological point of view, the challenge is to understand why the fMRI time series can be effectively represented as a train of discrete impulses, a view of BOLD time series also supported by studies performing blind de-convolution of spontaneous activity (Petridou et al., 2013). Electrophysiological experiments reveal that Local Field Potentials (LFP) are spatiotemporally distributed as power law avalanches (Beggs and Plenz, 2003): most frequently, spontaneous LFP increases span a limited spatial area, however, at certain (discrete) points in time, LFP might extend up to the size of the tissue under study (an event termed avalanche). If LFP avalanches are, indeed, distributed following a scale-free power law, then macroscopic events (i.e., in the centimeter scale) should be observed, which would be sufficient to elicit a measurable hemodynamic response (considering the correlation observed between LFP and BOLD signals, see Logothetis et al., 2001). Indeed, spatio-temporal avalanches of activity can also be observed with fMRI, following the same statistical laws as the electrophysiological avalanches (Tagliazucchi et al., 2012a). Large amplitude macroscopic LFP increases were reported in the monkey cortex (Thiagarajan et al., 2010) and termed coherence potentials. These largescale events are also stereotypical (in the words of the authors, much like action potentials at the single-cell level) and thus fulfill all the theoretical requirements for the electrophysiological underpinnings of the events in the spatio-temporal fMRI pointprocess.

# Contributors to the Resting State fMRI Signal

One of the main limitations of fMRI compared to other noninvasive neuroimaging techniques (EEG, MEG) is its limited temporal resolution. This limitation not only stems from the relatively slow acquisition of whole-brain volumes (i.e., long TRs, in the order of seconds) but also from the coupling between neural activity and the signal measured by fMRI. This coupling blurs temporally localized activity into a temporally extended response (given by the HRF). Therefore, improvement in fMRI sampling rates will only result in a better-sampled HRF, with no gain in the measurement of underlying neural activity, unless the distortion caused by the HRF can be inverted.

Our results suggest that the fMRI resting state signal comprises a temporal succession of well-localized events. The identification of these events has been shown to match a formal de-convolution of fMRI time series (Tagliazucchi et al., 2012a; Petridou et al., 2013). This inversion of the HRF blurring can allow to capitalize on improvements in fMRI acquisition rates. While the contributors to the task-evoked fMRI signal have been thoroughly investigated, this remains to be done in the context of spontaneous brain activity; the possibility of reducing resting state fMRI signals to a few high-amplitude events and still estimate all pair-wise interactions represents an important first step in this direction, and suggests a focus for future studies on the electrophysiological basis of spontaneous fMRI fluctuations.

# Caveats and Limitations

Generally, this procedure should yield equivalent results for any dataset in which high amplitude events do not arise spuriously as artifacts and represent important information in the data. From a neurophysiological perspective, the fulfillment of these conditions has been already demonstrated for BOLD time series by means of inverting the Hemodynamic Response Function (HRF) convolution of neuronal sources (de-convolution). As discussed in the previous sections, LFP giving rise to metabolic changes reflected in the BOLD signal are temporally cluttered into avalanches of activity (Beggs and Plenz, 2003; Tagliazucchi et al., 2012a; Shriki et al., 2013), presumably underlying the high information content of BOLD signal high amplitude events.