A compact and realistic cerebral cortical layout derived from prewhitened resting-state fMRI time series: Cherniak's adjacency rule, size law, and metamodule grouping upheld

Introduction

Resting (i.e., task-free) fMRI is becoming an increasingly employed method by which to infer functional connectivity between brain areas. This body of research since the first paper by Biswal et al., published in 1995, has led to a multitude of presumed associations among various cortical areas and possible brain networks (Sporns, 2010). These results come from correlating time series of resting BOLD fMRI data after they are filtered down to very low frequencies (typically <0.1 Hz). These time series are typically nonstationary, i.e., they contain trends but also show autocorrelations for other reasons, commonly due to autoregressive (AR) and moving average (MA) processes. Commonly, all three components (trends, AR, and MA) are present in fMRI time series. Although these processes in resting fMRI time series may be of intrinsic value with regard to their biological origin, they will have to be removed prior to correlating the time series because, if left intact, the fundamental assumption of independence of errors (residuals) in the least square regression is violated and spurious correlations will be obtained. The removal of trends, AR, and MA dependencies can be done in various ways, including the Box-Jenkins AutoRegressive Integrated Moving Average (ARIMA) procedure (Box and Jenkins, 1970; Box et al., 2008), the linear transfer modeling approach (Liu, 1991; Pankratz, 1991), or a combination of linear or nonlinear detrending and fitting an ARMA model. The time series obtained as a result of eliminating these intrinsic processes are call “innovations,” to denote that they carry new information, unrelated to their history. Innovations are practically white noise (if the elimination process above is successful) and the overall procedure is called “prewhitening,” a term coined by Press and Tukey (1956, p. 220; see also Blackman and Tukey, 1959) to denote the fact that the spectrum of the innovations is white, i.e., prewhitened. (See the Appendix below for a brief discussion of the basic issues in time series analysis and their proper statistical treatment.) Locascio et al. (1997) drew attention to the pitfalls of correlating raw fMRI time series and outlined the correct procedures for computing valid correlations. Since 1995, we have applied prewhitening to fMRI data to detect task effects (Tagaris et al., 1995, 1997) and to perform brain network analyzes during cognitive processing (Georgopoulos et al., 2005, 2006; Lewis et al., 2006). In addition, we have used this method of correlating prewhitened neural time series of magnetoencephalographic (MEG) data to investigate task differences (Leuthold et al., 2005) and to assess the potential for such correlations to discriminate among, and classify, brain diseases in the resting state (Langheim et al., 2006; Georgopoulos et al., 2007, 2010; Engdahl et al., 2010).

In a recent paper, (Christova et al., 2011) we reported on the general results from a voxel-by-voxel correlation analysis of prewhitened resting fMRI BOLD time series. The main finding was that correlations are strongest within an area, followed by correlations between the same areas in the opposite hemispheres (homotopic areas), next between an area and other areas in the same hemisphere, and lastly between an area and areas in the opposite hemisphere (excluding the homotopic ones). These findings support the hypothesis that the strength of correlation between cortical areas reflects the density of overall anatomical connectivity, given the well-known strong local connectivity within an area, the strong interhemispheric connectivity between homotopic areas, the substantial, overall ipsilateral connectivity and the relatively sparse contralateral connectivity. Thus, we concluded that those results essentially follow the known pattern of the strength of anatomical connectivity, and that they could be inferred from this pattern under the reasonable assumption that there is an ongoing neural activity and communication among the various areas. In the present paper, we present the results of a more detailed analysis of our findings (Christova et al., 2011) obtained using a hierarchical tree modeling approach to reveal relations among specific areas. The layout of the cortical areas in the resulting tree (dendrogram) resembled closely the actual cortical layout of various areas and conformed to principles of efficient, compact layout of components proposed by Cherniak (Cherniak, 1994, 1995, 2012; Cherniak et al., 2004).

Materials and Methods

We analyzed a correlation matrix obtained from prewhitened resting fMRI BOLD series in 18 human subjects, as described in a previous paper (Christova et al., 2011). Details concerning subjects and methods are given in that paper and are summarized below.

Subjects

Eighteen healthy human subjects participated in these experiments as paid volunteers. They ranged in age from 21 to 44 years; 9 were men (32.9 ± 2.2 years, mean ± SEM; range: 25–44 years) and 9 were women (25.2 ± 1.1 years; range 21–32 years). All subjects participated in the study after providing informed consent, in adherence to the Declaration of Helsinki. The study protocol was approved by the respective Institutional Review Boards.

Task

The experimental task was simple, short, did not require a practice session and engaged the brain in a stable condition. Subjects lay supine within the scanner and fixated their eyes on a spot in front of them in the center of the screen. The absence of eye movement during this fixation period was verified by using an eye tracking system (ASL eye tracker, Applied Science Laboratories, Bedford, MA). Subjects were asked to remain still. Participants wore earplugs to reduce the scanner noise.

Image Acquisition

Blood oxygenation level dependent (BOLD) contrast functional images were acquired with a whole-body 3T MRI scanner (Magnetom Trio, Siemens, Erlangen, Germany) at the Center for Magnetic Resonance Research of the University of Minnesota using a gradient echo echo-planar imaging (EPI) (T2*) sequence with the following parameters: echo time (TE) = 23 ms; repetition time (TR) = 2 s; flip angle = 90°; in-plane resolution, 3 × 3 mm; slice thickness, 3 mm without inter slice gap. Whole-brain functional volumes (N = 203) of 38 axial slices covering the whole brain, cerebellum, and brain stem were obtained for each subject. A high-resolution anatomical T1-weighted 3D flash scan was obtained with the following parameters: TE = 4.7 ms; TR = 20 ms; flip angle = 22°; in-plane resolution = 1 × 1 mm; slice thickness = 1 mm; 176 slices in total.

Data Extraction

