“More Is Different” in Functional Magnetic Resonance Imaging: A Review of Recent Data Analysis Techniques

Correlation-based measures

Pearson's linear correlation is a bivariate linear measure of association with values ranging between +1 and −1. Linear correlations between fMRI time series are sometimes used synonymously with the term “functional connectivity.” In fact, the majority of seed-based functional connectivity studies are based on linear correlations. The basic strategy is to start with a user-specified seed voxel or region of interest (ROI) and test time courses in the remaining voxels for correlations with the seed. This strategy was successfully used in a number of studies, perhaps one of the most famous example being the study published by Biswal and colleagues (1995). Many other publications too numerous to list have followed since then.

Note that linear correlation can only detect linear dependencies between time series. Therefore, uncorrelatedness is not equivalent to independence because of the potential presence of nonlinear dependencies. Equivalence only holds in the special case when the two time series are jointly normal. Therefore, one should take care to avoid false conclusions about independence when linear correlations were used.

Spearman's rank correlation measures the linear correlation between ranks of values rather than the values themselves. Therefore, this measure is slightly less prone to errors of this type. Partial correlation can be used to measure the degree of linear association between two time series with the effect of other influences regressed out. Smith and colleagues (2011) found that partial correlation gave high sensitivity and performed comparatively well.

Frequency-based measures

As noted by Salvador and colleagues (2005), interregional dependencies can often be more readily observed in the frequency domain than in the time domain. The two most commonly used measures of association in the frequency domain are spectral coherence and wavelet coherence. In both cases, two time courses are considered to be coherent if their power spectra correlate for a given frequency window. Both measures ignore phase shifts, which makes them insensitive to a fixed lag between two time series.

Their main difference is that spectral coherence assumes stationarity of the signals, whereas wavelet coherence does not. Therefore, spectral coherence is insensitive to changes in the coupling between the two time series. It was used for fMRI data analysis in Salvador and colleagues (2005), Müller and colleagues (2003, 2007), and Lohmann and colleagues (2010).

Wavelet coherence was used in fMRI in many studies (e.g., Achard et al., 2006). Chang and Glover (2010) used it to investigate the nonstationary effects of the coupling. Some care must be taken because much of the power in fMRI signals resides in very slow fluctuations around 1/30 Hz. In other words, only about two of these very slow fluctuations occur within 1 min of measurement time. Therefore, to be able to reliably observe effects of nonstationarity, the time series must be long enough to cover several occurrences of these very slow waves.

Granger causality

Perhaps, the best known measure of directed connectivity is Granger causality (GC). A time series X is said to Granger cause some other time series Y if knowledge of past values of X helps to predict future values of Y in a statistically significant way. It is assumed that the data can be adequately described by a linear autoregressive model of some user-specified order p, where p generally ranges between one and four. It was introduced to the fMRI community by Roebroeck and colleagues (2005). For a recent review and recent developments regarding GC, see for example, Deshpande et al., 2010.

The application of GC to fMRI data has been controversial. Smith and colleagues (2011) investigated GC in the context of the simulation study mentioned above. They found that measures based on temporal lags such as GC performed very poorly. Their main point of critique was that the variability of the BOLD response makes the interpretation of temporal effects prone to error. Inter-regional differences in the hemodynamic response function may, in fact, be due to differences in the vasculature rather than to neuronal effects. For an experimental evaluation of the variability of the BOLD response, see e.g., Neumann et al., 2003.

However, rebutting this critique, Schippers and colleagues (2011) claimed that GC fared not so badly when used at the group level, rather than at the individual subject level. They reported that the directional information found by GC was “accurate in the vast majority of the cases.” Smith and colleagues (2012a) responded by noting that there may be a systematic bias not only at the single-subject level, but also at the group level, so that the conclusion drawn by Schippers and colleagues should be taken with care.

This debate was recently taken up by Seth and colleagues (2013). They show that even though GC is surprisingly robust to a wide variety of changes in hemodynamic response properties, it is still prone to error in the presence of heavy measurement noise or severe downsampling. Since both effects are typical of fMRI data, it would be inappropriate to give an all-clear signal for GC.

