Variance and Autocorrelation of the Spontaneous Slow Brain Activity

Slow (<0.1 Hz) oscillatory activity in the human brain, as measured by functional magnetic imaging, has been used to identify neural networks and their dysfunction in specific brain diseases. Its intrinsic properties may also be useful to investigate brain functions. We investigated the two functional maps: variance and first order autocorrelation coefficient (r 1). These two maps had distinct spatial distributions and the values were significantly different among the subdivisions of the precuneus and posterior cingulate cortex that were identified in functional connectivity (FC) studies. The results reinforce the functional segregation of these subdivisions and indicate that the intrinsic properties of the slow brain activity have physiological relevance. Further, we propose a sample size (degree of freedom) correction when assessing the statistical significance of FC strength with r 1 values, which enables a better understanding of the network changes related to various brain diseases.


Introduction
Spontaneous fluctuations of blood oxygen level-dependent (BOLD) signals, as measured by functional magnetic resonance imaging (fMRI), are not simply caused by random noise but represent brain functions. Functional connectivity (FC) analysis [1] of these signals, pioneered by Biswal et al. [2], has revealed various brain networks that are related to specific functions [3][4][5] and their relationship with brain diseases [6,7]. The investigation of intrinsic properties of spontaneous BOLD fluctuations also revealed the other aspects of the brain function. Garrett et al. showed the relationship between standard deviation of BOLD signals and chronological age [8] and Baria et al. showed distinct spatial distributions of BOLD signals that reflect regional functional complexity [9]. These studies suggest that the intrinsic properties of BOLD signals provide novel information about regional differentiation of the brain.
In this study, we measured the intrinsic properties of spontaneous BOLD fluctuations using the two parameters, variance and autocorrelation coefficient, both of which were extracted by the autocorrelation function. The autocorrelation function of a time series is the correlation between its past and present states; thus, a high correlation indicates that the series state does not so change over time. The autocorrelation function has been used to extract periodicity in unitary neuronal activity [10][11][12] because conventional frequency analysis such as fast Fourier transform cannot be directly applied. Application of the autocorrelation function to continuous time series data such as BOLD signals is useful in extracting two distinct properties, i.e., variance and the first-order autocorrelation coefficient (r 1 ). The autocorrelation function is usually shown after normalization so that the correlation at various lags is between 21 and 1 (that is, the correlation coefficient). The value at lag zero before normalization divided by the sample number of the time series corresponds to the variance of the data. The value of r 1 can be calculated by the correlation at lag zero and lag 1. We show here that the spatial distributions of these values are different from each other and the distributions have functional relevance. Further, r 1 is useful for the correction of the degree of freedom during the assessment of FC strength between the data from two time series.

Participants
Twenty eight healthy subjects (13 women and 15 men, mean age: 34.5+/27.3 years) were recruited for this study. All were right-handed according to the Edinburgh Handedness Inventory [13] with a mean score of 92+/210.5 and gave informed consent prior to the study. This study protocol was approved by the Ethics Committee of Wakayama Medical University and we obtained written informed consent from all participants involved in this study.

MRI Data Acquisition
Each subject's brain structural and resting state functional images were acquired on a 3 Tesla MRI (PHILIPS, The Netherlands) using a 32-channel head coil (SENSE-Head-32CH). High-resolution three-dimensional T1-weighted anatomical images were collected with the following parameters: TR = 7 ms, TE = 3.3 ms, FOV = 220 mm, Matrix scan = 256, slice thickness = 0.9 mm, flip angle = 10u. Functional data were acquired using a gradient-echo echo-planar pulse sequence sensitive to BOLD contrast [14] with the following parameters: TR = 3000 ms, TE = 30 ms, FOV = 192 mm, Matrix scan = 64, slice thickness = 3 mm, flip angle = 80u. Three runs, each of which comprised 107 volumes (for 5 min 21 s), were administered to each subject. During acquisition, the subjects were instructed to stay awake with their eyes closed.