All analyzes were performed on the BOLD time series signal acquired per individual voxel in the whole brain of each subject. Coordinates in Talairach space for each voxel, as well as the BOLD intensity for each voxel, were extracted using Brain Voyager QX (v.1.10, Brain Innovation B.V., Maastricht, The Netherlands). Slice scan time correction was performed using sinc interpolation based on the information about TR and interleaved order of slice scanning. Three dimensional motion correction was performed to correct for small head movements, if present, by spatially aligning all volumes of a subject to the first volume using rigid body transformations. The estimated parameters of translation and rotation were inspected and did not exceed 3 mm or 2°. The 3-D volumes were then aligned with the corresponding 3-D anatomical volumes and normalized to standard Talairach space (Talairach and Tournoux, 1988). Matlab (R2008b, Mathworks, Natick, MA, USA) programs were implemented to enable BOLD time series extraction from the volume time course and anatomical mask available from Brain Voyager. For each subject, 203 functional images were acquired continuously, yielding a sequence of 203 BOLD signal values per voxel; of these, the first 3 volumes were discarded, leaving a time series of 200 BOLD values for analysis. Because the coefficient of variation is higher in the vicinity of large vessels and outside of the brain (Kim et al., 1994), we analyzed only voxels with coefficient of variation of no more than 5%.

A high-speed database server called Talairach Daemon (Lancaster et al., 1997, 2000) was used for automatic brain segmentation of individual brains in Talairach space. Talairach coordinates of each voxel were used to search the Talairach Daemon database (www.talairach.org, v.2.4.2) for the Talairach label using the “single point” search option. All voxels of the gray matter of the following 26 areas of the cerebral cortex in the left and right hemispheres were analyzed (for a total of 52 areas): correlations were obtained between the following areas of the left and right hemispheres: Precentral Gyrus, Superior Frontal Gyrus, Middle Frontal Gyrus, Inferior Frontal Gyrus, Paracentral Lobule, Medial Frontal Gyrus, Postcentral Gyrus, Superior Parietal Lobule, Inferior Parietal Lobule, Angular Gyrus, Supramarginal Gyrus, Precuneus, Superior Occipital Gyrus, Middle Occipital Gyrus, Inferior Occipital Gyrus, Cuneus, Lingual Gyrus, Superior Temporal Gyrus, Middle Temporal Gyrus, Inferior Temporal Gyrus, Transverse Temporal Gyrus, Fusiform Gyrus, Cingulate Gyrus, Anterior Cingulate, Posterior Cingulate, Parahippocampal Gyrus. The mean x,y,z coordinates of each area and its average volume (across the 18 subjects) are given in Table 1.

TABLE 1

Table 1. Mean x-, y- and z-Talairach coordinates and mean volumes (mm³) of 52 areas of the cerebral cortex of 18 subjects.

Data Preprocessing: Prewhitening the Raw Bold Time Series

Initial inspection of the BOLD time series from many voxels revealed that they were non-stationary with respect to the mean and highly autocorrelated. (The variance did not vary much along the series.) Since we were interested in calculating correlations between these time series, it is required, from first principles (Box and Jenkins, 1970; Box et al., 2008), that individual series be rendered stationary and nonautocorrelated for their correlation (i.e., not spurious). For that purpose we prewhitened each series using an ARIMA (15, 1, 1) model which yielded practically white noise innovations, i.e., stationary and nonautocorrelated residuals (see Christova et al., 2011 for details).

Hierarchical Clustering Analysis

Next, we analyzed a correlation matrix obtained from the prewhitened resting fMRI BOLD series. We computed correlation coefficients between 52 cortical areas (see above), as follows. For a pair of areas, e.g., A and B, all voxel innovations time series in A were correlated with all voxel innovations time series in B. The correlation coefficients were z-transformed (Fisher, 1958; Christova et al., 2011), averaged, and the means converted back to correlations. Because correlations could take negative values, a constant (=2) was added to these means to transform them to a proximity measure. Finally, the data were averaged across subjects to obtain a 52 × 52 area proximity matrix which was subjected to a hierarchical tree modeling analysis. For that purpose we used the Hierarchical Cluster procedure of the IBM-SPSS statistical package (version 20) using the squared Euclidean distance as measure and the average between-groups linkage Unweighted Pair-Group Method using arithmetic Averages (UPGMA) as the clustering method. The result was visualized as a dendrogram in which the various cortical areas were segregated in various groups. Finally, standard statistical methods were used where needed (Snedecor and Cochran, 1989).

Permutation Testing of Cherniak's Compactness

We further analyzed the placement of cortical areas in the tree in the context of Cherniak's compact component placement theory (Cherniak, 1994, 1995, 2012; Cherniak et al., 2004). In general terms, this theory states that in a component-interacting system, an optimal strategy to optimize compactness is to place individual components as close to each other as possible, according to the strength of their interaction. Cherniak has posited three basic postulates of this theory, as follows. (1) The adjacency rule: “if components are interconnected, then they are positioned contiguously to each other, other things being equal” (Cherniak, 1995, p. 523); (2) the size law: “the larger the proportion of a total optimal system that the evaluated subsystem is, the better its optimization” (Cherniak, 2012, p. 365); and (3) the metamodule grouping: “if a set of connected components is optimally placed, then a set of metamodules, each consisting of a subset of those components in the same positions, is also optimally placed” (Cherniak et al., 2004, p. 1084). In our application, the component placement corresponds to the location of a brain area in the derived tree. We tested all three of the postulates above by a procedure comprising the following steps.

(1) We computed M = 51 successive pairwise distances D (given a total of 52 areas) between brain areas in the derived tree:

$\begin{array}{cc}{D}_{tree}^i={({(L{R}_a-L{R}_b)}^2+{(P{A}_a-P{A}_b)}^2+{(I{S}_a-I{S}_b)}^2)}^{1/2}& (1)\end{array}$

where i is an index of successive distance (i = 1, M); LR, PA, IS are Talairach coordinates: Left-Right(negative→ positive), Posterior-Anterior (negative→positive), and Inferior-Superior (negative→positive); and a, b are successive brain areas in the tree.

(2) We calculated the average distance

$\begin{array}{cc}{\overline{D}}_{tree}=\frac{\text{1}}{M}{\displaystyle \sum _{i=1}^M{D}_{tree}^i}& (2)\end{array}$

(3) We randomly permuted N times the placement of these areas in the tree and recomputed the distances D^k_perm as above (where k denotes a specific permutation, k = 1 to N) and their average ${\overline{D}}_{perm}^k$.

(4) We compared each ${\overline{D}}_{perm}^k$ to ${\overline{D}}_{tree}$ and counted how many times ${\overline{D}}_{perm}^k<{\overline{D}}_{tree}$.

(5) We computed the grand average