A possible solution to this problem was proposed by Deshpande and colleagues (2009). They note that a version of GC called “directed transfer function (DTF)” may be appropriate for fMRI data provided it is used in a way that is compatible with the true temporal resolution of the hemodynamic response, namely, at a scale of 6–12 sec. Using a this much reduced temporal resolution, they found that the primary motor area had a causal influence on the supplementary motor area, premotor area, and cerebellum.

Another problem inherent in bivariate directed measures such as pairwise GC is that misleading information may result if the system at hand has interdependencies among more than two nodes (see e.g., Blinowska et al., 2004). In this case, bivariate measures are often inadequate and should be replaced by multivariate methods. Figure 2 shows an example where a trivariate model would be required to disambiguate the two cases.

OPEN IN VIEWER

FIG. 2.

These two cases cannot be disambiguated using bivariate Granger causality.

Note, however, that the number of coefficients to be estimated for a multivariate autoregressive linear model is quadratic in the number of nodes (Harrison et al., 2003). For instance, if the model order p is two and the number of nodes d is five, then the number of coefficients to be estimated is k=p×d×d=50. In the case of fMRI, this may be problematic. The number of data points required for a reliable estimation of a multivariate autoregressive model should be many times larger than k, so that very long time series are needed. In fMRI, time courses generally only contain a few hundred data points so that this requirement may be hard to meet. Nonetheless, Deshpande and colleagues (2011) recently presented a multivariate GC analysis using 33 ROIs, and Tang and colleagues (2012) applied a highly multivariate GC to every voxel within two ROI.

Measures based on information theory

Principal component analysis

Principal component analysis (PCA) is a classical approach to matrix factorization. It was originally invented by Pearson (1901). It is used to obtain maximally decorrelated components.

Mathematically, it can be described as follows. Let V be a real-valued p×q matrix whose rows are measurement vectors, for example, fMRI time series. The columns of V must be centered, that is, their means must be subtracted. Let C be the p×p covariance matrix of V, and E the matrix of eigenvectors of C. Since E is orthogonal, EE^T equals the identify matrix I, and hence \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}V = EE^TV\end{align*} \end{document}

Let the eigenvectors in E be sorted such that the leftmost columns belong to the largest eigenvalues, and let W denote the matrix consisting of the r leftmost columns of E. Then, an approximation to the matrix V can be obtained by \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}V \approx WW^TV\end{align*} \end{document}

Thus, with H=W^TV, we have a matrix decomposition of the form V≈W H. The first principal components contain most of the variation in the data.

PCA was first applied to fMRI data by Friston and colleagues (1993) who aimed at identifying spatially decorrelated activity patterns. Nowadays, PCA is often used as a data reduction step to be followed by other methods such as ICA.

Independent component analysis

While the aim of PCA is to obtain maximally decorrelated components, ICA is intended to produce maximally independent components (Bell and Sejnowski, 1995). In this sense, ICA can be seen as a special case of blind source separation (BSS). Since very good reviews about ICA exist, we keep this section short and refer the reader to other references (e.g., Calhoun et al., 2009; Kiviniemi et al., 2003; Margulies et al., 2010).

In ICA, we assume that the data can be represented as a mixture of independent sources so that the data matrix V can be approximated as V≈W S, where W is a mixing matrix usually assumed to be square and invertible, and S is a matrix of source signals that we try to recover. The sources can then be reconstructed by inverting the mixing matrix so that the final outcome of an ICA decomposition is \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}S \approx W^{- 1}V\end{align*} \end{document}

where W ⁻¹ is sometimes referred to as the “unmixing matrix.” Many procedures for estimating the unmixing matrix exist (see e.g., Hyvärinen and Oja, 2000).

ICA was introduced into fMRI data analysis by McKeown and colleagues (1998), and later became one of the principal tools for analyzing resting-state fMRI data (Beckmann et al., 2005). In a highly cited study, Smith and colleagues (2009) reported that ICA components at the resting brain match up well with patterns of activations found in the BrainMap database of functional imaging studies (Laird et al., 2005). They concluded that networks in the brain are continuously active even when at rest.

Methodological extensions to the original ICA were, for instance, proposed in Beckmann and Smith (2005), who described a tensorial ICA that allows investigating groups of subjects rather than single-subject data. Most recently, decompositions into the temporal rather than spatial components were discussed in Smith and colleagues (2012b).