MRI Data Analysis
Preprocessing of functional MRI (fMRI) data was conducted using SPM8 (http://www.fil.ion.ucl.ac.uk/spm) and in-house software developed with MATLAB (Mathworks, Natick, MA, USA). The first 5 volumes of each fMRI acquisition run were discarded to allow for T1-equilibration effects leaving 102 consecutive volumes per session. Rigid body translation and rotation were used to correct head motion, and spatial normalization was achieved by 12-parameter affine transformation to the International Consortium for Brain Mapping Echo-Planar Imaging template in SPM8. Each image was resampled to 2-mm isotropic voxels and spatially smoothed using an 8-mm full width at half maximum Gaussian kernel. Similarly normalized and resampled structural images were then used to extract time series data for the cerebrospinal fluid (CSF), white matter (WM), and gray matter (GM),which were used to reduce non-physiological  noise in the time series of BOLD signals (see below). Each subject's three tissue images (CSF, WM, GM) were generated using SPM8 with a probability threshold of 90%.
Exclusion of the signals unrelated to brain function (i.e., brain tissue fluctuations due to head motion, cardiac activity, and respiration) was done using CompCor [4,15] and global signal regression [16]. Briefly, CompCor includes the following steps: identification of voxels showing the highest temporal variation (top 2%), principal component analysis (PCA) of these voxels and voxels within CSF and WM, identification of the PCA components accounting for a significant proportion of the variance in the data, and exclusion of the identified signal time course for each voxel using linear regression. Temporal (band-pass) filtering (ranging from 0.01 to 0.1 Hz) removed constant offsets and linear trends over each run. The 102 preprocessed images from each session were concatenated into single four-dimensional (time and 3 spatial data) images and, thus, the data from 3 sessions for each subject were used for the following analysis.
The autocorrelation function (not normalized by dividing the values at lag 0, which is also known as autocovariance function) for each voxel of the functional GM volumes was calculated with the custom Matlab command (Fig. 1) and the variance (v) and firstorder (lag 1) autocorrelation coefficient (r 1 ) were calculated for each voxel using the following equations: where C(k) is autocorrelation at lag k of N sample data (x 1 , x 2 , …, x N ) (N = 102, in this study). Then, the mean image for 3 sessions for each subject was calculated and standardized with the mean and standard deviation among all the voxels' data. Each subject's v and r 1 Z-score maps were used to extract voxels whose data were significantly different from mean zero, as assessed by a random effect one-sample t-test. The significance level was set at p,0.05 corrected for multiple comparisons with the false discovery rate (FDR) at 5% [17]. To assess the effect of lag, we created r 2 (lag 2) and r 3 (lag 3) maps and compared with r 1 map.
The three dimensional presentations of the functional maps were created using MRIcron [18], which was also used to estimate  Table 2. High r 1 cortical regions: Brodmann's area (BA); Zscore (Z, mean (SD)). the anatomical localization of the regions of interest. At each region of interest, we extracted the mean values of v and r 1 within a radius of 4 mm for each subject and assessed the regional difference of these values.

Estimation of Effective Sample Size by the Autocorrelation Function
The cross-correlation coefficient is commonly used to evaluate the magnitude of FC between two voxels in a functional volume [2,19]. The statistical significance of the correlation between two random time series can be assessed by the t -test: where r denotes the correlation and N is the number of time points. Functional data, however, are not random [20,21], i.e., the observed N samples are not independent from each other. Thus, the size of N has to be replaced by the effective sample size for correlation test, which is estimated from the following equation [22]: where N9 is the effective sample size, and r1 and r1' are the respective first order autocorrelation coefficients of the two time series. Fig. 2 shows the effective sample sizes calculated using various values of r1 and r1' when N = 102. We tested the effect of correcting the sample size on FC analysis. The default mode network (DMN) was assessed using the FC between the ventral part of the posterior cingulate cortex (PCC) (Montreal Neurological Institude (MNI) coordinates: 2, 258, 28) [23] and other brain regions. Pearson's correlation coefficient between the ventral PCC and other brain regions was calculated for each subject. For the seed data, we used the mean time course of the voxels less than 4 mm from the center of the ventral PCC. We created 2 correlation image datasets; one was made with the The result of one-sample t-test is shown excluding non-significant voxels (p.0.05 with FDR corrected). The distribution pattern is similar to that for r 1 (Fig. 4). doi:10.1371/journal.pone.0038131.g005 values above the correlation threshold at p,0.05 when the degree of freedom was 100 (i.e., N22). The other one was made with the values above the correlation threshold at p,0.05, but the degree of freedom was corrected using the two first order autocorrelation coefficients according to Equation (5).
For the group analysis, setting a threshold at each subject data may not be appropriate because the data would not follow a Gaussian distribution even after Fisher's r to z transform. Thus, we first converted the cross-correlation coefficient values to the tvalues using Equation (4) with N' for each voxel pair calculated by Equation (5). These t values were then converted to Z-scores with the same N' to perform one-sample t-test as done in v and r 1 . We used the t-to-z transform algorithm proposed by Hughett [24]. The result was compared with the Z-score map created with the conventional procedure that uses Fisher's r to z transform. Fig. 3 shows the result of the one-sample t-test for the variance as the t-value map excluding non-significant voxels (p.0.05, FDR corrected). As seen, high variance of the spontaneous BOLD signals was found in the area of the PCC, precuneus, lateral parietal cortex, parahippocampus, superior temporal gyrus, and cerebellum. Notably, the precuneus, PCC, and medial occipital cortex were divided into several regions by the distribution of the variance. In the frontal lobe, the ventromedial prefrontal cortex (vmPFC) showed a higher level of variance than the other frontal areas. MNI coordinates and Z-scores for the high variance regions are shown in Table 1. In contrast, the insula, posterior temporal cortex, primary sensorimotor areas, and ventral striatum had low    Fig. 4 shows the result of the one-sample t-test for r 1 values as the t-value map excluding non-significant voxels (p.0.05, FDR corrected). Cortical areas generally showed high r 1 values except the insula and primary sensorimotor area. Relatively high r 1 values were seen in the lateral parietal lobe, lateral prefrontal lobe, PCC, and vmPFC. Interestingly, these distributions correspond to the regions in which myelin developed last [25] and with the lowest myelin content [26]. Notably, the distribution of the r 1 values in the PCC and precuneus is different from the distribution of the variance. The MNI coordinates and Z-scores for the high r 1 regions are shown in Table 2. Significantly low r 1 values were observed in the amygdala, hippocampus, and cerebellum, in contrast to the distribution of the variance (see Fig. 3). Fig. 5 shows the result of the one-sample t-test for r 2 values as the t-value map excluding non-significant voxels (p.0.05, FDR corrected). As seen, the distribution pattern is nearly identical to the map for r 1 (see Fig. 4). No significant voxels were found for r 3 Z-score values.

Posterior Cingulate Cortex and Precuneus
Recent studies showed that the PCC and precuneus have distinct FC properties and that these areas could be divided into several subdivisions. Because the v and r 1 maps showed significantly higher values in these regions, it is of interest to check if these two values are different among the proposed subdivisions of the PCC and precuneus. Fig. 6 shows the mean Zscore maps across all subjects for both values excluding nonsignificant voxels (p.0.05, FDR corrected). On the basis of previous studies [23,27], we seeded 2 and 4 locations in the PCC and precuneus, respectively. Their relative locations are shown in Fig. 6 and the MNI coordinates are shown in Table 3. Fig. 7 shows the mean (+/2SEM) v and r 1 Z-score values among the subjects for the PCC subdivisions. Remarkable differences in both values were seen between the two subdivisions. The ventral part of the PCC showed significantly higher v and r 1 values than the dorsal part PCC (p,0.0001, paired t-test).
Three subdivisions in the precuneus showed significantly different variance from each other (Fig. 8). The sensorimotor region had the lowest variance, and while visual region had the highest. The variance of the cognitive region variance was not significantly different from that of the transitional zone. In contrast, these regions had similar r 1 values and the subdivisions could not be differentiated by r 1 . Fig. 9 shows the distribution of r 1 for the gray matter voxels in one subject's brain. The range from 0.3 to 0.7 indicates that the effective sample size for the assessment of the cross correlation coefficient could vary from 59 to 91 when the original sample size is 102 according to Equation (5).

Effective Sample Size
The effect of sample size correction on the selection of significant (p,0.05) FC voxels is shown in Fig. 10. As seen, the area was much reduced and 46.1% of voxels in the whole brain were excluded by the correction of the sample size. The mean (+/2SD) percent of the excluded voxels across the subjects was 43.4+/24.7%. Fig. 11 shows the results of one-sample t-test for the Z-scores without (A) and with temporal autocorrelation correction (B) as the t-maps excluding non-significant voxels (p.0.05, FDR corrected). The difference (corrected map minus conventional map that uses Fisher's r to z transform) of these results is shown in Fig. 11C. Both procedures revealed the same distribution pattern but the t values were generally higher for the corrected map than the uncorrected map. Further, we checked the effect of sample size correction on the Z-score map by comparing the map created by the same procedure without sample size correction, that is the map was created by the same r to t and t to z transform using the  original sample size (n = 102). The distribution of the uncorrected map was the same as that for the corrected map. Fig. 11D shows the difference between sample size corrected map and sample size uncorrected (otherwise the same as sample size corrected map) map. As seen, the t values in DMN regions were higher for the corrected map than the uncorrected map.

Discussion
We created two human brain maps using the distribution of the variance (v) and r 1 , both of which were calculated from the autocorrelation function for each voxel. These maps showed distinct spatial distributions and had functional relevance because the precuneal and PCC subdivisions had significantly different values.
Although spontaneous fluctuations of BOLD signals has been clearly shown to have functional relevance in a number of FC studies, signal variance (or standard deviation) itself in a certain time period has not been paid much attention until recently [8,28]. One possible reason is that BOLD signal variance could be largely affected by non-physiological noise such as heart beat, respiration, and head movements. Thus, proper preprocessing of BOLD time series data is the cardinal step for the assessment of variance as shown by Garrett et al. [8].
The physiological relevance of the higher variance seen in the superior temporal gyrus, lateral parietal lobe, vmPFC, cerebellum, and parahippocampus (Fig. 3) is not known. Variance is thought to be related to efficient neural processes [8,29]. If its distribution was found to be related specifically to the resting state, the functional relationship of these regions to the DMN would be an interesting subject of a future study. One of the outstanding results in this study is that the subdivisions in the precuneus and PCC had significantly different distributions of the variance (Figs. 7 and 8). These subdivisions were proposed on the basis of previous FC studies [23,27]. Our results indicate that these regions have distinct intrinsic activity that might be important for their specific functions using the related neural networks. Further, the unique distributions seen in the precuneus and PCC (Fig. 6) suggest that these regions could be divided into even more subdivisions, which was unexpected from the findings of previous studies [27,30].
Previously proposed parameter, amplitude of low-frequency fluctuation (ALFF) [31,32] would give the same results as the variance map when the same frequency range and proper standardization procedure was used. Detail analysis using ALFF and other frequency analysis may reveal further detail functional localizations in the precueus and PCC. We used variance instead of ALFF in this study simply because variance was lag zero autocovariance (see Methods).
The autocorrelation coefficient (r 1 ) distribution map revealed a distinct difference between the two subdivisions of the PCC (Fig. 6). Namely, the ventral part of the PCC had higher r 1 and v values than the dorsal part, which provides further support for the fractionation of the PCC in addition to its distinct cytoarchitectonics [33] and differential activation and FC during a cognitive task [23]. The results indicate that future DMN studies should take these subdivisions of the PCC into account. The distribution of r 1 Figure 10. Effect of sample size correction. The top images show the distribution of the cross correlation coefficients (p,0.05, t-test for each paired voxels' data) between the ventral PCC and the other brain voxels without sample size correction (i.e., N = 102). The bottom images show the distribution of the voxels with the cross correlation coefficients that are significantly different from zero (p,0.05). The effective sample size (N') (see text) was calculated for each pair of voxels with their autocorrelation coefficients and each pair's N' was used to assess the significance of the crosscorrelation coefficient. For this subject, ,46% of voxels were revealed not to be significant after sample size correction. doi:10.1371/journal.pone.0038131.g010 in the precuneus and PCC was quite different from that of v (Fig. 6), which suggests that these values reflect distinct properties of the neural activity in the same region.
Interestingly, the map showing the high r 1 regions (Fig. 4) has some resemblance to DMN regions [7], low frequency power distribution [9,31,34], cortical hubs [6], low myelination regions [25,26], and amyloid beta deposition in Alzheimer's disease [35,36]. Similarities between low frequency distribution and DMN [31], between cortical hubs and amyloid beta distribution [6], and between amyloid beta deposition and low myelination regions [35] have been reported. It is reasonable that the distribution of r 1 reflects low frequency power distribution rather than high frequency distribution because the phase is larger for lower frequency oscillations than for higher frequency oscillations. Although the physiological relevance of the other similarities is unknown, the clinical application of r 1 maps may become useful for the diagnosis of Alzheimer's disease.
In most biological time series data, the first-order autocorrelation coefficient is the highest because the present value must be affected by the most recent past value if any relationship between past and present value exists. Indeed, using lags other than 1 was not successful because the distribution pattern for lag 2 was similar to that for lag 1 (Fig. 5). The analysis using lag 3 did not show any significant voxels probably because most of the coefficient values at lag 3 were not significantly different from the values for random process. The results also validate the procedure for sampling size correction using only first-order autocorrelation coefficient.
Sampling size (or degree of freedom) correction for the statistical assessment of FC strength is important for the time series with short acquisition time and short TR [19]. Further, the temporal correlation of BOLD signals that affects the degree of freedom can vary between subjects even in the same region, which must be taken into account when an FC study is carried out with different subject groups. This is because the difference in FC between the two groups could be due to the different temporal correlations. Temporal correlation is related to the intrinsic activity in a region but FC represents a cross-correlation between different regions. Thus, these measures detect distinct brain functions. We showed that the effective sample size estimated by each voxel's r 1 value was useful to exclude ,43% voxels, which were determined to have a significantly high correlation with the seed voxel before sample size correction. This sample size correction may also be useful to compare FC maps with different sample sizes.
The fact that the regions involved in DMN showed relatively high r 1 values (Fig. 4) raises a doubt about its robust functional connectivity, because functional connectivity strength is exaggerated by the high local r 1 values. We found it was not the case for DMN. The connectivity pattern was the same and the strength could be higher even after the sample size correction (Fig. 11). Thus, high temporal autocorrelation process in each DMN region might be caused by strong functional connectivity in the network. The usefulness of sample size correction was shown by the difference between sample size corrected map and uncorrected map (Fig. 11D). Even though each subject's Z-scores in the DMN regions should have been reduced by the sample size correction, group analysis shows that t-values were higher for the corrected map than the uncorrected map. This indicates that the variance of FC strength among subjects was reduced by the sample size correction and suggests that sample size correction increases the possibility to detect physiologically meaningful FC.
In conclusion, we showed that the distinct distribution patterns of the two parameters, v and r 1 . Detailed human brain mapping with these parameters might be useful for the identification of the estimated 150-200 cortical areas [26], because these maps seem to Figure 11. Functional connectivity strength seeding at vPCC. The results of one-sample t-test are shown as t-value maps excluding nonsignificant voxels (p,0.05 with FDR corrected) (A: the result for the conventional analysis; B: the result after the sampling size correction). Both patterns are nearly identical but the latter t-values are higher than the former as shown by the difference map (C). The difference between the sample size corrected map and sample size uncorrected map created by the same r to t and t to z transform) is shown in D. The t-values in the DMN regions for the sample size corrected map were higher than that for the uncorrected map though the Z-score distribution was the same (not shown). doi:10.1371/journal.pone.0038131.g011 represent distinct brain properties that have not been shown by other parameters. Further, we proposed a sample size correction method with r 1 , which will be important for future FC studies on brain diseases.