$\begin{array}{cc}{\overline{\overline{D}}}_{perm}=\frac{\text{1}}{N}{\displaystyle \sum _{k=1}^N{\overline{D}}_{perm}^k}& (3)\end{array}$

(6) Finally, we computed the “compactness index” for the complete set as the ratio

$\begin{array}{cc}{c}^{set}=\frac{{\overline{\overline{D}}}_{perm}}{{\overline{D}}_{tree}}& (4)\end{array}$

as a measure of gain in compactness of area placement in the tree solution, as compared to the long-term average yielded by the permutations; the higher the c, the higher the gain in compactness. This procedure was applied to the tree derived as above from all 52 cortical areas of both hemispheres (after rectifying the LR coordinates) (M = 51), and, in addition, to trees derived separately for the left and right hemispheres. In the latter cases, M = 25, given 26 areas in a hemisphere.

The analysis above tested the presence of the adjacency principle on the whole set. The postulates of the size law, and metamodule grouping were tested by applying the same procedure to random subsets of areas in the tree, as follows. Let L be the size of a subset of the tree, i.e., the number of areas in a subset, L ≤ M. The measures ${\overline{D}}_{tree}^{subset}$, ${\overline{D}}_{perm}^{subset}$, ${\overline{\overline{D}}}_{perm}^{subset}$, c^subset were computed for each subset whose size L was varied systematically; the origin of a subset was randomly selected in each permutation from the full range of allowable tree locations of serial order 1 to M-L. Since L ≤ M, the comparison of c^subset to c^set would inform us whether the size law holds, which would predict that c^subset < c^set. On the other hand, the frequency of occurrence of ${\overline{D}}_{perm}^{subset}<{\overline{D}}_{tree}^{subset}$ would inform us whether the metamodule grouping holds, which would predict that subsets of the tree would still be compact, as compared to those with random permutations, i.e., that the frequency of occurrence of the inequality above would be close or equal to zero. In this analysis, we used the following parameters: N = 1000 and 10 < L ≤ M.

Results

The dendrogram obtained is shown in Figure 1. The following can be seen. First, the tree consists of several partially overlapping clusters. Second, a section of the tree at the lowest level (red line) indicates the close clustering of interhemispheric (homotopic) areas. And third, there is a systematic progression of areas clustered, from anterior to posterior. We quantified this relation as follows. For each area, we calculated the average Talairach coordinates (Talairach and Tournoux, 1988) of all voxels and all subjects. Then we performed a multiple linear regression analysis between the serial order of the location of an area in the tree and its Left-Right (LR; negative→positive), Posterior-Anterior (PA; negative→positive), and Inferior-Superior (IS; negative→positive) coordinates, as follows:

$\begin{array}{cc}Treeorder=a+bLR+c\text{PA}\text{+}d\text{IS}& (5)\end{array}$

FIGURE 1

Figure 1. Dendrogram derived by hierarchical tree clustering of 52 cortical areas (see text for details).

Where a-d are regression coefficients. The equation obtained was:

$\begin{array}{cc}Treeorder=11.13+0.013(LR)-0.345(PA)+0.213(IS)& (6)\end{array}$

Overall, this was highly statistically significant (F_(3,48) = 53.8, P < 0.001) and with excellent fit (R² = 0.771). As expected (since the homotopic areas clustered closely together), the LR coefficient was not statistically significant (t₅₁ = 0.432, P = 0.668). The most significant was the PA coefficient (t₅₁ = −12.4, P < 0.001), followed by the IS coefficient (t₅₁ = 4.3, P < 0.001). Figure 2 illustrates these results by plotting the predicted tree order against the observed one. These results indicate an orderly progression in the tree from anterior and inferior to posterior and superior. A fronto-parietal cluster is contained within the blue box in Figure 1 and its component areas are shown in Figure 3.

FIGURE 2

Figure 2. Scatter plot of the sequence of cortical areas in the dendrogram of Figure 1 (from top to bottom) against the predicted sequence using the average Talairach coordinates (per area) as the independent variables (see text for details). Red and black denote left and right hemispheres, respectively. Homotopic areas cluster closely as doublets of red and black circles.

FIGURE 3

Figure 3. Areas of the fronto-parietal cluster, enclosed in the blue box in Figure 1, are plotted on the T1-weighted structural images of one subject (left hemisphere). Three axial slices are shown. Cortical areas of the left hemisphere are presented with the following colors: red, Cingulate Gyrus; violet, Precentral Gyrus; cyan, Postcentral Gyrus; yellow, Paracentral Lobule; green, left Inferior Parietal Lobule; blue, left Precuneus; and dark red, left Superior Parietal Lobule.

Permutation Testing

Adjacency rule

The adjacency rule states that “If areas a and b are connected, then a and b are contiguous” (Cherniak, 1994, p. 2421); in the present case, we take “connected” to mean “interacting synchronously,” as quantified by the zero-lag correlation of the prewhitened resting fMRI time series. The adjacency rule can be tested by counting how many times the average consecutive distances between areas in the permuted tree would be smaller than the average distance obtained from the derived tree, i.e., ${\overline{D}}_{perm}^k<{\overline{D}}_{tree}$. We found the following with respect to the 52 areas and N = 1000 permutations, there was not a single case where ${\overline{D}}_{perm}^k<{\overline{D}}_{tree}$; this was true also for N = 1,000,000. The same results were obtained for the separate analyzes performed for each hemisphere. These findings demonstrate that the tree layout of cortical areas is much more compact than any random permutation would yield, up to N = 1,000,000.

Size law

The size law states that subsets of a compact layout would be less compact than the whole, i.e., that c^subset < c^set. We evaluated this by computing those variables for subsets of L = 10, 11, 12,…, 51, 52 for the 52-area bilateral tree, and L = 10, 15, 20, 25 for the 26-area left and right hemisphere tree. For this analysis, 1000 subsets of each size above were evaluated using random starts along the tree (up to M − L, see “Materials and Methods”); the area placement in a subset was permuted 1000 times. This yielded 1000 values of c^subset per subset size L. Figure 4 plots the geometric mean of c^subset (average compactness index) against subset size. (Since c^subset is a ratio, it was log-transformed before averaging). It can be seen that, as predicted by the size law, c^subset is smaller for smaller subset sizes, and increases in a practically linear fashion as the subset size increases. A similar increase was observed for individual hemispheres.

FIGURE 4

Figure 4. Compactness index is plotted against subset size (see text for details).