Both PCA and ICA require specification of a model order, that is, the number of components to be expected. Since the intrinsic dimensionality of the data is unknown, the selection of an appropriate model order is a difficult problem. For more details, see for example, Abou-Elseoud and colleagues (2010), Li and colleagues (2007), and Yourganov and colleagues (2011).

In addition to the difficulty of specifying model order, the selection and interpretation of meaningful components may often be quite problematic. Furthermore, Daubechies and colleagues (2009) have noted that two ICA algorithms most commonly used for fMRI (InfoMax and FastICA) seem to select for sparsity rather than for independence (Daubechies et al., 2009). In other words, ICA tends to lead to decompositions in which components are sparse in the sense that many of their voxels are close to zero [see Daubechies et al. (2009) for a discussion of the terms “independence/sparsity/separability”]. As a consequence, the authors suggest developing decomposition methods that directly target sparsity instead of independence. In the following paragraph, we introduce an alternative decomposition technique, which may perhaps be useful in this regard.

Non-negative matrix factorization

An alternative to ICA is non-negative matrix factorization (NMF) (Lee and Seung, 1999; Paatero and Tapper, 1994). In contrast to ICA, the components obtained by NMF are not designed to be independent. Rather, the aim is to obtain a decomposition into additive parts. NMF produces a decomposition of the n×m data matrix V such that \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}V \approx WH\end{align*} \end{document}

with W an n×r -matrix, and H an r×m-matrix. The r columns of W are called basis functions, and the columns of H represent encodings. All entries of the matrices V, W, H are non-negative so that the decomposition can be viewed as a weighted positive summation of factors that produce V. This distinguishes NMF from other decomposition techniques such as PCA or ICA where negative values are allowed. From their experimental results, Lee and Seung (1999) concluded that NMF seems to yield naturally sparse decompositions. Later studies, however, suggested that sparsity may have been an accidental byproduct in this early study, and other NMF algorithms have thus been proposed which directly target sparsity (e.g., Kim and Park, 2007).

NMF can be seen as a form of BSS with non-negativity constraints. Non-negativity is a sensible requirement as the networks that are detected by decomposition techniques can only be present or absent, but not negatively present. Thus, negative weights for networks are difficult or even impossible to interpret in this context. Note that the non-negativity requirement still permits anticorrelations between networks.

In the context of fMRI, NMF permits a spatiotemporal factorization. In this case, the columns of V are time courses measured at m voxels so that the decomposition yields r spatial and n temporal components. The order r of the decomposition must be specified beforehand. As with ICA, determining a suitable order may be difficult.

Surprisingly, while NMF has been widely used in other domains of science, its application to fMRI data has not been widespread so far. Only very few publications exist so far (Ferdowsi et al., 2011; Lohmann et al., 2007; Potluru and Calhoun, 2008; Zhang et al., 2011). Given the large potential that this method might have, we have included it in this review.

Classifiers employed for decoding

Decoding techniques employ classifiers, which, put simply, are functions that partition a set of data into two or more classes (Langford, 2006). The input data space X consists of examples (or observations) \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\vec{x} \in X , \vec{x} = [x_1 , x_2 , \ldots , x_v]$$ \end{document} . The v elements of the examples are called features. For instance, let X be the data space of weather recordings, where \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$x \in X$$ \end{document} is an example, which corresponds to one weather measurement of all weather stations. The individual elements of \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\vec {x}$$ \end{document} (i.e., the features) then are the data value of a single weather station at the given time.

The output space of the classifier are the labels Y={−1, 1}, given a two class paradigm. Crucially, there exists a distribution D: X×Y. Using this notation, we can define a classifier as a function \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$f (\vec{x}) = y$$ \end{document} , which maps an example \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\vec {x} \in X$$ \end{document} to a label \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$y \in Y$$ \end{document} . For this, the classifier first has to be trained on a subset of the data, which yields a decision boundary. Once trained, the classifier can then be used to predict the labels of an unseen subset of the data to assess the generalizability of the learned relationship between features and labels.

Classifiers differ in the way the decision boundary is derived. Linear classifiers achieve the classification decision by computing the dot product between the example \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\vec{x}$$ \end{document} with a weight vector \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\vec {w}$$ \end{document} the latter is derived by the training. In a two class paradigm, if the product \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\vec {w} \cdot \vec{x} > c$$ \end{document} , the classifier decides for the first class, while the classifier determines the second class for \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\vec {w} \cdot \vec{x} < c$$ \end{document} (where c is a constant). By usage of the dot product, each feature influences the decision only by its weight, while interactions across features have no direct influence (Pereira et al., 2009).