Metamodule grouping

The metamodule grouping conjecture states that “If a set of connected components is optimally placed, then a set of metamodules, each consisting of a subset of those components in the same positions, is also optimally placed.” (Cherniak et al., 2004, p. 1084). We tested this conjecture by counting how many times the condition ${\overline{D}}_{perm}^{subset}<{\overline{D}}_{tree}^{subset}$ occurred, i.e., how many times the average distance between consecutive areas in a permuted subset of the tree would be smaller than the observed distance in the original sequence. We did not find any such instance in the testing of any subset of any size in the case of the 52-areas bilateral tree and in the case of the 26-areas left and right hemispheres. This finding shows that subsets of the tree layout are always more compactly placed than when randomly permuted, as predicted by the metamodule grouping conjecture.

Discussion

Methodological Considerations

Our study differs fundamentally from all other studies of resting fMRI in four major respects. First, our analysis of functional connectivity was based on correlations computed between prewhitened resting fMRI time series, in contrast to other studies which have correlated time series without paying due consideration to their stationarity (i.e., constant mean over a long time; Bartlett, 1978) or the presence of AR (dependence on previous values) or MA (dependence on variation of previous values) processes, all of which would render the correlations thus obtained invalid, ranging from inaccurate at best and totally spurious at worst, depending on the strength of these intrinsic factors (see “Introduction”). Although ways and methods for eliminating within-series dependencies may be available in statistical packages aimed to analyze fMRI data, it is not commonly stated what use, if any, has been made of these methods, and with what result or effectiveness in achieving that goal. Serious concerns regarding this problem and the adequacy of various approaches implemented by existing fMRI analysis packages have been raised recently (Lenoski et al., 2008; Monti, 2011; Eklund et al., 2012).

In general, the treatment of this issue in published studies takes one of three forms: (1) In some studies, it is explicitly stated that filtered (to low frequencies) but non-prewhitened (i.e., detrended plus ARMA) time series have been used for correlations; in this case, the computed correlations are in doubt. (2) In other studies, the use of packages is mentioned (e.g., SPM, FSL, etc.) but without details as to what modeling on the time series was done and whether this was effective in eliminating dependencies; in this case, the computed correlations are in doubt too due to the lack of detailed information about this issue. (3) Finally, in other studies nothing is mentioned about this problem, so there is no way to tell what has been done, and, therefore these correlations are also in doubt. To ameliorate this uncertain situation, we propose that the Durbin-Watson statistic (Durbin and Watson, 1951) be reported together with every correlation analysis. This statistic evaluates the autocorrelation function of the residuals from the correlation analysis, and small values indicate spuriousness of the correlation (see also Draper and Smith, 1981, pp. 151–169). Alternatively, the Box-Ljung test on these residuals could be reported (Ljung and Box, 1978), which also evaluates the presence of significant autocorrelations. Both of these tests are readily given in standard statistical packages (e.g., IBM-SPSS). The advantage of this proposal is that no details of the kind of prewhitening procedure need to be given (e.g., nonlinear detrending, etc.): we only need to know whether the computed correlation is potentially spurious or not, and these statistics (with their associated statistical significance) will provide that information. We discussed this issue in detail, and illustrated results from our data, in a previous publication (Christova et al., 2011).

Second, our dataset comes from extensive, voxel-by-voxel correlations, specifically from more than a billion pairs of BOLD time series, instead of correlating averaged times courses from regions of interest (ROIs). No spatial smoothing was applied to the BOLD data, hence the analysis is uncontaminated from potentially artificial correlations imputed by the commonly applied substantial spatial smoothing. Although such smoothing is appropriate for certain analyzes (e.g., the detection of task effects in ROI-based analyzes) it is inappropriate for voxel-by-voxel correlations. To our knowledge, our correlation analysis (Christova et al., 2011) is the most extensive of its kind on record.

Third, our approach differs from others in that we analyzed the whole correlation matrix of 52 available cortical areas (26 per hemisphere), unlike other studies which have typically focused on certain “seed” areas. Although this seeding is valid, the results obtained are by necessity restricted to that particular seed, and a general cortical functional connectivity cannot be obtained by stitching together findings from separate seed studies to arrive at a global picture. In contrast, we aimed from the beginning at the whole picture.

Finally, we analyzed the data using a multivariate clustering approach (hierarchical tree clustering) which enabled us to identify some general principles of cortical interactions and layout as well as derive more specific clusters of cortical areas. As is the case for every multivariate analysis (e.g., factor analysis [FA], multidimensional scaling [MDS], etc.) hierarchical tree clustering does not yield a unique solution. For example, FA will yield different solutions based on the application of a factor rotation, and the kind of rotation; and MDS will yield different solutions depending on the choice of the level of measurement (e.g., metric vs. nonmetric). The choice of parameters in such multivariate analyzes typically comes from considerations of the data themselves and the objectives of the analysis; for example, Young et al. (1995) argued convincingly for the use of nonmetric MDS for analyzing cortical connectivity. In our analysis, there were two parameters involved. The first is the distance measure for which we used the squared Euclidean distance as the appropriate measure for continuous data. The second parameter concerns the method for combining clusters for which we used the average between-groups linkage (UPGMA, see “Materials and Methods”). This method was found by Milligan (1980) to be the top performer in recovering clusters (at 99.8% correct rate), among seven other methods he investigated in hierarchical cluster analysis. As a result, this is the standard method of choice, and the default method in the IBM-SPSS statistical package we used.

Interhemispheric Integration

We observed a pervasive close clustering of homotopic (left-right) areas (Figure 1, red line). This is in keeping with the well-known strong anatomical connectivity between homotopic areas and obviously underlies interhemispheric integration (Gazzaniga, 1970). Although this fact has usually being taken for granted, interhemispheric integration has not been a prominent feature of postulated resting state networks (van den Heuvel and Hulshoff Pol, 2010).

Clustering of Cortical Areas

It can be seen in Figure 1 that cortical areas cluster in a concatenated manner. Although there is a general progression from anterior to posterior, there is also fine structure in the clustering pattern, reflecting, most probably, cooperative relations. As an example, consider the cluster highlighted by the blue box in Figure 1 and illustrated in Figure 3. It consists of several areas (always integrated across hemispheres) which have been involved in visuomotor coordination. The most closely integrated sub-cluster includes the pre- and post-central areas of both hemispheres. Also of interest is the cluster consisting of paracentral, precuneus, and superior parietal lobules, all areas involved in visuomotor coordination and anatomically interconnected (Ferraina et al., 1997; Marconi et al., 2001; Margulies et al., 2009). This cluster exemplifies the proximity arrangement of cortical areas from different lobes (e.g., frontal and parietal).