Examples for linear classifiers are Fisher's linear discriminant analysis, linear support vector machines (SVM), or pattern correlation classifiers. Nonlinear classifiers, such as the Gaussian naive Bayes, nonlinear SVMs, or artificial neural networks, provide a more flexible decision boundary, however, are more prone to overfitting of the training data (i.e., adjusting the decision boundary to potentially noisy characteristics). A comparison of multiple classifiers in the context of fMRI can be found in Misaki and colleagues (2010) and Pereira and colleagues (2009). For more information on SVM, see Schölkopf and Smola (2002).

Strategies for feature reduction

As the number of voxels in an fMRI data set is very large (about 50,000 for 3 T data and 500,000 or more for 7 T data), the dimensionality of the feature space (the classifier's input space) most commonly is reduced. Moreover, the number of available samples (e.g., scans) typically is very low, which increases the chance that the classifier overfits the data (O'Toole et al., 2007).

The reduction of the feature space in the case of fMRI usually increases performance, for instance, it appears straight forward that it may well be beneficial to exclude a voxel containing both negligible levels of information and high levels of noise. In addition, feature selection methods increase the interpretability of the decoding result, as they indicate which locations of the brain actually contain relevant information. As there exists a variety of different strategies for the selection of features, we want to briefly sketch out the most common approaches used in the context of fMRI.

Multivariate regression

Instead of classifying the fMRI data into different classes using a discrete label space (for instance, Y={ −1, 1}), it is also possible to regress the data to a continuous response variable of an interval, for example, Y=[0, 1]. As example for Y, consider an fMRI experiment where the subject continuously reports how positive his/her thoughts are. Using linear regression, it is then possible to find a set of voxels, which best describes the subject's rating. The best estimate of the rating variable Y then is given by \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}\widehat{Y} = x_1 \beta_1 + x_2 \beta_2 + \ldots + x_v \beta_v = X \beta\end{align*} \end{document}

where x_i is the time series of the i-th voxel (given a total number of v voxels) and v β-coefficients are determined. In other words, the matrix X contains v columns where each column is an fMRI time series, and β is a vector of length v, so that the time series of all voxels combined are used to explain the continuous response variable Y.

The standard method is an ordinary least squares approach, where the sum-squared approximation of the error becomes minimal (Carroll et al., 2009). Alternatively, a regularization can be introduced by applying a constraint on the β-coefficients of the regression model. For instance, L ₁ norm regularizations can be applied (introducing a constraint on the sum of the absolute values of the coefficients). The LASSO method (Tibshirani, 1996) uses such a constraint and is shown to result in sparse solutions. Further regularization constraints are proposed by Hoerl and Kennard (1970), Zou and Hastie (2005), Carroll and colleagues (2009), and Ryali and colleagues (2010). Alternatively, multivariate regression can be performed using support vector regression (Formisano et al., 2008).

Note that multivariate regression is strongly underdetermined, as in general, there are vastly more voxels than time points in fMRI data. Therefore, the results of this procedure depend heavily on the choice of constraints. For instance, sparsity constraints may be problematic as it is not clear to what extent the brain function really is sparse. As noted earlier, whole-brain participation is found even for very simple visual paradigms (Bandettini et al., 2012).

Encoding approaches

Arguably, the weak point of the classification and regression methods is the limitation of neurophysiological interpretability of the data, while classification methods allow claiming that there exists a difference in the evoked response patterns of two conditions, it is hardly possible to gain insight into the underlying mechanisms driving the formation of these patterns. Encoding approaches are complementary and attempt to fill the gap left by decoding methods. Principally, encoding methods first explicitly construct a model, which implements the experimental stimuli as input and then predict patterns of BOLD responses. Importantly, the quality of the model can be evaluated and compared, for instance, by assessing how much of the variance in the data can be captured by the model. Encoding approaches mostly find application in prediction of early sensory (mostly visual) areas (Naselaris et al., 2011; Nishimoto et al., 2011).