Overall clustering of cortical areas based on correlations (i.e., without focus on a specific seed) has been reported previously (Cordes et al., 2002; Salvador et al., 2005; van den Heuvel et al., 2008, 2009). However, we cannot discuss them, since these correlations were apparently computed from non-prewhitened fMRI time series (as far as we can tell from the section “Materials and Methods” of these papers), and since no information was given regarding the nonspuriousness of the computed correlations (e.g., by providing the Durbin-Watson statistic).

Resting State Networks

Our results indicate an orderly functional organization of cortical areas adhering to known anatomical facts and the results of modeling and optimization studies (Young, 1992; Cherniak, 1994, 1995, 2012; Young et al., 1995; Cherniak et al., 2004). No special roles or prominence or peculiarities were found regarding areas putatively involved in such networks. We cannot comment on the existing literature in that topic, given that published resting fMRI correlational studies come from correlating nonstationary and/or autocorrelated time series yielding inaccurate correlations. The net effect of this inaccurate analysis is that claims based on such studies, and postulated brain networks derived from them, cannot be evaluated. The reason that resting state networks inferred from such correlations may be found to be in accord with known anatomical connectivity stems from the fact that an existing, true relationship is part of the equation but its relative contribution is undetermined due to the presence of nonstationarity and other factors inducing autocorrelation in the individual series (e.g., AR and MA processes).

To illustrate this point, consider the following example. The wind is blowing from the north with various gusts, and the branches and leaves of two adjacent small trees move accordingly in a direction from north to south and with various intensities, depending on the moment's gust. A recording of their motion provides two time series that consist of trends in the trees' motions induced by the wind. Now, a high correlation between these time series is found but this does not mean that the trees are connected: the correlation is due to the common wind. An analysis of the trends in the individual series will reveal very similar trends in their motion (they will differ because the trees are not identical in their physical properties), and an inference could be made about the characteristics of the wind that produced them but not about the connectivity of the trees.

However, in another variant of this example, consider that the trees are indeed tied loosely by a rope (in a functional sense). The wind blows and the tree motions are recorded. In this case, the correlation between the two time series will reflect both the influence of the wind and the constraint of the connectivity of the trees by the rope. A correlation between the two series could not really provide accurate information either for the wind or for the strength of connectivity between the trees, simply because the motion recorded is the result of both factors. Also, a child could come by and shake one of the trees: this will produce an additional ripple in the time series of the motion of that tree, and so on for other perturbations. To find the truth for the “true” connectivity between the two trees, the influence of the wind and other factors should be removed first, and this is what a prewhitening of the series would accomplish: “clean” the series by removing the trend (and other dependencies), and only then correlating them to find out the presence and strength of their functional connectivity. (Conceptually, this is similar to multiple regression analysis where the influence of certain factors are “partialed out” or “regressed out” to find the correct relation between the dependent and specific independent variable.) Finally, all of these considerations apply to frequency-domain approaches as well (Blackman and Tukey, 1959; Jenkins and Watts, 1968; Granger, 1969; Priestley, 1981).

Cortical Layout

The global organization of brain areas obtained by the hierarchical tree clustering analysis (Figure 1) resembled very much the actual layout of the cortex. The correspondence of the order of areas in the tree to the orderly progressing Talairach coordinates from anterior to posterior, primarily, and from inferior to superior, secondarily, is remarkable (Figure 2). This result is in accord with results and outcomes of modeling and optimization studies (Young, 1992; Cherniak, 1994, 1995, 2012; Cherniak et al., 2004). More specifically, our results can be interpreted in the context of Cherniak's component placement idea (Cherniak, 1994, 1995, 2012; Cherniak et al., 2004), as follows. Cherniak has proposed that the cortical layout reflects an underlying cost minimization principle, namely the cost of axon wiring between cortical areas. Although plausible, this cannot be tested rigorously in the whole cortex due to the lack of quantitative neuroanatomical information about the density of connections between areas (in monkeys or humans), a point that has been repeatedly made clear (see, e.g., Young et al., 1995; Caspers et al., 2011, 2012). Older (e.g., Caminiti et al., 1985; Barbas and Pandya, 1989) and recent studies (e.g., Barbas et al., 2005; Dancause et al., 2007; Petrides and Pandya, 2007) with an emphasis on quantifying anatomical connectivity in the monkey are usually focused on a specific set of areas (e.g., prefrontal, premotor, parietal, etc.) and, in addition, they differ substantially in the underlying tracing methodology to allow pooling of results. Be that as it may, although Cherniak's suggestion about axonal wiring cost optimization cannot be tested thoroughly, his idea about a compact cortical layout is very interesting and can be tested outside (i.e., irrespective of) anatomical considerations. The point is that, although the concept of “layout” implies structural/spatial considerations, a specific layout can be derived based on data other than structural/anatomical ones. For example, sections in a library are defined based on the content of groups of books; floors or large spaces in thematic buildings (e.g., Ministry of Education, headquarters of a company, hospitals, etc.) are allocated to similar operations (e.g., inpatient/outpatient operations, personnel department, fiscal department, etc.) within which elements (individuals, small groups of people, etc.) are doing similar operations, etc. This typically holds for aggregates of units; e.g., various research laboratories are placed near each other in larger “research areas,” as opposed to administrative offices, which are placed in “administrative areas,” etc. The cerebral cortex reflects the outcome of long-time evolutionary pressures, and its layout seems to reflect overall spatial placements of areas conforming to the strength of their zero-lag correlations. It is very possible that this interaction is associated with anatomical connectivity, or with other aspects of neural communication as well. What is remarkable, however, is that this cortical layout was derived from purely functional considerations, since the Talairach coordinates of cortical areas did not enter in any of the calculations that led to the derivation of the tree.

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.

Acknowledgments

This work was supported by United States Public Health Service grant NS32919 and the United States Department of Veterans Affairs.

References

Ashby, F. G. (2011). Statistical Analysis of fMRI Data. Cambridge, MA: MIT.