At this point, we want to note that a possible conceptual weakness of encoding models lies in the construction of the model, where features are computed from the experimental stimuli. Presently, at best, these features are extracted from the experimental stimuli by neurally inspired mechanisms (Nishimoto et al., 2011). It remains unclear to what extent these mechanisms reflect the actual underlying neuronal computations—a problem shared by almost all models proposed in cognitive neuropsychology. However, the data-driven nature of encoding approaches offers significant promise for more robust generalizations regarding classification strategies endogenous to the human brain.

Random networks

Random networks are useful as a benchmark against which more informative topologies can be tested. They are constructed by randomly connecting pairs of nodes according to some prespecified probability p. The average path length L in a random network is quite small and grows only slowly with network size.

Small worlds

A more interesting topology is found in so-called small-world networks. The principal idea was first published in a seminal article by Watts and Strogatz (1998). They observed that many real-life networks showed not only short paths between any two nodes, but also a strong local cohesiveness or cliquishness expressed by a parameter called clustering coefficient. This coefficient measures the extent to which friends of some node are also friends with each other.

The authors showed that such topologies may, for instance, arise by the following process. Starting from an initial regular network in which every node is connected to a fixed number of neighbors, some of these connections are randomly reassigned so that nodes that are far apart have a chance of becoming connected. These long-range connections provide shortcuts through the network and cut down on the average path length.

A number of studies have shown that small-world properties are present in fMRI data of the human brain (e.g., Achard et al., 2006; Varoquaux et al., 2012), for a review, see for example, Bassett and Bullmore (2008). However, as Zalesky and colleagues (2012) pointed out in a recent study, some care must be taken when using linear correlation as a measure for network analysis. The correlation coefficient is transitive, that is, if A and B are each strongly correlated with C, then A is also correlated with B, regardless of whether or not there is a direct connection between A and B. As a result, linear correlation shows a local cohesiveness suggestive of small-worldness, which may not show up with other connectivity measures. Linear correlation may therefore paint a biased picture of small-worldness. Zalesky and colleagues recommend to benchmark network measures against some null networks to counterbalance this problem.

Scale-free networks

The degree of a node is the number of its connections, that is, the number of nodes to which it is linked. The degree distribution is the probability distribution of these degrees over the whole network. Of particular interest are networks whose degree distribution is scale free (Barabasi and Albert, 1999).

Let p_i (k) denote the probability that node i has degree k, and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$P (k) = \frac {1} {n} \sum \nolimits_ {i = 1} ^n p_i (k)$$ \end{document} node degrees averaged over the entire network. Then, for scale-free networks, we have \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}p (k) = ck^{- \gamma}\end{align*} \end{document}

with c a proportionality constant so that if degrees are scaled, the distribution does not change: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}P (\alpha k) = c (\alpha k) ^{- \gamma} = c \alpha^{- \gamma} P (k) \propto P (k)\end{align*} \end{document}

In many real-like networks, the parameter γ is in the range between two and three. Scale-free properties have been observed in several fMRI studies (Eguíluz et al., 2005; He, 2011; Tomasi and Volkow, 2011).

Barabasi and Albert (1999) proposed a mechanism by which such networks can be constructed. The idea is to start with a small collection of initial nodes and expand by addition of new nodes such that new nodes attach preferentially to sites that are already well connected. This is called the preferential treatment principle, or “The rich get richer” principle. Networks that are constructed by this principle have a scale-free degree distribution and obey a power law.

The most important property of such distributions is that they are fat-tailed so that there exist some very few nodes with extremely high degrees. These nodes act as hubs and play central roles in these networks. These hubs provide shortcuts through the networks so that the average path lengths become very small. Scale-free networks are therefore sometimes called “ultra-smallworld” (Cohen and Havlin, 2003).

Note that a network may have small-world characteristics without having major hubs. On the other hand, scale-free networks may lack local cohesiveness so that they are not always small worlds. See Figure 3 for an illustration of these concepts.

OPEN IN VIEWER

FIG. 3.

A random network (left) and a scale-free network (right). Node degree is color coded such that red corresponds to high degrees, and blue to low degrees. Note that the scale-free network has a few distinct hubs, whereas the random network has a more even distribution of degrees. The corresponding degree distributions are shown below.

Node selection