Barbas, H., and Pandya, D. N. (1989). Architecture and intrinsic connections of the prefrontal cortex in the rhesus monkey. J. Comp. Neurol. 286, 353–375.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Barbas, H., Hilgetag, C. C., Saha, S., Dermon, C. R., and Suski, J. L. (2005). Parallel organization of contralateral and ipsilateral prefrontal cortical projections in the rhesus monkey. BMC Neurosci. 6, 32. doi: 10.1186/1471-2202-6-32.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Bartlett, M. S. (1978). Stochastic Processes, 3rd Edn. Cambridge: Cambridge University Press.

Biswal, B., Yetkin, F. Z., Haughton, V. M., and Hyde, J. S. (1995). Functional connectivity in the motor cortex of resting human brain using echo-planar MRI. Magn. Reson. Med. 34, 537–541.

Pubmed Abstract | Pubmed Full Text

Blackman, R. B., and Tukey, J. W. (1959). The Measurement of Power Spectra from the Point of View of Communications Engineering, New York, NY: Dover.

Box, G. E. P., Hunter, W. G., and Hunter, J. S. (1978). Statistics for Experimenters. New York, NY: John Wiley.

Box, G. E. P., and Jenkins, G. M. (1970). Time Series Analysis: Forecasting and Control. San Francisco, CA: Holden-Day.

Box, G. E. P., Jenkins, G. M., and Reinsel, G. C. (2008). Time Series Analysis: Forecasting and Control, 4th Edn. Hoboken, NJ: Wiley.

Caminiti, R., Zeger, S., Johnson, P. B., Urbano, A., and Georgopoulos, A. P. (1985). Cortico-cortical efferent systems in the monkey: a quantitative spatial analysis of the tangential distribution of cells of origin. J. Comp. Neurol. 241, 405–419.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Caspers, S., Amunts, K., and Zilles, K. (2012). “Posterior parietal cortex: multimodal association cortex,” in The Human Nervous System, 3rd Edn. eds J. K. Mai and G. Paxinos (San Diego, CA: Academic), 1036–1055.

Caspers, S., Eickhoff, S. B., Rick, T., von Kapri, A., Kuhlen, T., Huang, R., Shah, N. J., and Zilles, K. (2011). Probabilistic fibre tract analysis of cytoarchitectonically defined human inferior lobule areas reveals similarities to macaques. Neuroimage 58, 362–380.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Cherniak, C. (1994). Component placement optimization in the brain. J. Neurosci. 14, 2418–2427.

Pubmed Abstract | Pubmed Full Text

Cherniak, C. (1995). Neural component placement. Trends Neurosci. 18, 522–527.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Cherniak, C. (2012). Neural wiring optimization. Progr. Brain Res. 195, 361–371.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Cherniak, C., Mokhtarzada, Z., Rodriguez-Esteban, R., and Changizi, K. (2004). Global optimization of cerebral cortex layout. Proc. Natl. Acad. Sci. U.S.A. 101, 1081–1086.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Christova, P., Lewis, S. M., Jerde, T. A., Lynch, J. K., and Georgopoulos, A. P. (2011). True associations between resting fMRI time series based on innovations. J. Neural. Eng. 8, 046025.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Cordes, D., Haughton, V., Carew, J. D., Arfanakis, K., and Maravilla, K. (2002). Hierarchical clustering to measure connectivity in fMRI resting-state data. Magn. Reson. Imaging 20, 305–317.

Pubmed Abstract | Pubmed Full Text

Dancause, N., Barbay, S., Frost, S. B., Mahnken, J. D., and Nudo, R. J. (2007). Interhemispheric connections of the ventral premotor cortex in a new world primate. J. Comp. Neurol. 505, 701–715.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Draper, N. R., and Smith, H. (1981). Applied Regression Analysis, 2nd Edn. New York, NY: Wiley.

Durbin, J., and Watson, G. S. (1951). Testing for serial correlation in least squares regression. II. Biometrica 38, 159–177.

Pubmed Abstract | Pubmed Full Text

Eklund, A., Andersson, M., Josephson, C., Johannesson, M., and Knutsson, H. (2012). Does parametric fMRI analysis with SPM yield valid results? – an empirical study of 1484 rest datasets. Neuroimage 61, 565–578.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Engdahl, B., Leuthold, A. C., Tan, H.-M., Lewis, S. M., Winskowski, A. M., Dikel, T. N., and Georgopoulos, A. P. (2010). The synchronous neural interactions test as a functional neuromarker for post-traumatic stress disorder (PTSD): a robust classification method based on the bootstrap. J. Neural Eng. 7, 066005.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Ferraina, S., Garasto, M. R., Battaglia-Mayer, A., Ferraresi, P., Johnson, P. B., Lacquaniti, F., and Caminiti, R. (1997). Visual control of hand-reaching movement: activity in parietal area 7m. Eur. J. Neurosci. 9, 1090–1095.

Pubmed Abstract | Pubmed Full Text

Fisher, R. A. (1958). Statistical Methods for Research Workers, 13th Edn. Edinburgh, Scotland: Oliver and Boyd.

Friston, K. J., Holmes, A. P., Poline, J.-B., Grasby, S. C., Wiiliams, S. C. R., Frackowiak, R. S. J., and Turner, R. (1995). Analysis of fMRI times-series revisited. Neuroimage 2, 45–53.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Gazzaniga, M. S. (1970). The Bisected Brain. New York, NY: Appleton-Century-Crofts.

Georgopoulos, A. P., Christova, P. S., Lewis, S. M., Slagle, E., and Ugurbil, K. (2005). Co-processing of mental rotation and memory scanning in the brain: a voxel-by-voxel analysis of 4 Tesla fMRI activation. Soc. Neurosci. Abstr. 640.16.