Graph-based algorithms require specification of network nodes. In the context of fMRI-based analysis, this usually amounts to some form of presegmentation or parcellation of the imaging data so that each parcel or ROI represents a node. Several authors have proposed and investigated various strategies, and they all agree that this remains a challenging problem (e.g., Smith et al., 2011; Zalesky et al., 2010). No universally accepted parcellation method exists so far, making this point one of the main issues of graph-based analysis techniques.

Zalesky and colleagues (2010) provide a quantitative evaluation of several methods for ROI identification. They found that the overall small-world or scale-free topology remained intact regardless of the spatial scale at which ROIs were defined. However, the quantitative markers of small-worldness changed considerably depending on the spatial scale.

Perhaps, the most frequently used approach is the automated anatomical labeling (AAL), which provides an atlas-based predefined parcellation into 90 anatomical regions (Tzourio-Mazoyer et al., 2002). While this method is easy to use and therefore quite popular, its suitability is often quite problematic. Its main drawback is that the parcellation may be too crude and may fail to adequately represent individual neuroanatomy.

The number of regions of interest may also play a significant role and was recently discussed by Craddock and colleagues (2012). They reported a tradeoff between interpretability and accuracy with regard to the number of ROIs. Fewer than 200 ROIs seemed to make identification of the reported regions easier, while 600 or 1000 ROIs were more accurate in representing functional connectivity patterns. They propose a data-driven parcellation method based on spectral clustering.

Mumford and colleagues (2010) present an interesting algorithm for node selection called “weighted voxel coactivation network analysis”. The main idea is to assign two voxels to the same ROI if their fMRI time series are highly correlated and the voxels they are correlated with are also highly correlated among themselves. This roughly corresponds to the local clustering coefficient explained above.

Centrality mapping

In the following, we will review techniques for assessing hubness or centrality of network nodes. The degree of a node is the sum of its connections or in binary graphs its number. Nodes with high degrees can be assumed to have important roles within their network. Many studies have investigated changes in the node degree in response to changes in the brain state or as a consequence of pathologies (Buckner et al., 2009; Tomasi and Volkow, 2011).

For lack of space, we will not discuss other graph-based measures, and refer the reader to other reviews (e.g., Bullmore and Sporns, 2009) and a test–retest evaluation (Braun et al., 2012). For a comprehensive assessment of a range of centrality measures, see for instance, Zuo and colleagues (2012).

Note that preprocessing strategies may heavily influence results obtained by any of these methods, see also the section Preprocessing strategies. A commonly used procedure in this context is to first parcellate the data space according to some regime as described in the section Node selection. In principle, any of the measures of association described in the section “Measures of association” can be used to define links in the network. However, the most commonly used measure is the pairwise linear correlation between time courses averaged within each parcel. Correlation values are often thresholded so that the connection strength is ignored (Zuo et al., 2012).

Eigenvector centrality mapping

An alternative to degree centrality is eigenvector centrality. It was originally introduced in the field of social network research by Bonacich (1972). A later variant of eigenvector centrality plays a major role in Google's PageRank algorithm [Langville and Meyer (2006)]. Much like node centrality, it favors nodes that have high correlations with many other nodes (see Fig. 4). However, in contrast to node centrality, it specifically favors nodes that are connected to nodes that are themselves central within the network. It is this feature that mimics the preferential treatment effect typical in scale-free networks. As discussed earlier, this effect is caused by positive feedback mechanisms that are essential ingredients of complex systems.

OPEN IN VIEWER

FIG. 4.

Illustration of centrality measures and replicator dynamics. Nodes A and C are most central in both degree and eigenvector centrality (EC). In degree centrality, nodes B, D, E, are all equal because they all have the same number of connections. However, in eigenvector centrality, B is stronger than D, E because B has connections to two strong nodes, whereas D, E, have only connections to one strong and one weak node. Replicator dynamics identifies the clique structure A, B, C after 100 iterations and sets the other two nodes to zero.

Eigenvector centrality mapping (ECM) was proposed for fMRI data analysis in (Lohmann et al., 2010), and later applied in several other fMRI studies. In its original formulation, the ECM is applied at the voxel level so that each node in the network corresponds to one voxel, and no prior parcellation is required. The mathematical background is as follows. Let A denote an n×n similarity matrix with non-negative values only. The eigenvector centrality x_i of node i is defined as the i-th entry in the eigenvector belonging to the normalized largest eigenvalue of A. Let λ be the largest eigenvalue and x the corresponding eigenvector, then \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}Ax = \lambda \ x \ {\rm or \ equivalently} , x = \frac {1} {\lambda} Ax , {\rm and} \ x_i = \mu \sum \nolimits_ {j = 1} ^n a_ {ij} x_j\end{align*} \end{document}