Georgopoulos, A. P., Karageorgiou, E., Leuthold, A. C., Lewis, S. M., Lynch, J. K., Aslam, Z., Carpenter, A. F., Georgopoulos, A., Hemmy, L. S., Koutlas, I. G., Langheim, F. J. P., McCarten, J. R., McPherson, S. E., Pardo, J. V., Pardo, P. J., Parry, G. J., Rottunda, S. J., Segal, B. M., Sponheim, S. R., Stanwyck, J. J., Stephane, M., and Westermeyer, J. J. (2007). Synchronous neural interactions assessed by magnetoencephalography: a functional biomarker for brain disorders. J. Neural Eng. 4, 349–355.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Georgopoulos, A. P., Tan, H.-M., Lewis, S. M., Leuthold, A. C., Winskowski, A. M., Lynch, J. K., and Engdahl, B. (2010). The synchronous neural interactions test as a functional neuromarker for post-traumatic stress disorder (PTSD): a robust classification method based on the bootstrap. J. Neural Eng. 7, 16011.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Georgopoulos, A. P., Christova, P. S., and Lewis, S. M. (2006). Functional network analysis of brain areas using fMRI in four cognitive tasks. Soc. Neurosci. Abstr. 262.16.

Granger, C. W. J. (1969). Investigating causal relations by econometric models and cross-spectral methods. Econometrica 37, 424–438.

Granger, C. W. J. (2004). “Time series analysis, cointegration, and applications.” in Nobel Prize Lecture 2003, The Nobel Prizes 2003, ed T. Frängsmyr (Stockholm, Sweden: Nobel Foundation), 360–366.

Granger, C. W. J., and Newbold, P. (1977). Forecasting Economic Time Series. New York, NY: Academic.

Jenkins, G. M., and Watts, D. G. (1968). Spectral Analysis and its Applications. Oakland, CA: Holden-Day.

Kim, S.-G., Hendrich, K., Hu, X., Merkle, H., and Ugurbil, K. (1994). Potential pitfalls of functional MRI using conventional gradient-recalled echo techniques. NMR Biomed. 7, 69–74.

Pubmed Abstract | Pubmed Full Text

Lancaster, J. L., Rainey, L. H., Summerlin, J. L., Freitas, C. S., Fox, P. T., Evans, A. C., Toga, A. W., and Mazziotta, J. C. (1997). Automated labeling of the human brain: a preliminary report on the development and evaluation of a forward-transform method. Hum. Brain Mapp. 5, 238–242.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Lancaster, J. L., Woldorff, M. G., Parsons, L. M., Liotti, M., Freitas, C. S., Rainey, L., Kochunov, P. V., Nickerson, D., Mikiten, S. A., and Fox, P. T. (2000). Automated talairach atlas labels for functional brain mapping Hum. Brain Mapp. 10, 120–131.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Langheim, F. J. P., Leuthold, A. C., and Georgopoulos, A. P. (2006). Synchronous dynamic brain networks revealed by magnetoencephalography. Proc. Natl. Acad. Sci. U.S.A. 103, 455–459.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Lenoski, B., Baxter, L. C., Karam, L. J., Maisog, J., and Debbins, J. (2008). On the performance of autocorrelation estimation algorithms for fMRI analysis. IEEE J. Sel. Top. Sig. Proc. 2, 828–838.

Leuthold, A. C., Langheim, F. J. P., Lewis, S. M., and Georgopoulos, A. P. (2005). Time series analysis of magnetoencephalographic (MEG) data during copying. Exp. Brain Res. 164, 411–422.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Lewis, S. M., Jerde, T. A., Goerke, U., Ugurbil, K., and Georgopoulos, A. P. (2006). Functional network analysis of superior parietal lobule (SPL) in two spatial cognitive tasks using fMRI at 7 Tesla. Soc. Neurosci. Abstr. 262.8.

Liu, L.-M. (1991). Use of linear transfer function analysis in econometric time series modeling. Statistica Sinica 1, 503–525.

Ljung, G. M., and Box, G. E. P. (1978). On a measure of a lack of fit in time series models. Biometrika 65, 297–303.

Locascio, J. J., Jennings, P. J., Moore, C. I., and Corkin, S. (1997). Time series analysis in the time domain and resampling methods for studies of functional magnetic resonance brain imaging. Hum. Brain Mapp. 15, 168–193.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Marconi, B., Genovesio, A., Battaglia-Mayer, A., Ferraina, S., Squatrito, S., Molinari, M., Lacquaniti, F., and Caminiti, R. (2001). Eye–hand coordination during reaching. I. anatomical relationships between parietal and frontal cortex. Cereb. Cortex 11, 513–527.

Pubmed Abstract | Pubmed Full Text

Margulies, D., Vincent, J. L., Kelly, C., Lohmann, G., Uddin, L. Q., Biswal, B. B., Villringer, A., Castellanos, F. X., Milham, M. P., and Petrides, M. (2009). Precuneus shares intrinsic functional architecture in humans and monkeys. Proc. Natl. Acad. Sci. U.S.A. 106, 20069–20074.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Milligan, G. W. (1980). An examination of the effects of six types of error perturbation on fifteen clustering algorithms. Psychometrika 45, 325–342.

Monti, M. M. (2011). Statistical analysis of fMRI time series: a critical review of the GLM approach. Front. Hum. Neurosci. 5:28. doi: 10.3389/fnhum.2011.00028

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Pankratz, A. (1991). Forecasting with Dynamic Regression Models. New York, NY: Wiley.

Petrides, M., and Pandya, D. N. (2007). Efferent association pathways from the rostral prefrontal cortex in the macaque monkey. J. Neurosci. 27, 11573–11586.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Pierce, D. A. (1979). R² measures for time series. J. Am. Statist. Assoc. 74, 901–910.

Poldrack, R. A., Mumford, J. A., and Nichols, T. E. (2011). Handbook of Functional MRI Data Analysis. Cambridge, MA: Cambridge University.

Press, H., and Tukey, J. W. (1956). “Power spectral methods of analysis and their application to problems in airplane dynamics,” in Flight Test Manual, NATO, Advisory Group for Aeronautical Research and Development (Berkeley Heights, NJ: Bell System Monograph No. 2606), 1–41.

Priestley, M. B. (1981). Spectral Analysis and Time Series. San Diego, CA: Academic.

Salvador, R., Suckling, J., Coleman, M. R., Pickard, J. D., Menon, D., and Bullmore, E. (2005). Neurophysiological architecture of functional magnetic resonance images of human brain. Cereb. Cortex 15, 1332–1342.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Snedecor, G. W., and Cochran, W. G. (1989). Statistical Methods. Ames, IA: Iowa State University.

Sporns, O. (2010). Networks of the brain. Cambridge, MA: MIT.