with proportionality factor μ=1/λ so that x_i is proportional to the sum of similarity scores of all nodes connected to it. Thus, nodes that are well connected to other well-connected nodes are favored.

Eigenvector centrality is related to PCA in that both methods are based on eigenvector decompositions of similarity matrices. However, the difference is that the first principal component is really a projection of the data along the direction of the first eigenvector, not the eigenvector itself (see description of PCA in the section “Principal component analysis”). Furthermore, PCA differs from ECM in that it is based on linear correlations as a similarity metric. However, linear correlations may be negative so that the first principal component is not uniquely defined because of possible multiplicities of eigenvalues. Other metrics such as spectral coherence, for example, are not customary with PCA.

Many algorithms for computing eigenvectors of symmetric matrices are known. In the present context, it suffices to find the eigenvector belonging to the largest eigenvalue. For this special case, the power iteration method (Golub and van Loan, 1996, 405ff) is one of the most efficient.

Replicator dynamics

Let us now discuss another approach for analyzing the network structure. Here, we want to identify networks of brain activity such that all members that belong to the same network interact with each other. In such a network, every element is close to every other element. Such networks can be detected using a concept well known in theoretical biology called “replicator dynamics” (RD). RD describes the growth of populations consisting of several species that interact with each other [Hofbauer and Sigmund (1988) and Schuster and Sigmund (1983)]. It is described by the following equation known as the replicator equation: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}x_i (k + 1) = x_i (k) \frac {(Ax (k)) _i} {x (k) ^T Ax (k)}\end{align*} \end{document}

where x_i , i=1,…n are membership values for each of the n elements, and the n×n matrix A is a similarity matrix whose elements must be non-negative, and k denotes the iteration step. Initially, at the first iteration k=0, all elements i=1,…n have the same weight \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$x_i = \frac {1} {n} .$$ \end{document} . However, as time progresses, elements that interact most closely with many other high-weighted elements gain in weight. Interaction with low-weighted elements becomes less and less profitable. Eventually, a small set of closely interacting elements will have received a large weight, while the remaining elements become negligible. The process terminates when a stationary value \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$\overline {x}$$ \end{document} is reached that maximizes x^T Ax. For more information on the mathematical background, see for example, Hofbauer and Sigmund (1988) and Schuster and Sigmund (1983).

RD was first applied to fMRI data in Lohmann and Bohn (2002), and later, also in Müller and colleagues (2007) and Neumann and colleagues (2005, 2006). Recently, an extension to group level analysis was proposed by Ng and colleagues (2012). RD shares some similarities with ECM. Both RD and ECM highlight networks in which interaction among highly weighted elements is favored, and RD was shown to yield scale-free networks (Santos et al., 2006).

See Figure 4 for an illustration of centrality measures and RD.

Dynamic causal modeling

Dynamic causal modeling (DCM) (Friston, 2011; Friston et al., 2003) is a technique designed to investigate the influence between brain areas using time series data measured by fMRI or EEG/MEG. DCM is a model-based method in which a set of user-specified directed graphs represent possible arrangements of functional connections between predefined brain areas. Based on a model fit and parsimony, the best model is chosen and taken to be the winner. DCM is typically applied to small networks consisting of a handful of nodes.

At the heart of DCM is the following system of differential equations, which expresses temporal dynamics of functional connections represented by their corresponding matrices: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}\dot {z} = \left(A + \sum \nolimits_{j = 1}^m u_j B^j \right) z + Cu\end{align*} \end{document}

with matrix A representing intrinsic connectivity among n brain regions, matrices B^j , j=1,…m representing connectivity modulated by m experimental manipulations, and matrix C representing influences of inputs on the neural activity. The brain regions whose connections we want to investigate are usually represented by peak activations obtained by a standard GLM-based analysis. The transformation from neuronal states to the observed hemodynamic response is modeled via the Balloon-Windkessel model (Buxton et al., 1998). A dynamic causal model is generative in that it generates synthetic time courses, which are subsequently matched against the measured data. Several extensions of this method have been proposed (Li et al., 2011; Penny et al., 2010; Stephan et al., 2008). For a review, see Friston (2011).