Tagaris, G. A., Kim, S.-G., Ugurbil, K., and Georgopoulos, A. P. (1995). “Box-Jenkins intervention analysis of functional MRI data Society of Magnetic Resonance,” in Third Scientific Meeting, and European Society for Magnetic Resonance in Medicine and Biology, 12th Annual Meeting (Nice, France, August 19–25, 1995), 814.

Tagaris, G. A., Richter, W., Kim, S.-G., and Georgopoulos, A. P. (1997). Box-Jenkins intervention analysis of functional magnetic resonance imaging data. Neurosci. Res. 27, 289–294.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Talairach, J., and Tournoux, P. (1988). Co-planar Stereotaxic Atlas of the Human Brain. New York, NY: Thieme.

van den Heuvel, M. P., and Hulshoff Pol, H. E. (2010). Exploring the brain network: a review on resting-state fMRI functional connectivity. Eur. J. Neuropsychopharmacol. 20, 519–534.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

van den Heuvel, M. P., Mandl, R. C. W., and Hulshoff Pol, H. E. (2008). Normalized cut group clustering of resting-state fMRI data. PLoS ONE 3:e2001. doi: 10.1371/journal.pone.0002001

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

van den Heuvel, M. P., Mandl, R. C. W., Kahn, R. S., and Hulshoff Pol, H. E. (2009). Functionally linked resting-state networks reflect the underlying structural connectivity architecture of the human brain. Hum. Brain Mapp. 30, 3127–3141.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Worsley, K. J., and Friston, K. J. (1995). Analysis of fMRI time-series revisited – again. Neuroimage 2, 173–181.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Young, M. P. (1992). Objective analysis of the topological organization of the primate cortical visual system. Nature 358, 152–155.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Young, M. P., Scannell, J. W., O'Neill, M. A., Hilgetag, C. C., Burns, G., and Blakemore, C. (1995). Non-metric multidimensional scaling in the analysis of neuroanatomical connection data and the organization of the primate cortical visual system. Philos. Trans. R. Soc. Lond. B Biol. Sci. 348, 281–308.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Yule, G. U. (1926). Why do we sometimes get nonsense-correlations between time-series? – a study in sampling and the nature of time-series. J. Roy. Stat. Soc. 89, 1–63.

Appendix

On Correlating Time Series

It is an undisputed fact that correlation (or regression) between any two time series is valid if, and only if, the series are stationary, i.e., if the mean of the series is constant throughout the time course and thus independent of the location of the time instant where it is calculated along the series (Priestley, 1981; Box et al., 2008; Jenkins and Watts, 1968; Bartlett, 1978). If this assumption does not hold, i.e., if the series are nonstationary, then spurious correlations are obtained. This issue had been recognized in the early twentieth century (Yule, 1926) but was brought forcefully to the foreground in the 1970s in the field of econometrics (Granger and Newbold, 1977). It ultimately led to a Nobel Prize awarded in 2003 to Clive Granger for (among other achievements) demonstrating that “the statistical methods used for stationary time series could yield wholly misleading results when applied to the analysis of nonstationary data.” (From Nobel Prize citation: http://nobelprize.org/nobel_prizes/economics/laureates/2003/press.html). Granger also discussed this issue in his Nobel Prize Lecture (Granger, 2004). The solution to the problem is to render the series stationary which can be achieved by differencing (i.e., taking successive differences in the series) or by the application of various detrending procedures. The hallmark of nonstationarity is the presence of autocorrelations in the series which can be detected by calculating the autocorrelation and partial autocorrelation functions (ACF and PACF, respectively). Suitable detrending will remove autocorrelations due to trends. However, a suitably detrended series may still show significant autocorrelations due to the presence of other processes, typically AR and MA processes (Box and Jenkins, 1970). The presence of an AR process indicates a dependence of a given value on previous values (irrespective of trend), whereas the presence of a MA process indicates a dependence of a given value on the variation (“random shock”) of previous values. Such dependencies need to be removed before correlation or regression analyzes are performed to ensure that the correlation obtained is truly due to the relation between the two time series and does not simply reflect, or is contaminated by, the time history of the series themselves. In a similar approach, as that applied for detrending, AR and MA dependencies are detected and removed. When this preprocessing is complete, i.e., when trends are removed and AR and/or MA dependencies eliminated, the resulting series are stationary and nonautocorrelated, which means that they are now ready to be correlated. Because there are no dependencies on previous values, the new “clean” series are called “innovations.” The standard way to accomplish this task consists of a three stage process (Box et al., 2008): (1) identify the sources of dependencies in the series (“model identification”), (2) calculate the coefficients for these dependencies (“model estimation”), and (3) take the residuals (innovations). As the final diagnostic check, the innovations time series should be stationary and non-autocorrelated. As mentioned above, there are three major sources of dependencies, namely trends (“Integrated” series), AR, and MA processes. The three stage process above used to identify, estimate and remove these dependencies is called ARIMA, from the initials (AR Integrated MA) (Box and Jenkins, 1970). Since the innovations series are essentially white noise, the data preprocessing above is called “prewhitening,” a term coined by John Tukey in 1956 (Press and Tukey, 1956) in the context of spectral analysis of time series. Detailed exposition, discussion and remedies for these problems can be found in any textbook on time series analysis (e.g., Priestley, 1981; Box et al., 2008), as detailed in a previous paper (Christova et al., 2011). It is unclear why this fundamental point has been so pervasively ignored in functional neuroimaging, while it has been commonplace in other fields, notably econometrics, for at least the past three decades. An obvious reason is that correlating autocorrelated time series requires special care of which the general neuroscientist is unaware. For example, a fundamental assumption of least squares linear regression is that the error terms be independent, and violation of this assumption invalidates statistical testing of the significance of regression slope, since its standard error is wrongly estimated. Although this is standard textbook knowledge (see, e.g., Snedecor and Cochran, 1989), it is not generally appreciated that regressing (or correlating) autocorrelated time series violates exactly this fundamental assumption, leading to inflated correlations (Pierce, 1979). In fact, in two early influential papers on fMRI data analysis, it was explicitly assumed that the errors in regression analysis between time series are independent, an erroneous assumption (Friston et al., 1995; Worsley and Friston, 1995). The detrimental effect of non-independent errors in regression analysis has been lucidly exposed by Box et al. (1978). It is interesting that the uncertainty as to how to correlate time series is clearly present in most recent textbooks of fMRI data analysis (Ashby, 2011; Poldrack et al., 2011).