The validity of DCM was recently debated in a series of articles (Breakspear, in press; Friston et al., in press; Lohmann et al., 2012; Lohmann et al., in press). The main points of critique were related to combinatorial explosion, the validity of the model selection procedure, and problems with respect to model validation.

The reasons underlying the problems of DCM may be seen in view of our introductory remarks about complex systems. Complex systems are characterized by a large number of constituents, and complexity arises out of the multitude of nonlinear interactions among those many agents. Small systems with few constituents cannot be complex, even though they may be highly complicated (Baranger, 2012; Rickles et al., 2007). As was recently shown, even very simple paradigms may activate large portions of the brain (Gonzalez-Castillo et al., 2012), so that insufficient model size of DCMs may in fact be one of the key issues to be investigated in this context.

In this context, it is important to note the difference between generative and discriminative models. Generative models such as DCMs should provide a full probabilistic description of all variables in the model, whereas discriminative models only need to provide a model for some target variables conditional on the observed variables (Ng and Jordan, 2001). In other words, generative models are much more ambitious, and we would expect a generative model of a brain function to reproduce at least some of its most salient features reflecting complexity such as large-scale spatial patterns of activation. We cannot expect a model with a handful of nodes to yield a realistic picture of cognitive and, hence, complex functions.

Bayesian networks

Bayesian networks have been proposed as a tool to estimate interdependencies between brain regions from fMRI time series or coactivation patterns. This idea rests upon the concept of conditional probabilities. Consider, for example, the activation of two functional brain regions A and B. If activation of region B depends on the activation of region A, then the conditional probability of activation of B given A will be greater than the conditional probability of activation of A given B: P (B|A)>P (A|B). Conditional probabilities can thus be utilized to estimate a pairwise directionality between brain regions (Patel et al., 2006).

Bayesian networks extend this notion, representing a set of random variables and their full probabilistic dependence structure. More formally, a Bayesian network is a directed acyclic graph G that comprises

• a set of nodes or vertices V in a one-to-one correspondence with a set of random variables \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document}$$X = \{X_ \nu : \nu \in V \} $$ \end{document}

Each variable X_i is assigned a conditional probability distribution P (x_i| Pa(X_i )), where Pa(X_i ) denotes the set of parents of X_i . A variable X_j is said to be a parent of X_i in G if there is a direct link pointing from X_j to X_i . If X_i has no parents, then this probability distribution reduces to the unconditional probability distribution P(x_i ). A Bayesian network then represents the joint probability distribution \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*}P (x_1 , \ldots , x_k) = \prod_{i = 1}^k P (x_i \mid {\bf Pa} (X_i)).\end{align*} \end{document}

The structure of a Bayesian network and hence, the interdependencies between brain regions can be learned from observational data, in our context from the covariance matrix of fMRI time series. For a discussion and comparison of the various network learning techniques see, for example, Chickering and colleagues (1995), Pearl (2000), Spirtes and colleagues (2000), and Steyvers and colleagues (2003). One caveat of structure learning is that some of the links typically remain undirected due to the statistical equivalence of several networks underlying the observed data (Pearl, 2009). Despite this limitation, Smith and colleagues (2011) demonstrated in an extensive simulation assessing the suitability of different network modeling approaches for the detection of network connections and their directionality that Bayesian network techniques scored among the highest when applied to good quality fMRI data.

This technique is typically applied to small networks with few numbers of nodes. However, since it is not a generative modeling approach, lack of complexity is not an issue, so that here we do not have the problem inherent in DCM.

Importantly, the use of Bayesian networks for inferring directed functional brain networks can be extended to data obtained from several independently performed fMRI experiments (Laird et al., 2005; Neumann et al., 2010, 2011). On this so-called meta-analysis level, coactivation patterns for several brain regions rather than the covariance structure of individual fMRI time series serve as observational data. This way, inference on possible brain network structures can be grounded on data collected across multiple sources.

Such meta-analytic techniques are among the big data research directions that we mentioned in the introduction. “More is different” is particularly true for the technique described in (Neumann et al., 2010) because Bayesian network estimation requires very large amounts of data. It becomes unreliable otherwise.