Abstract

It is known that the brain's resting-state activity (RSA) is organized in low frequency oscillations that drive network connectivity. Recent research has also shown that elements of RSA described by high-frequency and nonoscillatory properties are non-random and functionally relevant. Motivated by this research, we investigated nonoscillatory aspects of the blood-oxygen-level-dependent (BOLD) RSA using a novel method for characterizing subtle fluctuation dynamics. The metric that we develop quantifies the relative variance of the amplitude of local-maxima and local-minima in a BOLD time course (amplitude variance asymmetry; AVA). This metric reveals new properties of RSA activity, without relying on connectivity as a descriptive tool. We applied the AVA analysis to data from 3 different participant groups (2 adults, 1 child) collected from 3 different centers. The analyses show that AVA patterns a) identify 3 types of RSA profiles in adults' sensory systems b) differ in topology and pattern of dynamics in adults and children, and c) are stable across magnetic resonance scanners. Furthermore, children with higher IQ demonstrated more adult-like AVA patterns. These findings indicate that AVA reflects important and novel dimensions of brain development and RSA.

Introduction

Understanding the organizing features of resting-state brain activity is one of the dominant themes in neuroscientific investigation of human cognition. Since Biswal et al.'s (1995) pioneering work demonstrating correlated spontaneous low-frequency blood-oxygen-level-dependent (BOLD) oscillations between the left and right motor cortices during rest, many studies have identified large, distributed, and functionally related brain regions whose low-frequency BOLD signals similarly oscillate in the absence of an explicit task. Some of the most prominent examples of these resting-state networks are the ventral and dorsal attention networks (Fox et al. 2006), a default-mode network (DMN) including midline and inferior parietal regions (Buckner et al. 2008), networks within auditory and visual cortices (Cordes et al. 2000), and a basal-ganglia network (Robinson et al. 2009). The topology of these networks exhibit “small-world” features (Achard et al. 2006). The size of the networks and the strength of the correlations that bind them have been linked to task performance (Hampson et al. 2006), disease (for a review, see Zhang and Raichle 2010), development (Fair et al. 2008), and genetic background (Glahn et al. 2010).

More recently, and conjointly with some of these developments, a somewhat divergent line of research has begun to show that resting-state activity (RSA) contains information that is independent of what can be observed by studying low-frequency BOLD oscillations or the networks that they maintain. Specifically, several studies have shown the variance of BOLD fluctuations within short time frames, which cannot be driven by low-frequency components, contain information important for understanding RSA (e.g., Fair, Schlaggar et al. 2007; Hasson et al. 2009; Garrett et al. 2010) as well as non-random patterns of BOLD spikes during rest (Morgan et al. 2008; Petridou et al. forthcoming; Tagliazucchi et al. 2010, 2012). Developing on this framework, we present a novel approach for studying dynamic components of BOLD RSA. As we show, this method affords new understanding of BOLD dynamics during rest, differentiates between modes of RSA not identified to date and is correlated with both age and IQ.

In what follows we briefly review empirical evidence highlighting the importance of BOLD signal features that extend beyond low-frequency–driven connectivity. We then present the logic of the current method, which is based on studying peak asymmetry in BOLD resting-state signals and its application to 3 resting-state datasets. After presenting the results, we discuss their implications for understanding RSA and relate them to recent approaches for understanding brain activity in terms of dynamical systems.

Information Contained in Higher Frequency Fluctuations

Recent work strongly supports the view that there is more information in RSA than low-frequency oscillation. As detailed below, sporadic bursts of strong activity that have typically been considered as noise has been shown to carry meaningful information. Furthermore, studies that have employed procedures where splices from larger sections of time series are concatenated (thus reducing the impact of low-frequency oscillations) have demonstrated that variance within these splices carries meaningful information and is associated with functional connectivity. For example, Garrett et al. (2010) showed that the variance of the BOLD signal within short 20 s fixation epochs correlates with chronological age. Fair, Schlaggar et al. (2007) found that established resting-state connectivity patterns are obtained by concatenating short BOLD time series collected between tasks blocks (30 s, 17 s splices). Similarly, Hasson et al. (2009) found that BOLD activity during short rest epochs (16 s) was synchronized across cortical regions, indicating that this covariance was not driven by low-frequency components.

These findings are supported by several studies showing that an important aspect of RSA consists of spontaneous bursts of BOLD activity whose profiles are similar to hemodynamic responses evoked by external stimuli. Patterns of spontaneous BOLD spikes during rest have been documented in the primary visual cortex (Wang et al. 2008). Spikes were defined as time points where the BOLD magnitude exceeded 2 standard deviations from the mean. Timing of these spikes was correlated with activity in a distributed network, including visual association regions, pre- and post-central gyri, and inferior/middle temporal cortices. Hunter et al. (2006) used a similar criterion for studying spontaneous BOLD spikes in the primary auditory cortex during silent intervals and found that the spikes correlated across the left transverse temporal and superior temporal gyri (STG). Tagliazucchi et al. (2010) identified spontaneous BOLD peaks in the pre-central gyrus (PCG) and reported that these peaks form part of a hemodynamic component highly similar to an evoked response, with an onset-to-offset interval of around 11 s and a peak at 4 s. Furthermore, activity peaks in left PCG coincided with activity peaks in right PCG, the supplementary motor area (SMA), and several other regions. Thus, peak-activity points appear to be part of organized activity across distributed regions.

An extension of this work (Tagliazucchi et al. 2012) suggests that peak-activity points carry a disproportionally large role in driving resting-state connectivity patterns. Specifically, the authors focussed on seed regions well known for their involvement in separate resting-state networks. For each seed region, activity peaks were located, and brain regions where peaks temporally coincided with the seed region's peaks were identified. Although this analysis included as few as 8 peaks from an entire 10-min resting-state scan (i.e., around 90% reduction in data), it still managed to identify resting-state connectivity profiles precisely matching those obtained by conventional independent component analysis methods. The finding suggests that connectivity patterns associated with networks demonstrating low-frequency oscillations could also be driven by spontaneous bursts of activity in those regions (see also Morgan et al. 2008, for similar developments). In related work, Petridou et al. (forthcoming) applied a technique for identifying onsets of BOLD events in the absence of stimulus timing information to BOLD resting-state data. They identified spontaneous event-like hemodynamic responses in 4 networks, and showed that the amplitude of these responses was similar to that found for task-dependent activity, both findings consistent with the work of Tagliazucchi et al. (2010, 2012). Importantly, removing the impact of these spontaneous activities from the BOLD time series reduced correlation patterns in all networks by a considerable proportion (as much as 50%) and reduced the power in the low-frequency range by as much as 60% for certain networks. Together, these works provide further evidence that short spontaneous bursts of BOLD are important drivers of resting-state correlations.

The aforementioned studies show that a central feature of BOLD RSA, at least for certain systems, is its nonstationary property. This premise is independently supported by research that has documented nonstationary aspects of the BOLD resting-state signal and begun explaining its possible sources. A recent functional magnetic resonance imaging (fMRI) study employed a wavelet coherence analysis to quantify the substantial dynamic modulation of amplitude and/or coherence of slow intrinsic fluctuations across the DMN (Chang and Glover 2010). Furthermore, several studies using combined electroencephalography (EEG)/fMRI methods have shown that RSA is temporally structured as a cascade between quasi-stable spatiotemporal “microstates” of brain activity lasting ∼100 ms (e.g., Britz et al. 2010). This has led to the suggestion that these microstates may index rapid transitions between thought modes or distinct states of consciousness (Britz et al. 2010; Musso et al. 2010; Van De Ville et al. 2010). Furthermore, activity associated with these microstates is scale-free, exhibiting similar transition patterns across a range of temporal scales including those sampled by BOLD (Britz et al. 2010; Van De Ville et al. 2010).

Patterns of spontaneous BOLD activity outside the low-frequency range can indicate information processing or may alternatively be indices of noise. Several lines of recent theoretical and empirical work have specifically focussed on the potential functions of noise in neural systems. Some have suggested that the brain at rest is a dynamical system near the point of criticality (Durstewitz and Deco 2009; Chialvo 2010). A dynamical system at the point of criticality can shift, with little energy, to a wide set of different spatiotemporal activity patterns. The neurobiological significance of this possibility is 3-fold. First, it may explain the brain's ability to quickly react to the vast range of exogenous inputs that may perturb rest, possibly by organizing into different network configurations (Sporns 2011). Secondly, it could explain why during rest neural systems easily shift between activity configurations (McIntosh et al. 2010). Different types of human behaviors are thought to reflect transitions between states (phases) of dynamical systems (Spivey et al. 2009). In addition, neural systems, including the rabbit olfactory system (Freeman 2001), and those supporting decision-making in humans (Deco et al. forthcoming) are thought to demonstrate nonlinear dynamics. On another approach, apparent noise does not necessarily indicate information processing. Instead, internal noise within neural systems can usefully serve to enhance sensitivity to near threshold percepts, a processes referred to as stochastic facilitation (McDonnell and Ward 2011). Adding noise to a dynamical system model can result in realistic simulations of BOLD activity. Recent simulations have shown that when a network's nodes are modeled as dynamical systems, together with other important structural constraints (Ghosh et al. 2008; Deco et al. 2009), their interactions produce large-scale synchronized activity with 2 features of resting-state networks: fluctuations with maximal power in low frequencies and negative temporal correlations between different networks. Importantly, an above-minimum level of Gaussian noise is necessary for obtaining both these features (Deco et al. 2009). Recent magnetoencephalography (MEG) studies (see McIntosh et al. 2010) have quantified the degree of non-regularity in task-evoked MEG signals using multi-scale entropy (MSE, Costa et al. 2002), and have shown that greater disorder in the time course of evoked responses is associated with both aging and more accurate and consistent task performance.

To summarize, recent work has begun to identify the importance of nonstationary elements observed in the BOLD signal. It has shown that spontaneous activity has a temporal profile similar to that of evoked responses. The co-occurrence patterns of activity peaks identify separate resting-state networks, and these activity peaks may be an important driver of observed resting-state connectivity patterns.

At the present stage, however, nonstationary and high-frequency BOLD fluctuations are not well understood, and their potential for discriminating between functionally disparate resting-state dynamics has been ignored. Crucially, the ability to characterize different activity generators by studying dynamics has been overlooked. The current approach marks an initial effort to resolve these unknown attributes of resting-state dynamics.

Current Approach

Our current approach departs broadly from prior research into the dynamics of resting-state data in 2 important ways. First, our approach does not rely on spike-like burst activity to define resting-state dynamics (spikes are removed prior to analysis). We summarize in a single parameter one dynamic of an entire time series. Secondly, the subtle dynamic features quantified by our approach allow for discrimination between 2 types of previously undocumented resting-state patterns. Conceptually, our novel metric (see Fig. 1) is the relative variance of the amplitude of local-maxima and local-minima across the resting-state BOLD time series. We refer to this metric as the amplitude variance asymmetry (AVA).

Figure 1.

Schematic of an AVA analysis of a single-voxel resting-state BOLD time series. For each voxel, deflection points corresponding to local-maxima (peaks, in red) and local-minima (pits, in blue) are identified. The variance of the identified peaks' amplitude is compared with the variance of the identified pits' amplitude. If there is a reliable difference in variance, the set with the lower variance is considered the reference static state. (A) An example resting-state time series where the variance of the pits is lower than the variance of the peaks; referred to as a “floor-mode” baseline pattern. (B) The logic of hypothesis testing for AVA on the single-subject and -group level.

Figure 1.

Schematic of an AVA analysis of a single-voxel resting-state BOLD time series. For each voxel, deflection points corresponding to local-maxima (peaks, in red) and local-minima (pits, in blue) are identified. The variance of the identified peaks' amplitude is compared with the variance of the identified pits' amplitude. If there is a reliable difference in variance, the set with the lower variance is considered the reference static state. (A) An example resting-state time series where the variance of the pits is lower than the variance of the peaks; referred to as a “floor-mode” baseline pattern. (B) The logic of hypothesis testing for AVA on the single-subject and -group level.

The distribution of local-maxima (peaks) and local-minima (pits) values in a time series speaks directly to the dynamics of the series. In a completely regular oscillation, e.g., a sinusoidal wave pattern, the variance of the local-maxima and the variance of the local-minima are identical, since all local-maxima have the same value (variance = 0) and all local minima have the same value. Similarly, in a random series, the variance of the peaks and the variance of the pits are highly similar, giving an AVA not significantly different from 1. Consider, however, a system that tends to settle around a certain metabolic state some of the time but departs from it toward different states associated with higher activity levels. The activity of this system over time (see Fig. 1) may show AVA because the variance of the set of local-maxima points would be higher than that of the local minima points. Time series where the variance of the peaks reliably exceeds that of the pits indicates what we refer to as a “floor-mode” resting-state pattern where the relatively more stable activity state is the one associated with lower activity levels. Conversely, time series could demonstrate a converse “ceiling-mode” resting-state pattern when the more stable state is the one with the higher activity level. In that case, the variance of the time series' pits would be reliably greater than that of the peaks.

Although the existence of such activity patterns has not been examined to date using fMRI, prior EEG and MEG research support the possibility that these dynamic activity patterns exist in BOLD data. In a pioneering work, Stam et al. (1999) found that a small proportion of EEG epochs contained signatures of nonlinearity. More recently, using MEG Mazaheri and Jensen (2008) showed that alpha activity has a skewed oscillation pattern where the peaks of alpha oscillations have greater variance than the pits of the oscillations. The theoretical importance of such findings is that amplitude variance asymmetries are not found in linear systems and therefore differentiate between oscillatory and nonoscillatory fluctuations (although the absence of amplitude variance asymmetries does not rule out a nonlinear system; Stam et al. 1998). Note that the discussion above pertains to whether AVA reflects linear dynamics in a system generating a response, not to potential nonlinear relations between neural activity and the BOLD signal (see Deshpande et al. 2006, for initial evaluations of nonlinear interactions during the resting state).

As we demonstrate below, AVA profiles in BOLD time series represent a novel, simple, and reproducible metric that offers unique insight into RSA. Due to its formal characteristics, the method departs from prior work on describing BOLD fMRI time series features in that it provides information independent of what can be obtained by entropy or variance calculations and complements those in several ways. Importantly, AVA differs from any method that returns the same assessment for a time series y and its inverse [−1 × y], as we explain below and detail in the Discussion section.

Methods

Three different human resting-state datasets were investigated to characterize the AVA of each BOLD time series. One dataset was acquired locally using a 4-T scanner and healthy adult volunteers (University of Trento, Italy). The other 2 datasets were taken from a public source (The ADHD-200 Sample, http://fcon_1000.projects.nitrc.org/indi/adhd200/) and consisted of 3-T data: one dataset consisted of healthy adult volunteers (New York University, USA), and the other dataset consisted of data from apparently normal children and children diagnosed with attention-deficit hyperactivity disorder (ADHD, University of Oregon, USA). In what follows we describe each of these datasets, the participant cohort from which it was collected, and detail the analyses performed on each.

Children Cohort (Oregon)

Children Participants

Children data were obtained from the Oregon Health and Science University ADHD database available as part of the ADHD-200 Sample. Recruitment, screening, and informed consent procedures for these participants have been reported previously (Fair, Posner et al. 2010; Fair, Bathula et al. 2010).

A total of 79 children (aged 7–12, mean age 9.16 ± 3.93 years, 42 males) were included in the study. Thirty-seven of them (mean age 8.77 ± 3.93 years, 26 males) were diagnosed as having ADHD based on evaluations with the Kiddie Schedule for Affective Disorders and Schizophrenia (KSADS-I) administered to a parent; parent and teacher Connors' Rating Scale—3rd Edition; and a clinical review by a child psychiatrist and neuropsychologist who had to agree on the diagnosis. Intelligence was also evaluated for each child with a 3-subtest short form (Block Design, Vocabulary, and Information) of the Wechsler Intelligence Scale for Children—Fourth Edition.

Children were excluded if they did not meet criteria for ADHD or non-ADHD groups (i.e., children deemed sub-threshold by the clinicians were excluded). Children were also excluded if a history of neurological illness, chronic medical problems, sensorimotor handicap, autistic disorder, mental retardation, or significant head trauma (with loss of consciousness) was reported, or if they had evidence of psychotic disorder or bipolar disorder on the structured parent psychiatric interview.

The 42 typically developing children (mean age 8.90 ± 1.19 years, 17 males) composed a control group. They were evaluated for conduct disorder, major depressive disorder, or history of psychotic disorder, as well as for presence of ADHD.

MRI Acquisition

Images were acquired using a 3-T MRI scanner (Siemens Magnetom, Trio Tim) with a whole-body radiofrequency (RF) transmit and a 12-channel receiver RF coil (Fair et al. 2010). One structural T1-weighted image was obtained per child (3D MPRAGE, orientation = sagittal, TR = 2300 ms, TE = 3.58 ms, TI = 900 ms; flip angle = 10°; 256 × 256 matrix, voxel size = 1 × 1 × 1 mm3, total scan time = 9 min 14 s). For each child, 3 BOLD-weighted resting-state fMRI runs of 78 volumes lasting 195 s were acquired, leaving a total of 585 s of resting-state data available for analysis (TR = 2.5 s, TE = 30 ms, flip angle = 90°, FOV = 240 mm, slice thickness = 3.8 mm, in-plane resolution = 3.8 mm2, 36 slices to cover the full brain). Participants were instructed to stay still, keep their eyes open, and fixate on a standard fixation cross in the center of the display.

Adult Cohort (NYU Test–Retest Dataset)

NYU Participants

A resting-state dataset from 25 healthy adult participants (mean age 29.44 ± 8.64 years, 10 males) was obtained from The ADHD-200 Sample public distribution; the New York University Child Study Center cohort. Participants were scanned 3 times as part of a study that examined the test–retest reliability of resting-state functional connectivity measures (Shehzad et al. 2009). All participants had no history of psychiatric or neurological illness as confirmed by psychiatric clinical assessment. Informed consent was obtained prior to participation.

MRI Acquisition

Images were acquired with a 3.0-Tesla scanner (Siemens Magnetom Allegra) using a whole-body RF transmit and birdcage receive RF coil (see Shehzad et al. 2009, for details). A T1-weighted anatomical image was acquired for each participant (3D MPRAGE, orientation = sagittal, TR = 2530 ms; TE = 3.25 ms; TI = 1100 ms; flip angle = 7°; 176 slices; FOV = 256 mm, voxel size: 1.3 × 1.0 × 1.3 mm3, total scan time = 8 min 7 s). Three resting-state BOLD fMRI scans were obtained for each participant. Each scan consisted of 197 contiguous echo planar imaging (EPI) functional volumes (TR = 2 s; TE = 25 ms; flip angle = 90°, 39 slices, matrix = 64 × 64; FOV = 192 mm; acquisition voxel size = 3 × 3 × 3 mm3). Scans 2 and 3 were conducted within the second scan session, 45 min apart, and were 5–16 months (mean 11 ± 4) after the first scanning session. During the scan, participants were instructed to rest with their eyes open while the word “Relax” was centrally projected in white, against a black background.

Local Cohort (UNITN)

UNITN Participants

Thirty-one right-handed participants (mean age = 28.8 ± 7.4, 20 males) with no history of psychiatric or neurological diseases gave informed consent to participate in the study, which was approved by the ethical committee of The University of Trento.

MRI Acquisition

Images were acquired with a 4-T MRI scanner (Bruker Medical, Ettlingen, Germany) using a birdcage-transmit, 8-channel receiver head coil (USA Instruments, Inc., OH, USA). Structural images were acquired using a 3D MPRAGE sequence optimized for grey–white matter contrast, with TE/TR = 4.18/2700 ms, flip angle = 7°, 1 × 1 × 1 mm3, GRAPPA iPAT = 2, total scan time: 5 min 36 s. Two structural volumes were acquired for each participant and averaged to allow an accurate image processing within a surface-based analysis pipeline. Single-shot EPI BOLD functional images were acquired using the point-spread-function distortion correction method (Zaitsev et al. 2004). A total of 115 EPI volumes were acquired during the functional run (TR/TE = 1500/40 ms, 25 interleaved slices parallel to AC/PC, voxel size = 4 × 4 × 4 mm3, slice skip factor = 0.2, Ernst flip angle = 67°). Cardiac and respiration data were collected during the functional scan using the MRI vendor's physiological recording equipment. These were sampled at 50 Hz with a photoplethysmograph and respiration belt and saved to a data file. The 2.5 min of rest scan was the first part of a longer scanning session, and during this run participants were asked to fixate on a centrally presented visual cross.

We also collected data from 5 phantom sessions that were of the same length as the functional runs, using the same scanning parameters. The phantom consisted of a spherical (170 cm diameter) glass container with silicon oil. We analyzed these data to determine whether thermal and systemic noise fluctuations in the scanner may lead, in and of themselves, to biases in the AVA measures described below (see the Image analysis: AVA Maps section).

Image Preprocessing

The processing of the locally obtained data and those obtained from publically available repositories was largely the same, with 2 exceptions. First, processing of publically available datasets was conducted in the volume domain, with the datasets registered to Talairach space. This was done because some of the publically available structural scans contained artifacts rendering them nonsuitable for cortical parcellation, while others contained missing cortical or subcortical voxels due to the anonimization procedure used when making the data publically available. The locally obtained data were processed through a surface-domain analysis using FreeSurfer (1999), which performs registration to a cortical surface template. Secondly, for locally obtained data we also obtained concurrent measurements of cardiac and heart rate responses during the scan and the effects of these were removed prior to the analysis. For the publically available datasets, such data were not available, and a time series obtained from regions with cerebral spinal fluid was used as a proxy.

Preprocessing of Oregon and NYU Datasets

Normalization of each participant's EPI images to common space was conducted in a single transformation based on matrices generated by 3 prior steps: 1) aligning each structural image to the first image in the functional scan (generating an Aligned2EPI transformation matrix), 2) normalizing the structural image to Talairach space, and 3) combining the transformations calculated in Steps 1 and 2 to accurately describe the transform from functional data to Talairach space. The detailed procedure was as follows. In Step 1, a participant's structural image was skull stripped using FSL's BET utility (Smith 2002) and separately aligned (FLIRT v. 5.5, Jenkinson and Smith 2001) to each of the 3 resting scans available for each participant using a 7 degrees of freedom transformation (3 rotation, 3 translation, 1 global scaling where all axis are equally scaled). In Step 2, the original structural images were normalized to Talairach space using a 12 degree-of-freedom (“full affine”) transformation to the TT_icbm452 template (http://www.loni.ucla.edu/Atlases/Atlas_Detail.jsp?atlas_id=6). Then the inverse of the matrix generated in Step 1 was calculated, and this matrix was concatenated with the matrix generated in Step 2. This concatenated transformation matrix effectively describes the transform from the native resting-state EPI images to Talairach space guided by the anatomy of each subject's T1 images in a single transformation. This transform was then applied to the original EPI images with FSL's applyxfm procedure. Each step was verified manually for each participant's dataset.

The original anatomical scans were then segmented into grey matter, white matter, and cerebrospinal fluid (CSF) probability maps using FMRIB's Automated Segmentation Tool (FAST; Zhang et al. 2001). The CSF masks were then spatially registered to each participant's resting-state fMRI by applying the transformation matrixes created by the affine FLIRT transformations of the original anatomicals to each subject's resting-state image, using nearest neighbor interpolation. The CSF maps, now aligned to the resting-state fMRI images, were binarized at a threshold of 50%; every voxel with a probability of being CSF ≥50% was assigned a value of 1, with all other voxels being assigned a value of 0.

The resting-state functional images were preprocessed using AFNI (Cox 1996). The first 5 volumes were discarded from the analysis to allow for T1 stabilization effects. Slice-timing correction (3dTshift), motion correction (3dvolreg), and despiking (3dDespike, http://afni.nimh.nih.gov/pub/dist/doc/program_help/3dDespike.html) were performed on all of the images. The despiking algorithm performs time series analysis on a voxel-based basis in 2 main steps. First, a smooth curve is fit to the time series. A time series of the residuals is derived and points that depart from the fitted line by more than 2.5 standard deviations are identified as spikes. Secondly, spikes are replaced by interpolating their value with a new one that keeps the deviations from the fit within a maximum range of 4 standard deviations. The CSF mask was used to obtain a mean time series from CSF regions. This mean time series and the 6 motion parameters were included as regressors of no interest and their effect removed from time series using linear regression (3dDeconvolve). Following preprocessing, motion correction, and the removal of artifacts observed in the CSF signal, the time series of all voxels were exported to text files and analyzed using the “R” software (R_Team 2008).

Preprocessing of UNITN Dataset

Inter-subject registration to common space was performed using the procedures implemented in the FreeSurfer software package (v. 4.5, Fischl, Sereno, Dale 1999, Fischl, Sereno, Tootell et al. 1999). First, the 2 anatomical images of each participant were co-registered and averaged to increase image singal-to-noise ratio (SNR). The left and right hemispheres of each participant's volume were inflated to a surface representation. Intersubject averaging was implemented by aligning the cortical folds of each participant into the default template provided with FreeSurfer (“fsaverage”), which is specified in a spherical coordinate system (Fischl, Sereno, Tootell et al. 1999). The resulting anatomical representations were imported into AFNI's surface mapping module (Saad et al. 2004).

The first 15 resting-state volumes were discarded from the analysis to allow for T1 stabilization effects. AFNI was used for preprocessing and physiological noise correction (PN correction). Despiking of the time series was carried out as part of the PN correction workflow. To account for the potential impact of physiological noise on the BOLD time series, the physiological data were down-sampled to the acquisition rate of volume slices (16 Hz), and phase shifted to match the timing of each slice's acquisition using AFNI's retroTS.m procedure. From these 2 cardiac and respiration recordings, we derived 13 slice-based regressors (4 for the cardiac series and its harmonics, 4 for the respiratory series and its harmonics, and 5 for respiration variation of time and its harmonics, Birn et al. 2006). The variance explained by these 13 regressors was removed from the BOLD time series using the RETROICOR procedure (Glover et al. 2000) as implemented in AFNI. To determine whether physiological noise cleaning resulted in a statistically significant reduction in variance, we conducted permutation tests as follows: for each of the 31 participants' BOLD time series, we implemented the RETROICOR procedure using physiological time series from all participants apart from the one for which the BOLD series was acquired (i.e., physiological data of Participants 2–31 were applied to BOLD data of Participant 1, etc). For each permutation we derived the increase in temporal SNR (decrease in variance) induced by the procedure. This permutation process generated a sampling distribution for what impact PN correction can induce by chance. The permutation test indicated that the PN correction procedure removed above-chance variance for all but one participant, whose data were therefore not included in further analysis. For each participant, functional acquisitions were spatially aligned to a reference acquisition collected during the run, and data were corrected for slice-timing differences. We removed several sources of variance from the time series data via linear regression. These included i) 6 motion parameters estimated during the head motion correction and ii) linear, second-order, and third-order polynomial trends. Note that physiological-related variance was removed prior to this analysis stage, and for this reason no global-mean correction was performed (see Murphy et al. 2009; Scholvinck et al. 2010). No temporal or spatial smoothing was conducted as part of preprocessing of the volume data (“volume domain”). Instead, intersubject registration and spatial smoothing were carried out in the surface domain as detailed below.

Projection of functional time series to the surface domain was implemented as follows: The anatomical images were aligned to the functional data using semiautomatic procedures implemented in AFNI (Saad et al. 2009). Alignment was checked and manually adjusted if needed. The functional time series of each voxel (in 3D volumetric space) were then projected from the 3D EPI volumes onto the 2D anatomical surface. For each participant, the resulting dataset consisted of [100 time points × 196 000 vertices] per hemisphere, and it was in this 2D surface domain that all subsequent time series analyses were conducted. Spatial smoothing using a 7-mm smoothing kernel (Chung et al. 2005) was performed on the surface using AFNI's SurfSmooth utility to increase the signal-to-noise ratio of the time series. This method of group registration and data smoothing in the surface domain results in accurate reflection of individual data at the group level (Argall et al. 2005), ensures that spatial blurring avoids inclusion of white matter time series, and avoids averaging data across sulci (Desai et al. 2005).

Image Analysis: AVA Maps

The analysis to characterize AVA in the resting-state time series was the same for the 3 datasets. Time series were analyzed using the “R” software (R_Team 2008). The pastecs library (Ibanez et al. 2009) was used to identify the location of turning points in the time series. Turning points are deflections in the slope of the time series, and occur both when the slope changes from positive to negative (a local-maxima, i.e., “peak”) and from negative to positive (a local-minima, i.e., “pit”). The signal values for peaks and the pits were extracted, and the ratio of the variances calculated as: [σ2(peaks)/σ2(pits)]. We refer to this parameter as the variance ratio (VR) henceforth. Because these values are ratios, for purposes of the current analysis a value of 3/2 is as significant a departure from 1 as is a value of 2/3. For this reason, significance tests on a group level to determine whether the ratio departed from chance were performed after a natural log transform of the VR; this transformation assigns the same absolute value to a ratio and its inverse;, e.g., log(3/2) = log(2/3), changing only the sign. The logarithmic transform, log[σ2(peaks)/σ2(pits)], returns 0 when VR equals 1, returns positive numbers when VR > 1, and returns negative numbers when VR < 1.

In order to reduce the impact of high-frequency fluctuations (Morgan et al. 2008), before calculating the turning points and their variance, the time series were temporally smoothed using the following moderate 3-point temporal weighted moving average: 

$$\hbox{Signal}\_{\rm smoothed}\,\hbox{(}t\hbox{)} = 0.\hbox{25} \times \hbox{Signal}\,\hbox{(}t - \hbox{1)} + 0.\hbox{5} \times \hbox{Signal}\,\hbox{(}t\hbox{)} + 0.\hbox{25} \times \hbox{Signal}\,\hbox{(}t\hbox{ + 1)}.$$

This temporal smoothing was aimed at reducing the impact of very small deviations from local trends. Pilot analyses indicated that analyses of the raw time series resulted in qualitatively similar but statistically weaker patterns. This is due to the fact that the AVA analysis is not robust against very subtle fluctuations; very small deflections are formally identified as peaks or pits, but their inclusion skews the variance measure. For instance, in a vector such as 0, 1, 2, 1.8, 3, 4, 5, 4, 3, 2 that includes 1 distinct peak (5), the value “1.8” would be considered a pit, and “2” as a peak, resulting in strong increase of both σ2(peaks) and σ2(pits) due to what likely are instrument fluctuations. Prior work examining spontaneous co-occurrence of peaks (Morgan et al. 2008) performed a similar 3-point running average. They found that not implementing temporal smoothing strongly reduced the overlap of the results with “ground truth” results.

Group-Level Analysis

Group-Level Analysis of Oregon and NYU Datasets

Following the AVA analysis we converted the output back into AFNI file format with 3dUndump and then to neuroimaging informatics technology initiative (NIFTI) file format with 3dAFNI to NIFTI. For adults, we examined whether AVA differed from chance on the group level; this consisted of T-tests contrasting AVA to 0. For the NYU dataset, the AVA values were averaged across the 3 sessions of each participant, and for the children dataset, the AVA values were averaged across the 2 or 3 sessions available for each participant. For children we first evaluated whether there were reliable differences between the AVA of the clinical and typical populations. Given a null result for this test, the data of all children participants were treated as a single group, and AVA values were compared to 0 as for adults. Because IQ scores were provided as additional datum for each child, IQ data were included as an additional regressor in the latter analysis using an ANCOVA model. That is, the analysis of AVA for children reflects these scores after adjustment to the variance reflected by IQ. All analyses were thresholded at a single voxel controlled for multiple comparisons. We controlled for familywise error (FWE) using AFNI's AlphaSim which implements the algorithm detailed by Forman et al. (1995). This procedure estimates the critical cluster size where all voxels must pass the single-voxel threshold for a cluster not to be likely to occur by chance. The single-voxel threshold (uncorrected) was set at P < 0.005 and simulations indicated the required minimum cluster size. The simulations take as a parameter the intrinsic smoothing of the data.

Group-Level Analysis of UNITN Dataset

For the UNITN group data, the AVA analysis was conducted in common cortical surface space. Cluster-based thresholding was used so that the single-voxel threshold was set at P < 0.005 (uncorrected) and simulations were conducted to find the minimal cluster size that exceeds chance where all surface vertices passed this threshold. The algorithm underlying the analysis has been described previously for volume-based analysis (Forman et al. 1995), and it was implemented here in the surface domain using in-house software (M. Schiel, SurfAlpha, available by request).

Because physiological data were acquired for this group of participants and its impact successfully removed, we conducted a functional connectivity analysis to determine the nature of connectivity between clusters demonstrating reliable AVA effects. Functional connectivity maps were constructed for each subject as follows: for each cluster showing reliable AVA effects, the time series corresponding to all surface vertices within the cluster were averaged and a full pair-wise correlation matrix between the mean cluster time course and the mean time course of the other clusters was calculated for the participants.

To determine whether the functional connectivity maps were indicative of a single source of fluctuation or multiple ones, we subjected the pair-wise correlation matrix to a hierarchical clustering solution that generated a dendogram reflecting hierarchical connectivity for each participant (using the hclust function in the R clue package). We then subjected these individual-participant dendograms to an analysis method for determining “consensus cluster analysis” (Hornik 2005) implemented using the clue package in the statistical software R. As a first step, all individual-level dendograms were aggregated (cl_ensemble). These were then analyzed using a procedure that computes the consensus clustering of these individual-level dendograms (cl_consensus). The resulting consensus-based clustering solution optimally represents the ensemble of individual-level solutions (a centroid/consensus dendogram).

In addition, we used the clue package to compute the similarity between the individual dendogram solutions by constructing a dissimilarity matrix capturing how different each participant's clustering solution is from all other participants;, i.e., dendogram dissimilarity between each pair of participants (cl_dissimilarity). If the internal organization of participants' networks were similar, their dissimilarity on this metric should be relatively low. To quantify this issue, we used a permutation-based procedure to determine whether the mean dissimilarity between the participants' dendograms was below that expected by chance. On the single-participant level, a permutation consisted of creating a surrogate dataset as follows: a) randomly pairing the region labels to the regions' time series, b) constructing a pair-wise correlation matrix and deriving a dendogram for each participant as described above. From these surrogate data, we then calculated a single “mean dissimilarity” value for the surrogate group (by averaging the [N2N]/2 dissimilarity values). That is, each simulation on the group level consisted of constructing matched “null” connectivity matrices for each person, applying hierarchical clustering to those, and estimating the dissimilarity of the clustering solutions on the group level. This procedure was repeated 1000 times, and the 50th lowest value (defining a type I error at an alpha level of P < 0.05) was taken to be a chance value. This value was compared with the actual value from the group data.

Correlation Between Time-Series' Variance and AVA Patterns

To determine whether differences in AVA could be explained by differences in overall variance within a voxel's time series, we performed a voxelwise correlation analysis between maps of the average (across each subject's 3 resting-state scans) standard deviation of each participant's BOLD time series, and the AVA of each voxel. This analysis was performed on the NYU dataset as it consisted of 3 scans per participant, and averaging AVA and variance values across scans allowed better estimation of the relation between the 2 measures. If there is a relation, then high variance in a voxel could be associated with high imbalances in the variance of peaks and pits, e.g., 10/2 or 2/10. For this reason in this analysis, the absolute value of the log of the VR was used. Standard deviation maps were created with AFNI's 3dTstat from the same preprocessed time series that were used to create the AVA maps. Each subject's standard deviation maps were then transformed into Talairach space. The correlation between participants' AVA value and STDEV values were computed on a single-voxel level to evaluate whether across participants, higher AVA values were associated with higher temporal STDEV in the BOLD time series. This resulted in an output specifying the correlation between the two values, and the reliability of the correlation for each voxel. This result was thresholded at a single-voxel level of P < 0.005 and controlled for multiple comparisons (familywise error corrected at P < 0.05 using cluster-based thresholding).

Results

We investigated the information in spontaneous BOLD signal transients during resting state by analyzing the amplitude variance asymmetry measure (AVA; log[σ2(peaks)/σ2(pits)]) at the group level. Note that if no information is contained in this pattern, then the mean value for AVA should not depart significantly from 0 (i.e., [σ2(peaks)/σ2(pits)] should not depart from 1). We first evaluated the AVA metric using phantom data to confirm that scanner noise does not bias the results. We then evaluated AVA maps for healthy adults, revealing a dominant “floor-mode” pattern in regions that constituted separate networks as revealed by subsequent functional connectivity analyses. We verified that the results are robust and consistent across 2 different datasets of healthy adults measured with very different acquisition parameters. We also verified that the results are not a measure of time series variance. Finally, we tested the metric on a third dataset of ADHD children and age-matched controls. We found that the AVA metric did not distinguish between ADHD and healthy children, that it strongly differentiated children from adults, with several areas showing significant positive correlations between AVA values and IQ.

UNITN Results

Given the novelty of the AVA technique, we confirmed its validity using phantom data. We analyzed time series acquired during 5 phantom sessions. First, we found that scanner noise did not induce more asymmetry than would be expected by chance: only 5% of the voxels were identified as showing reliable effects on the single-scan level, quantified on the single-voxel level for differences between σ2(peaks) and σ2(pits) using Levene's test for the homogeneity of variance across conditions (Fox 1997). This corresponds exactly to the expected theoretical value. Secondly, of the 5% of phantom voxels characterized as “active,” 50% showed upward transients from baseline (a “floor-mode” profile;, i.e., VR > 1) and 50% showed downward transients (a “ceiling-mode” profile; VR < 1). As shown below, patterns found in humans departed strongly from those induced by systemic scanner noise.

In the human data, we identified cortical regions demonstrating AVA by analyzing the time series of each voxel to identify local minima (pits) and maxima (peaks), quantifying the relative variance within the sets of peaks and pits (see Fig. 1 for schematic) and then analyzing these ratios at the group level (see Methods).

A group analysis conducted over the AVA values indicated several regions demonstrating departure from chance (see Fig. 2). These regions included, with good bilateral correspondence, lateral temporal cortices, premotor regions extending to the post-central sulcus, and occipital regions that largely excluded primary visual cortex and extended dorsally toward the medial intraparietal sulcus. In all identified clusters, the mean VR substantially exceeded the chance value of 1, corresponding to a “floor-mode” pattern (see Table 1 for location of clusters in Talairach space and corresponding Brodmann areas). In this dataset, we did not find any cluster where AVAs demonstrated a “ceiling-mode” pattern, i.e., a mean VR reliably below 1.

Table 1

Clusters showing significant AVA for the adult UNITN cohort

Cluster Voxels 1 mm3 Center of mass
 
Max Tstat Coordinates (max T)
 
AVA mean (SD, Extreme) Location (max T)
 
x y z x y z Talairach location BA 
3924 −26 −70 23 6.24 −31 −83 13 1.06 (0.07, 1.57) L. MOG (LVis) 19 
1606 −42 −36 22 4.14 −55 −48 20 1.14 (0.08, 1.49) L. STG (LTemp) 40 
1755 −45 −9 48 7.24 −50 −4 42 1.11 (0.06, 1.59) L. PCG (LPreM) 
2872 29 −64 28 5.55 19 −68 52 1.10 (0.08, 1.64) R. Precuneus 
1379 62 −25 22 7.90 57 −19 17 1.08 (0.08, 1.53) R. PCG (RPreM) 40 
1358 −44 24 4.24 −42 32 1.12 (0.07, 1.40) R. PCC (RPcc) 31 
1507 57 −27 −8 5.72 59 −26 −6 1.12 (0.07, 1.44) R. MTG (RTemp) 21 
1118 11 −43 7.39 18 −49 −3 1.07 (0.09, 1.78) R. Parahipp. G. (RFus) 19 
Cluster Voxels 1 mm3 Center of mass
 
Max Tstat Coordinates (max T)
 
AVA mean (SD, Extreme) Location (max T)
 
x y z x y z Talairach location BA 
3924 −26 −70 23 6.24 −31 −83 13 1.06 (0.07, 1.57) L. MOG (LVis) 19 
1606 −42 −36 22 4.14 −55 −48 20 1.14 (0.08, 1.49) L. STG (LTemp) 40 
1755 −45 −9 48 7.24 −50 −4 42 1.11 (0.06, 1.59) L. PCG (LPreM) 
2872 29 −64 28 5.55 19 −68 52 1.10 (0.08, 1.64) R. Precuneus 
1379 62 −25 22 7.90 57 −19 17 1.08 (0.08, 1.53) R. PCG (RPreM) 40 
1358 −44 24 4.24 −42 32 1.12 (0.07, 1.40) R. PCC (RPcc) 31 
1507 57 −27 −8 5.72 59 −26 −6 1.12 (0.07, 1.44) R. MTG (RTemp) 21 
1118 11 −43 7.39 18 −49 −3 1.07 (0.09, 1.78) R. Parahipp. G. (RFus) 19 

MOG, middle occipital gyri; STG, superior temporal gyri; PCG, precentral gyri; MTG, middle temporal gyri. Short abbreviations in parentheses match region labels in Figure 3. AVA Mean and AVA Extreme values refer to parameters extracted from group-level map, per each cluster. AVA SD refers to the standard deviation across participants of the mean AVA in each cluster.

Figure 2.

Regions showing AVA for the UNITN adult cohort. Data are projected to an inflated cortical surface. LH, left hemisphere (lateral, medial, posterior view). RH, right hemisphere. Warmer colors indicate a VR increasingly >1; in all areas VR exceeded 1.0 corresponding to a “floor-mode” fluctuation profile. Table 1 reports center of mass coordinates for clusters in Talairach space.

Figure 2.

Regions showing AVA for the UNITN adult cohort. Data are projected to an inflated cortical surface. LH, left hemisphere (lateral, medial, posterior view). RH, right hemisphere. Warmer colors indicate a VR increasingly >1; in all areas VR exceeded 1.0 corresponding to a “floor-mode” fluctuation profile. Table 1 reports center of mass coordinates for clusters in Talairach space.

We then evaluated the relative frequency of transition points for time series showing statistically significant AVA. For each participant, voxels where AVA was statistically significant were identified (P < 0.05 using Levene's test), and the number of peaks in these time series was extracted. We then created a histogram for each participant of the number of peaks in these voxels, and averaged the histograms across participants. The maximal bin contained 15–20 peaks per time series (see Supplementary Figure S1), indicating that, on average, peaks occurred every 8–10 s. To obtain a more detailed understanding of the temporal dynamics of the phenomenon we conducted an analysis of the peak-to-peak intervals in voxels that were statistically significant on the group level. This analysis is particularly important because the temporal filter that we applied reduced the relative power in frequencies >0.13 Hz, with a strong cutoff at around 0.2 Hz (see Supplementary Figure S2). To study peak-to-peak timings, for each time series we generated a histogram describing the distribution of peak-to-peak intervals in the time series. For each participant, these per-voxel-histograms were collapsed across voxels belonging to clusters where AVA was significant on the group level. This generated a histogram profile per cluster. Finally, these participant-level histograms (1 per cluster) were evaluated by plotting group-level histograms (see Supplementary Figure S3). This analysis revealed 2 important data patterns. First, in most clusters, the modal peak-to-peak timing was 6 s, and in a few it was 7.5 s. Secondly, in a considerable percentage of cases (ranging from 11% to 17% across clusters), the peak-to-peak interval was 4.5 s. These data suggest the phenomenon targeted is not one of rare spikes, but continuous fluctuations, where peaks can appear in relatively quick succession.

Given that the AVA RSA profile was identified via data-driven methods that were not based on identifying spatial co-activity patterns (contrasted with, e.g., seed-based connectivity analysis, independent component analysis (ICA), or 2DTCA; cf. Morgan et al. 2008), it was important to determine whether the regions identified formed functional networks. Figure 3A presents group-level connectivity maps for AVA regions, calculated by averaging the correlation maps of all participants. As shown, these regions were associated with medium-to-strong positive correlations (largely in the range of 0.4–0.8).

Figure 3.

Network features of regions showing AVA. (A) Group-level network connectivity patterns derived by averaging single-participant connectivity maps. Colored lines indicate mean correlation as calculated via Pearson's R (black: 0.4, yellow: 0.5, green: 0.6, blue: 0.7, red: 0.8). (B) Group-level clustering solution derived by finding a centroid solution that maximizes similarity to all single-participant clustering solution. Regions demonstrating AVA are clustered in a homotopic manner corresponding to different sensory systems. PCC, posterior cingulate; Fus, fusiform gyrus; vis, visual cortex cluster; PreM, premotor cluster; Temp, temporal cortex cluster.

Figure 3.

Network features of regions showing AVA. (A) Group-level network connectivity patterns derived by averaging single-participant connectivity maps. Colored lines indicate mean correlation as calculated via Pearson's R (black: 0.4, yellow: 0.5, green: 0.6, blue: 0.7, red: 0.8). (B) Group-level clustering solution derived by finding a centroid solution that maximizes similarity to all single-participant clustering solution. Regions demonstrating AVA are clustered in a homotopic manner corresponding to different sensory systems. PCC, posterior cingulate; Fus, fusiform gyrus; vis, visual cortex cluster; PreM, premotor cluster; Temp, temporal cortex cluster.

Most importantly, connectivity patterns among these regions were similar across participants. We derived the network structure on the single-participant level, using a hierarchical clustering solution applied to the connectivity data of each participant. For each participant, this clustering procedure generated a solution (dendogram) of regions sharing reliable levels of variance. A nonparametric permutation test showed that interparticipant similarity of the dendograms exceeded what would be expected by chance (P < 0.05; see Methods), indicating reliable resemblance between participants' individual network configurations. Having established above-chance similarity in the hierarchical solution of participants' network connectivity, from the single-participant dendograms we derived a group-level “consensus” clustering solution, which represents an average clustering solution reflecting common features of the individual dendograms (see Methods). Note that this consensus solution is informative—i.e., beyond being a “best fit” solution—given our finding that similarity between participants' dendograms exceeded chance probability.

As shown in Figure 3B, the group-level consensus solution demonstrated a clear partition into 3 clusters based on the type of sensory system: the left and right pre/postcentral regions formed 1 cluster, the left and right superior temporal regions formed a second cluster, and the third cluster consisted of bilateral visual cortices, the right fusiform, and right PCC.

NYU Adult Results

AVA patterns identified for this group of participants were qualitatively similar to those found in the UNITN dataset, with bilateral AVA patterns found in ventral precentral and postcentral regions, occipital regions extending dorsally toward inferior parietal areas, and lateral temporal regions. Neuroanatomical regions associated with such activity are shown in Table 2, and the distribution of these clusters activity is shown in Figure 4. As for the UNITN dataset, in all cases the AVA pattern was a floor-mode one with higher variance for peaks than pits of the BOLD response. The overall agreement between the NYU and UNITN results is noticeable, considering the data were acquired from 2 different groups of healthy adults, at different scanners with very different acquisition parameters. We examined the spatial overlap between the NYU and UNITN AVA results by thresholding both maps at a liberal threshold (AVA > 1, P < 0.05 uncorrected) and quantifying the overlap of the 2 sets using Dice's (1945) coefficient [(2 × set_intersection)/set_union], whose value was 0.59. The same analysis applied to sets of voxels showing AVA < 1 demonstrated a much lower spatial overlap of 0.02. These data indicate a relatively high degree of spatial reproducibility in floor-mode patterns across the adult groups, but low agreement in ceiling-mode patterns.

Table 2

Clusters showing significant AVA for adult NYU cohort

Cluster Voxels 1 mm3 Center of Mass
 
Max Tstat Coordinates (max T)
 
AVA Mean (SD, Extreme) Location (max T)
 
x y z x y z Taliarach location BA 
41 686 −66 27 9.01 −31 71 1.13 (0.04, 1.24) R. Paracentral Lobule 
8967 19 −58 −10 7.45 33 −72 −24 1.11 (0.04, 1.21) R. Uvula NA 
7232 56 −11 21 6.61 60 −17 32 1.13 (0.09, 1.30) R. Postcentral G 
5890 −55 −16 21 6.32 −61 −5 13 1.14 (0.09, 1.28) L. Precentral G. 43 
1891 −52 −31 42 5.90 −55 −33 52 1.11 (0.08, 1.17) L. Postcentral G. 40 
Cluster Voxels 1 mm3 Center of Mass
 
Max Tstat Coordinates (max T)
 
AVA Mean (SD, Extreme) Location (max T)
 
x y z x y z Taliarach location BA 
41 686 −66 27 9.01 −31 71 1.13 (0.04, 1.24) R. Paracentral Lobule 
8967 19 −58 −10 7.45 33 −72 −24 1.11 (0.04, 1.21) R. Uvula NA 
7232 56 −11 21 6.61 60 −17 32 1.13 (0.09, 1.30) R. Postcentral G 
5890 −55 −16 21 6.32 −61 −5 13 1.14 (0.09, 1.28) L. Precentral G. 43 
1891 −52 −31 42 5.90 −55 −33 52 1.11 (0.08, 1.17) L. Postcentral G. 40 

AVA Mean and AVA Extreme values refer to parameters extracted from group-level map, per each cluster. AVA SD refers to the standard deviation across participants of the mean AVA in each cluster.

Figure 4.

Regions showing AVA for the NYU adult cohort. Data are projected to inflated cortical surface. LH, left hemisphere (lateral, medial, posterior view); RH, right hemisphere. Warmer colors indicate a voxel ratio increasingly >1; in all areas VR exceeded 1.0 corresponding to a “floor-mode” fluctuation profile. Table 2 reports center of mass coordinates for clusters in Talairach space.

Figure 4.

Regions showing AVA for the NYU adult cohort. Data are projected to inflated cortical surface. LH, left hemisphere (lateral, medial, posterior view); RH, right hemisphere. Warmer colors indicate a voxel ratio increasingly >1; in all areas VR exceeded 1.0 corresponding to a “floor-mode” fluctuation profile. Table 2 reports center of mass coordinates for clusters in Talairach space.

As with the UNITN dataset we conducted an analysis of peak-to-peak timings for all voxels falling in clusters showing statistically significant AVA effects. Given the lower temporal resolution for this cohort (TR = 2.0), the temporal smoothing had a stronger effect. The modal peak-to-peak interval in all clusters was 8 s (25% of the cases), with 18% of the cases falling in the 6 s interval, and 20% of the cases falling in the 10 s interval (see Supplementary Figure S4).

Potential Relation of AVA Patterns to Variance of Time Series

Could AVA be a proxy for the amount of variance in a voxel's time series? From a formal perspective, this clearly need not be the case since a time series' variance could vary without affecting AVA. To illustrate, increasing the amplitude of a regular sine-wave-like activity increases the time series' variance, but consistently maintains a VR of 1 since in a sine-wave the variance of the local minima and maxima is identical. Thus, changes in variance do not necessarily indicate changes in AVA. Conversely, a voxel with a high AVA could be associated with lower time-series variance than a voxel with low AVA.

To examine the relationship between AVA and the variance of the time series, we conducted a voxelwise analysis of the NYU datasets. For every voxel we calculated the correlation between AVA values and the standard deviation of the voxel's time series. We found no brain region where this correlation was reliable, either at typical single-voxel threshold (P < 0.005, FWE corrected for multiple comparisons using cluster threshold) or at a less strict single-voxel level (P < 0.05, FWE corrected for multiple comparisons using cluster threshold).

Oregon Children Results

After the removal of low-quality datasets as indicated by Oregon's scanning logs (these had head motion exceeding voxel size in at least 1 direction), we obtained 38 datasets of nonclinical and 30 datasets of clinically diagnosed children. IQ scores differed reliably between the 2 groups, with an average of 107 (SD = 13.8) for the ADHD-diagnosed children and an average of 119 for the typically developing children (SD = 12.8), t(65) = 3.45, P = 0.001. We first evaluated whether there were any regions where AVA patterns differed in a statistically significant way for the typical and ADHD-diagnosed children. Because of the differences in IQ scores, the evaluation of differences in AVA patterns was done using analysis of covariance that included IQ as a covariate. This analysis did not reveal any statistically significant cluster where the 2 groups differed (P < 0.005 voxel threshold uncorrected, adjusted for multiple comparisons via cluster thresholding). Because we found no reliable differences, we combined the 2 datasets in further analyses.

As seen in Figure 5 and Table 3, AVA activation patterns in children were identified in both lateral and medial frontal regions, with distributions that differed remarkably from those found in adults. The clearest pattern evident in the children group is an extensive fronto-parietal distribution of regions showing a ceiling-mode AVA pattern, with a greater variance for the pits than peaks of the BOLD response. These were found in practically the entire frontal cortex, both laterally and medially. In addition, such patterns were found in the inferior parietal lobule extending to posterior STG and posterior cingulate gyrus extending to precuneus. Two regions showed a floor-mode AVA pattern for children, one covering sections of the lingual and parahippocampal gyrus, and the other covering dorsal parts of the pre- and postcentral gyri (these were located dorsally to regions found for adults).

Table 3

Clusters showing significant AVA for children Oregon cohort

Cluster Voxels 1 mm3 Center of mass
 
Max Tstat Coordinates (max T)
 
AVA Mean (SD, Extreme) Location (max T)
 
x y z x y z Taliarach location BA 
245 781 −2 27 13 −11.0 −45 41 −2 0.84 (0.10, 0.69) L Inferior Frontal G 10 
21 939 −48 −51 33 −9.35 −53 −41 36 0.85 (0.12, 0.73) L Supramarginal G 40 
15 756 −46 31 −7.42 11 −64 32 0.86 (0.11, 0.75) R Precuneus 
15 231 −1 −51 −35 −7.00 −14 −61 −36 0.88 (0.12, 0.60) L Cerebellar Tonsil NA 
14 775 51 −46 35 −7.51 54 −49 30 0.84 (0.13, 0.75) R Supramarginal G 40 
13 414 −61 −6.0 5.24 −18 −62 1.16 (0.26, 1.25) L Lingual G 19 
4220 59 −24 −10 −6.58 64 −40 −1 0.87 (0.15, 0.74) R Mid. Temp. G 21 
3823 −34 −34 59 4.79 −37 −41 60 1.15 (0.26, 1.22) L Postcentral G 
3745 −11 −31 66 5.11 −5 −37 71 1.14 (0.29, 1.19) L Paracentral Lobule 
Cluster Voxels 1 mm3 Center of mass
 
Max Tstat Coordinates (max T)
 
AVA Mean (SD, Extreme) Location (max T)
 
x y z x y z Taliarach location BA 
245 781 −2 27 13 −11.0 −45 41 −2 0.84 (0.10, 0.69) L Inferior Frontal G 10 
21 939 −48 −51 33 −9.35 −53 −41 36 0.85 (0.12, 0.73) L Supramarginal G 40 
15 756 −46 31 −7.42 11 −64 32 0.86 (0.11, 0.75) R Precuneus 
15 231 −1 −51 −35 −7.00 −14 −61 −36 0.88 (0.12, 0.60) L Cerebellar Tonsil NA 
14 775 51 −46 35 −7.51 54 −49 30 0.84 (0.13, 0.75) R Supramarginal G 40 
13 414 −61 −6.0 5.24 −18 −62 1.16 (0.26, 1.25) L Lingual G 19 
4220 59 −24 −10 −6.58 64 −40 −1 0.87 (0.15, 0.74) R Mid. Temp. G 21 
3823 −34 −34 59 4.79 −37 −41 60 1.15 (0.26, 1.22) L Postcentral G 
3745 −11 −31 66 5.11 −5 −37 71 1.14 (0.29, 1.19) L Paracentral Lobule 

AVA Mean and AVA Extreme values refer to parameters extracted from group-level map, per each cluster. AVA SD refers to the standard deviation across participants of the mean AVA in each cluster.

Figure 5.

Regions showing AVA for the children dataset. Data are projected to inflated cortical surface. LH, left hemisphere (lateral, medial); RH, right hemisphere. Warmer colors indicate VR increasingly >1 and cooler colors indicate a voxel ratio increasingly <1.0. In contrast to adults, children showed numerous regions where the voxel ratio was <1.0, corresponding to a “ceiling-mode” fluctuation profile. Table 3 reports center of mass coordinates for clusters in Talairach space.

Figure 5.

Regions showing AVA for the children dataset. Data are projected to inflated cortical surface. LH, left hemisphere (lateral, medial); RH, right hemisphere. Warmer colors indicate VR increasingly >1 and cooler colors indicate a voxel ratio increasingly <1.0. In contrast to adults, children showed numerous regions where the voxel ratio was <1.0, corresponding to a “ceiling-mode” fluctuation profile. Table 3 reports center of mass coordinates for clusters in Talairach space.

A single region demonstrated a reliable correlation between children's AVA values and IQ, located in the posterior section of the right inferior temporal gyrus/middle occipital gyrus (see Fig. 6). The correlation was highly reliable within each of the voxels in the region (P < 0.005), with a higher AVA associated with higher IQ scores. The correlation between IQ and AVA across the entire region was also statistically significant (Pearson's r = 0.50, P < 0.0001).

Figure 6.

Brain region showing a significant correlation between σ2(peaks)/σ2(pits) and IQ for children. Scatter plot shows correlation between the two variables in the voxel showing the maximal effect in that region. Familywise error corrected at P < 0.05 controlling for cluster extent and P < 0.005 for each voxel uncorrected.

Figure 6.

Brain region showing a significant correlation between σ2(peaks)/σ2(pits) and IQ for children. Scatter plot shows correlation between the two variables in the voxel showing the maximal effect in that region. Familywise error corrected at P < 0.05 controlling for cluster extent and P < 0.005 for each voxel uncorrected.

The peak-to-peak timing analysis showed that modal peak-to-peak interval in the AVA clusters occurred within 10 s. Around 16% of the instances fell in the 7.5 s bin, and 19% of the instances in the 12.5 s bin.

Direct Evaluation of Differences Between Adult and Children AVA Patterns

A direct contrast of AVA patterns found for adults and children revealed statistically significant differences in inferior parietal regions, middle and inferior frontal gyri, and ventral medial prefrontal cortex (vmPFC), bilaterally. In all cases, the AVA values were below unity for children, but near unity for adults. There were no reliable differences in the dorsal occipital regions. Figure 7 and Table 4 show regions where the AVA maps differed between the two groups.

Table 4

Clusters where AVA patterns differed for children and adults

Cluster Voxels 1 mm3 Center of mass
 
MaxTstat Coordinates (max T)
 
Location (max T)
 
x y z x y z Taliarach location BA 
7204 42 −38 7.58 44 −43 −1 L. Inferior Frontal G. 10 
5972 −33 −44 5.31 −34 −49 −2 R. Middle Frontal G. 10 
5845 54 42 33 7.08 56 41 38 L. Inf. Par. Lobule 40 
4929 −53 43 38 5.22 −54 39 43 R. Inf. Par. Lobule 40 
4522 −31 −10 5.07 −24 −11 L. Anterior Cingulate 24 
3508 42 −7 34 5.21 46 −11 33 L. middle frontal gyri, 
1737 −1 −41 18 5.39 −2 −39 19 R. Anterior Cingulate 32 
Cluster Voxels 1 mm3 Center of mass
 
MaxTstat Coordinates (max T)
 
Location (max T)
 
x y z x y z Taliarach location BA 
7204 42 −38 7.58 44 −43 −1 L. Inferior Frontal G. 10 
5972 −33 −44 5.31 −34 −49 −2 R. Middle Frontal G. 10 
5845 54 42 33 7.08 56 41 38 L. Inf. Par. Lobule 40 
4929 −53 43 38 5.22 −54 39 43 R. Inf. Par. Lobule 40 
4522 −31 −10 5.07 −24 −11 L. Anterior Cingulate 24 
3508 42 −7 34 5.21 46 −11 33 L. middle frontal gyri, 
1737 −1 −41 18 5.39 −2 −39 19 R. Anterior Cingulate 32 
Figure 7.

Regions where AVA values differed significantly between adults and children. Yellow indicates lower AVA values for children than for adults. Lower panel indicates VR values in each of the identified clusters for children and adults. Cluster numbers correspond to those in Table 4. Note that in most of these regions adults showed AVA values near unity, whereas children showed ceiling-mode effects.

Figure 7.

Regions where AVA values differed significantly between adults and children. Yellow indicates lower AVA values for children than for adults. Lower panel indicates VR values in each of the identified clusters for children and adults. Cluster numbers correspond to those in Table 4. Note that in most of these regions adults showed AVA values near unity, whereas children showed ceiling-mode effects.

Differences in AVA ratios between groups identify differences in dynamics but cannot tell what is the factor driving the difference between the groups. The finding of ceiling-mode patterns for children but near-unity values for adults in the regions identified here could be due to the fact that, e.g., children and adults have similar σ2(peaks), but that children have a larger σ2(pits), or conversely, that both groups have similar σ2(pits) but different σ2(peaks). To illustrate, an AVA value of 1 for adults and 0.5 for children can be derived from a 100:100 ratio for adults, and a 100:200 ratio for children, but these AVA values could also derive from a 100:100 ratio for adults versus a 50:100 ratio for children. To resolve this ambiguity, we extracted σ2(peaks) and σ2(pits) values in each of the clusters identified in the adults versus children contrast. The analysis showed that qualitatively, differences were mainly driven by higher σ2(pits) values for children (see Supplementary Figure S6). There were 3 clusters where σ2(pits) values were reliably or marginally higher for children, but σ2(pits) values were not. These included 1 cluster in the right inferior parietal lobule (IPL) (Cluster 4, t(60) = 2.43, P = 0.018) and 2 additional clusters (5 and 7) showing the same trend but with statistically marginal differences (P = 0.053, 0.056). There was a cluster with a different pattern, where the difference in AVA appeared to be driven by higher σ2(peaks) for adults (Cluster 6, left DLPFC; t(60) = 2.74, P < 0.01). For this specific analysis, we only included data from typically developing children. The exclusion of clinically diagnosed children was due to generally higher values of both σ2(peaks) and σ2(pits) for clinically diagnosed children (vs. the normally developing children group); this was found in all of the significant adults versus children regions, with several regions showing statistically significant effects for both parameters (Clusters 1, 2, and 5). AVA values did not differ between the 2 children groups in these regions.

In order to quantify the similarity between patterns of AVA across the 2 age groups, we thresholded each groups' AVA maps at a liberal threshold (P < 0.05 uncorrected) and conducted a Dice overlap analysis between each age group's significant clusters. We performed this analysis separately for AVA > 1 and AVA < 1. Dice's coefficient was 0.32 for the AVA > 1 sets, and 0.14 for the AVA < 1 sets. This indicates a moderate–low overlap in distribution of floor-mode patterns between adults and children, and lower overlap in ceiling-mode patterns.

Comparing Adult and Children Datasets in Functional Regions of Interest

Given the relatively low similarity in AVA patterns found for adults and children, we wanted to investigate whether children showed “adult-like” AVA patterns in regions showing significant AVA patterns for adults and, conversely, whether adults showed child-like AVA patterns in regions showing significant AVA patterns for children. To evaluate this issue, we defined clusters identified for the NYU adult cohort (see NYU Adult Results section) and the children cohort (see Oregon Children Results section) as functional regions of interest (fROIs) and examined activity levels in these regions for the converse group age.

Children AVA Patterns in Functional ROIs Defined From Adult Data

Of the 5 regions defined from the adult data, 2 showed a statistically significant floor-mode AVA pattern for children;, i.e., a pattern similar to that found for adults. These included 1 cluster in lateral occipital cortex (Cluster 1 in Table 2, t(67) = 4.1, P < 0.001) and another in the vicinity of the lingual gyrus (Cluster 2 in Table 2, t(67) = 3.2, P = 0.002). A third cluster in the left IPL (Cluster 5) showed a reliable ceiling-mode effect for children, t(67) = −3.6, P < 0.001, although it demonstrated a floor-mode pattern for adults.

To examine whether higher AVA values in these functional ROIs correlated with more mature cognitive abilities, we correlated children's AVA values in these adult ROIs with the children's IQ measures. Reliable positive correlations were found for the occipital ROIs (Clusters 1 and 2 in adult maps) where higher AVA was associated with higher IQ; Pearson's r = 0.331, P = 0.006, Pearson's r = 0.27, P < 0.05, respectively. This suggests that more adult-like AVA patterns in these regions are associated with more adult-like cognitive abilities as measured by IQ. Because in this cohort, typically developing children had a higher IQ overall, it could be that the correlation between BOLD and IQ reflected the effect of clinical state. To investigate this potential factor, we re-analyzed the correlation between AVA values and IQ in these 5 fROIS separately for the clinical and nonclinical children groups. For the nonclinical children, Clusters 1 and 2 showed marginally positive correlations with IQ (r = 0.27, P = 0.08, 0.09), in the direction reported when considering all children conjointly. For the clinically diagnosed children, none of the clusters showed a correlation that approached significance. This suggests that the correlation between IQ and AVA seen when collapsing both children groups was not driven by IQ differences between the clinical and nonclinical children. We conducted a similar analysis to see whether children's IQ was correlated with AVA values in the 9 regions identified in their own data, but none of these correlations approached significance.

Adult AVA Patterns in Functional ROIs Defined From Children Data

The converse analysis examined the adult activation patterns in regions defined by children's AVA patterns. This analysis revealed that within one functional ROI (Cluster 1 in children's results), adults showed a reliable ceiling-mode activity pattern of the sort found for children (t(24) = −2.11, P < 0.05). In 4 other functional ROIs (Clusters 2, 8, 6, and 9), adults demonstrated a reliable floor-mode AVA pattern (all P < 0.005). Note, however, that for the latter 3 regions children showed floor-mode AVA patterns as well, so this finding is not at all surprising; that is, adults demonstrated floor-mode AVA patterns in all 3 regions for which children demonstrate floor-mode patterns. However, in the left IPL (Cluster 2), children showed a reliable ceiling-mode pattern whereas adults showed a reliable floor-mode pattern. To summarize, while the whole-brain analysis for adults did not identify any region showing a reliable ceiling-mode AVA pattern, this pattern was found for adults when examining activity in regions defined from the children's results. In 3 other ROIs, adults and children showed floor-mode AVA patterns. Finally, in a large region encompassing the left IPL, adults showed a floor-mode AVA pattern whereas children showed a ceiling-mode pattern.

Discussion

The motivation for the current study was the idea that neural processing taking place during wakeful rest may be manifested in heterogeneous forms of RSA dynamics, and that these dynamics could vary depending on chronological age, clinical status, and intellectual maturity. Our examination was motivated by prior work showing that dynamics of the resting-state BOLD signal that are independent of low-frequency oscillations carry meaningful information (Hasson et al. 2009; Britz et al. 2010; Chang and Glover 2010; Garrett et al. 2010; Musso et al. 2010; Tagliazucchi et al. 2010; Van De Ville et al. 2010; Becker et al. 2011; Petridou et al. forthcoming, Tagliazucchi et al. 2012) and that BOLD variance is correlated with age (Garrett et al. 2010). Our first and most general main finding is that, during wakeful rest, the human brain demonstrates strong patterns of AVA activity that can be described by the relative stability of the amplitude of the local minima and maxima point of the BOLD signal.

We identified 2 types of AVA patterns: a “floor-mode” pattern showing greater variance for peaks of the BOLD time series than for its pits, and a “ceiling-mode” pattern demonstrating the opposite effect, with a greater variance for the pits of the time series. Secondly, we report that these AVA profiles show some differences but also several important commonalities between adults and children. Finally, we show that higher scores on an intelligence quotient measure in children are associated with greater similarity to adult-like activity patterns in certain regions.

Amplitude Variance Asymmetries as Indices to Stable Baseline States

The floor- and ceiling-mode resting-state fluctuation patterns that we find have not been previously documented. Beyond the potential functional significance of these profiles, which we discuss below, identifying these patterns allows addressing the issue of whether a relatively stable baseline exists in a given brain region and in doing so also quantify upward or downward departures from that baseline. Several methods have been developed for analyzing the BOLD signal to identify brain regions marked by transient, internally generated events (e.g., Tagliazucchi et al. 2010; Petridou et al. forthcoming), but to our knowledge, this is the first method to identify brain regions where departures from baseline, in the form of activations or deactivations, are exemplified by differences in variance between the sets of local-maxima and local-minima, thus providing more information about the nature of the fluctuation patterns themselves.

While perfusion-based neuroimaging allows measuring an absolute physiological baseline, BOLD data allow only for relative comparisons, and for this reason strong arguments have been made against assuming a functional baseline from BOLD data (Morcom and Fletcher 2007). Identifying AVA patterns in resting-state data contributes a partial solution to this debate as it allows identifying a functional baseline, albeit, by definition, only for cortical regions where deviations from baseline occur. The principle underlying this method is that the functional baseline is the steadier or more static state, and that departures from this stable state are marked by less stability (increased variance) as measured in the BOLD time series. When described from the viewpoint of dynamic systems outlined in Introduction, one can consider the baseline activity as a strong attractor, whereas departure instances from that attractor terminate in a more diverse range of endpoints.

The Nature of AVA

In contrast to methods that study RSA by identifying functional networks, the approach developed here does not assume temporal stationarity in the signal, but identifies regions with specific statistical activation features. It can therefore identify groups of cortical regions that show similar modes of operation without being synchronized or that change their degree of synchronization over time. However, because the AVA method does not rely on spatial constraints, it was important to determine the extent to which the identified regions form functional resting-state networks.

We performed such an analysis for the UNITN dataset for which we obtained independent physiological measurements that could be partialled out of the signal so to rule out alternative, physiologically related, explanations for connectivity. Our examination of functional connectivity between AVA regions documented the following results. First, overall functional connectivity within AVA regions was quite high. Secondly, by applying hierarchical clustering to single-participant connectivity maps, we found that the resulting clustering solutions were highly similar across participants, beyond what would be expected by chance. This suggests that connectivity patterns within AVA regions reflect intrinsic biological constraints that impose similar network structures across individuals, at a fine level of resolution. The above-chance similarity of the hierarchical clustering solution across individuals allowed us to construct a consensus solution that indicated aspects of the network structure shared across participants and to interpret these solutions in a meaningful way. It showed that the dynamics in AVA regions did not reflect a single source of spontaneous activity (though such cases are likely to exist, Scholvinck et al. 2010). Instead, these regions were partitioned into 3 sets according to sensory modality, with 1 set for bilateral premotor regions, one for bilateral temporal regions, and one for bilateral visual regions. This finding is consistent with research showing that homotopic regions demonstrate strong synchronization during wakeful rest (Stark et al. 2008). More importantly, this clustering analysis suggests that AVA patterns reflect the output of different sources of activity during rest, rather than a single source, consistent with recent theoretical models (Deco et al. 2011).

Potential Functional Significance of AVA

The current findings raise the question of what is the functional significance of activity in AVA regions and which of those activity patterns, if any, may be indexing cognitive activity during rest. A relatively surprising null result in this context was that AVA patterns were not found in regions associated with the DMN, a major network thought to mediate spontaneous thoughts, mind wandering, or intrinsic activity more generally (Raichle 2010). In what follows we discuss the possible sources of AVA patterns, and its absence from DMN regions.

For the adult data, the spatial distribution of AVA activity in visual and premotor regions resembles that seen in recent multimodal EEG/fMRI studies examining BOLD correlates of EEG bands. The bilateral premotor regions that we identify have been associated with the power envelope of rolandic alpha and beta EEG rhythms (10 and 20 Hz; Ritter et al. 2009), as well as oscillations in the gamma and theta rhythms during rest (Giraud et al. 2007). Equally important, in adults, the premotor regions showing AVA patterns do not overlap with motor regions where BOLD or EEG activity are modulated by actual movement of hand or foot (Yuan et al. 2010). This argues against the possibility that these transients were produced by participants' motion during the scan. Additionally, the distribution of these regions does not match that of sensorimotor networks identified via ICA during the resting states, as the latter also encompass medial aspects of the sensorimotor cortex (De Luca et al. 2006). We therefore suggest that the AVA pattern in ventral premotor regions does not index actual movement or maintenance of a general sensorimotor network. Instead, it may be that activity in these regions is associated with inhibition of movement, due to the experimental demand of participants to keep still during the study. Whether such inhibition of movement is associated with floor-mode or ceiling-mode activity patterns can be empirically examined in the future via controlled studies that examining AVA patterns in motor cortex during left/right motor tapping tasks or go/no-go tasks.

The large bilateral visual cluster we identify very strongly resembles regions where negative correlations between BOLD activity and global field power or relative power in alpha band have been repeatedly observed (Laufs et al. 2003, 2006; De Luca et al. 2006; Scheeringa et al. 2008, 2009, 2011; Sadaghiani et al. 2010; Yuan et al. 2010; Freyer et al. 2011). A recent study also suggests a positive correlation between BOLD activity and the phase of high-gamma in the visual cortex (Scheeringa et al. 2011). The negative correlation between BOLD and alpha power typically encompass a large extent of the dorsal visual stream, extending into the dorsal part of the intraparietal sulcus, and including also the lingual gyrus, strongly corresponding to our large bilateral visual cluster (an occipital–parietal pattern, Laufs et al. 2006; Sadaghiani et al. 2010).

This overlap between occipital and premotor regions showing AVA patterns and those identified by prior multimodal studies suggests that AVA patterns in these regions are related to departures from idling, often associated with desynchronization of the alpha rhythm (Freyer et al. 2011). Supporting this possibility is the finding that alpha band activity switches between a limited number of states signaling global phase synchronization, with states sometimes dwelled on for several seconds (Ito et al. 2007). If transients result from de-synchronization, they could be related to cognitive (information) processing during rest. However, for several reasons we are cautious to draw a strong link between the alpha activity per se and AVA. First, BOLD fluctuations during rest are best modeled as a combination of activity in different frequency bands rather than as being related to power in a particular band (Lu et al. 2007; Mantini et al. 2007; Scholvinck et al. 2010). The actual transfer function from neural activity to BOLD is likely a complex one, integrating over the entire frequency profile (Kilner et al. 2005; Rosa et al. 2010). Secondly, alpha power is correlated with a BOLD signal in fronto-parietal regions (Laufs et al. 2003), but we did not document an AVA profile in these regions, at least for adults.

The AVA patterns found in lateral temporal regions could be associated with various functions, including the generation of internal speech, imagery of auditory objects, or fluctuations in auditory attention. However, prior work suggests that modulation of auditory attention affects activity differently, with a strong impact on activity around primary and secondary auditory cortices (the supratemporal plane), including also SMAs, parietal regions, and prefrontal cortex (Voisin et al. 2006). This configuration is quite different from the one that we find in temporal regions. Engagement in auditory imagery of specific sounds activates aspects of STG and middle temporal gyri (MTG) that are more posterior than found here (Bunzeck et al. 2005). However, imagination of familiar songs during silence activates lateral temporal regions very similar to those that we identify (Kraemer et al. 2005). Thus, it is possible that AVA in lateral temporal regions may be associated with a process akin to “internal speech.”

When defining functional ROIs on the basis of children's data, we identified 1 region where adults showed a ceiling-mode activity pattern. This area encompassed both lateral and medial frontal regions. This mode of activity in frontal regions is consistent with findings from prior work documenting the relation between task and resting-state metabolic activity in this region using CMRO2 quantification (Lin et al. 2011). That work found that the left vmPFC demonstrated a unique CMRO2 deactivation pattern: it was characterized by task-induced deactivation for all participants, but interestingly, weaker deactivation for participants spending more time on task. That is, greater task involvement resulted in more baseline-like activity patterns. Taken together with the current finding of ceiling-mode effects in this region, this could suggest that, at least for some participants, vmPFC may often be in an “on” state of information processing during rest, and furthermore, that this state remains during strong task engagement. This conjecture is consistent with work showing that certain cortical regions demonstrate similar activity levels during rumination and highly demanding tasks (Dumontheil et al. 2010). However, the current study cannot empirically determine the precise source of difference between floor- and ceiling-mode activity patterns. The AVA is an unbiased method for identifying sources of spontaneous RSA. Cross-referencing AVA patterns with concurrent temporal data obtained via MEG or EEG will lead to an understanding of the types of neural fluctuations that differentiate the 2 patterns.

We expected the analysis of AVA patterns to identify DMN regions, but did not find this pattern. Our prediction was based on prior work (Golland et al. 2007), showing that several DMN regions, including the posterior cingulate and inferior parietal lobule, show activity patterns that are independent of environmental input. Specifically, that work demonstrated that BOLD fluctuations in these regions were not consistent across 2 presentations of the same movie segment;, i.e., these regions showed weak within-participant correlations, indicating that they were not driven by external input. That work also demonstrated that these “intrinsic” regions formed a functional network during rest and movie presentation. Our null findings for AVA patterns in the DMN are similar to those discussed previously (Laufs et al. 2006), noting that the DMN is not associated with BOLD correlates of rapid dynamics; specifically, BOLD fluctuations in DMN regions are not associated with either the occipital alpha or the fronto-parietal alpha patterns. Null effects however should always be interpreted with caution, as we discuss below (see Limitations and Directions for Future Work section).

Aspects of Development and Maturation

There were several findings relating AVA patterns to development and maturation. In several regions, there was a link between children's IQ and AVA patterns, and a direct contrast between the children and adults showed differences in AVA patterns in several regions. Before considering potential implications, we note that the results of the direct contrast between adults and children could be related to the fact that the children and adults were scanned at different MR facilities. Different scanning sites could vary on MR parameters that would affect sensitivity to AVA (Greve et al. 2010). These include differences in B0 distortions that could cause differential registration to common space and impact signal quality, or systematic differences in intensity that affect temporal SNR and thus sensitivity to AVA patterns. For these reasons, the direct comparison of adult and children data awaits a definitive replication at a single site. However, as we discuss below, the differences in patterns are consistent with prior work on differences between adult and children RSA generators.

The children showed a dominant ceiling-mode AVA pattern in practically the entire frontal cortex as well as inferior parietal regions bilaterally. These results may reflect well-documented differences in low-frequency oscillations in children versus adults. It is known that alpha, beta, and slow-wave power measured from frontal, parietal, and temporal regions weakens with development, and is accompanied by a reduction in gray matter volume (Whitford et al. 2007). Michels et al. (2011) compared EEG power in different bands between children (mean = 10 years of age) and adults, and showed that the strongest differences are found in dorsal pre- and post-central gyri as well as the medial and lateral frontal cortex. Interestingly, the areas where the differences in AVA patterns were strongest for adults and children are ones whose functional connectivity is stronger for adults than for children (Fair, Dosenbach et al. 2007). This may suggest that the reduction in spontaneous activity of the sort documented here gives rise to the long-range connectivity between IPL and dorsolateral prefrontal cortex that is more dominant in adults.

In the current study, some clusters where AVA signatures were found for adults were not found for children. It is possible, however, that this dissimilarity could be due to MR site-effects or to the statistical thresholds used. It is therefore important to not overlook the similarity in profiles found for adults and children: when adults' data were used as functional ROIs, signatures of adult-like AVA patterns were found in the children's data, and an overlap analysis showed moderate agreement between sets of regions showing floor-mode profiles for adults and children. Specifically, 2 occipital clusters where AVA was found for adults showed a floor-mode pattern in children when defined as functional ROIs, and in both of these clusters, higher IQ was associated with a higher σ2(peak)/ σ2(pits) ratio. This suggests that changes in AVA are linked to development or intellectual maturation and that AVA patterns in certain regions may show gradual changes over time. More generally, this overlap suggests that beyond developmental differences, there exist important similarities in AVA patterns between adults and children. On the basis of these findings, we suggest that the lack of AVA should not be treated as a hallmark of development. Instead, in fronto-parietal regions, development is associated with a shift toward unity, whereas other regions, e.g., occipital cortex, demonstrate an AVA above unity from childhood through adulthood. It may be that the transition from a ceiling-mode profile toward unity in fronto-parietal regions indicates that, for children, these regions mediate a tonic state of high activity with few departures from that state, potentially related to less-frequent activity bursts, or reduced power in lower frequencies. Evaluating these questions using EEG or MEG is an interesting question for future work.

The potential relation between development and the temporal features of neural activity has been the focus of several recent studies (see McIntosh et al. 2010, for recent review). Using MEG, these have quantified the degree of complexity in short response-locked epochs using MSE (Costa et al. 2002) and consistently found that activity in these epochs becomes more complex with age. This increase in complexity with age has generally been interpreted as a transition to more complex activity dynamics with development. To the extent that AVA departures from unity are interpreted as more complex patterns, the current findings indicate the converse; that in regions that differentiate adults from children, it is adults who show simpler dynamics. We note several essential differences between the prior and current work that could account for these apparent differences. First, in contrast to the aforementioned MEG work, the current study deals with endogenous, non-task-linked BOLD signals, devoid of any external driving input or decision-making. Thus, the actual cognitive process targeted is different. Secondly, the temporal length of the fMRI data analyzed in our work is 2–3 orders of magnitudes larger than the short epochs analyzed in that work, which often analyzed epochs shorter than 2 s. Therefore, conclusions derived from that developmental work on within-epoch signal complexity may not easily generalize over task, time scale, and type of measure and for this reason may not offer a direct basis for comparison. Furthermore, we do not see the 2 sets of results as necessarily divergent; it could be that MSE methods, if applied to resting-state fMRI data of appropriate length, would support the current findings. Our work is however consistent with fMRI studies, showing that increased BOLD variance is linked to behavioral indices: it is associated with better task performance (Garrett et al. 2011), and its sensitivity to task demands further correlates with age and task performance (Garrett et al. forthcoming).

Relation to Other Techniques for Analyzing Spontaneous Activity

Several features of the AVA analysis separate it from prior work on analysis of transient activity during rest. First, in contrast to several prior methods for examining variance structure of the BOLD signal (Morgan et al. 2008; Morgan and Gore 2009; Tagliazucchi et al. 2010), our work is not based on identification of (infrequent) spiking events and their spatial co-occurrence. Instead, it implements a holistic assessment of the entire time series after spike removal. That is, whereas spike analysis disregards time points that are not spikes, such data form valid inputs for our study and for this reason the data were purposefully de-spiked prior to analysis. Secondly, the method allows for nonstationarity in the time series that is manifested by slow shifts in mean and variance over time. Such slow shifts do not impact the sensitivity of the AVA analysis since local minima and maxima occur in a local context. As shown for the UNITN dataset, voxels demonstrating reliable AVA shifted between local minima and maxima quite often with modal peak-to-peak timings of 6 or 7.5 s in different clusters, and with a not-insignificant number appearing with 4.5 s of the prior peak. Thirdly, the results of the AVA analysis demonstrate high precision, in that the regions identified are structured in terms of anatomical systems. In contrast, the analysis of spike co-occurrence, while being able to reveal anatomical systems, may suffer from lower precision since it can identify extremely large sections of cortex where artifacts cause spike co-occurrence (Morgan and Gore 2009). Finally and most importantly, spike analyses are orthogonal to the issue of whether the spontaneous activity reflects floor- or ceiling-mode activity.

In this context it is important to consider more generally what type of additional information is offered by AVA in relation to other existing methods for quantifying time series complexity or disorder, and what are its relative applicative advantages. As mentioned above, MSE methods have been successfully used for quantifying the degree of disorder in highly complex signals, including MEG and EEG responses. These methods quantify the degree of disorder in different temporal scales and can potentially provide information regarding neural sources of signal fluctuations. However, to the best of our knowledge, MSE has not been applied to BOLD data, probably due to technical and interpretive hurdles: a) it necessitates very large time series of a sort not easily obtained in neuroimaging studies, so its applicability to shorter BOLD time series is unclear; b) MSE typically quantifies disorder in each temporal scale using sample entropy (SE). When applied to BOLD data, SE would partially reflect the degree of temporal autocorrelation in the signal, and could therefore be affected by factors impacting BOLD autocorrelation, including potential physiological impacts on the BOLD response.

By focusing on rapid fluctuations of local minima and maxima, AVA can be applied to short BOLD time series, as shown in its application to the 150 s UNITN data. This is an important factor when acquiring longer resting scans is not possible. AVA is also likely to accommodate impacts of smoothing, as differences in degree of smoothness due to non-neural factors may reduce the number of peaks/pits but will not necessarily affect their degree of variance asymmetry. Finally, a central aim of AVA is to characterize dynamics by identifying regions showing floor-mode or ceiling-mode activity patterns, as these may indicate different cognitive functions or particular baseline metabolic profiles that could be identified in future work. Entropy measures cannot address this issue as they will generate the same estimate to a time series y and to a time series [−1 × y].

Variance asymmetry also offers orthogonal information with respect to methods quantifying the degree of randomness in a signal, such as those evaluating the degree of diversity in signal values as quantified via Shannon's entropy (see de Araujo et al. 2003, for application to BOLD time series). Shannon's entropy increases with the uniformity of the distribution and with the number of unique elements in the input. Two inputs can have very different Shannon entropy values but the same AVA, as long as the relative variance of peaks and pits is similar in both. Other forms of entropy measures quantify the degree of sequential regularity or autocorrelation in the input (e.g., mutual information, sample entropy, conditional entropy). Inputs with sequential regularities such as sine-wave patterns depart from randomness due to their autocorrelation, but as mentioned, AVA may be identical for a random and completely sine-like pattern. To summarize, AVA does not measure the randomness of a series.

Recent work (Petridou et al. forthcoming) examining nonstationary aspects of RSA has identified several networks demonstrating spontaneous BOLD events during rest, whose magnitude appears to be similar to task-evoked responses. Using a paradigm-free sparse-regression method (SPFM; see Caballero-Gaudes et al. 2011, forthcoming, for technical development), the authors found spontaneous events in several resting-state networks and showed that these events account for nonstationary connectivity patterns within these networks. The networks identified in that work partially overlap with those identified by our AVA method, including occipital, pre- and postcentral gyri (note that Caballero-Gaudes et al. used a limited field of view). However, the nature of events captured by sparse regression appears quite distinct from the fluctuations associated with AVA. Several of the networks identified in that work showed as few as 2 events in a 10-min scanning session, and none exceeded 10 events for any participant. In contrast, the local-maxima and -minima constituting our data target a feature derived from more rapid fluctuations, with peaks appearing, on average, every 8–10 s in areas showing statistically significant AVA on the single-participant level. Unlike SPFM, AVA is not intended for the detection of single peaks within the data or their timing. Instead, the method collapses over time to obtain a single metric and can also be statistically significant even when the time series reflects numerous weak spikes. This aspect of AVA is complementary to SPFM as the latter models the BOLD time series using the most parsimonious positioning of assumed event timings. As noted by Caballero-Gaudes et al., SPFM performs optimally when intrinsic events are relatively infrequent, but less so if the BOLD time series reflects numerous generated events in tight procession (see, e.g., Skipper et al. 2009 for example of such cases during naturalistic movie observation). In such cases, combining the two methods can offer potential advantages as voxels demonstrating AVA patterns may be an interesting subset of voxel space for which it would be possible to implement appropriately calibrated paradigm-free analyses.

We finally consider the relation between AVA and the variance of the BOLD time series, a measure used in recent work to characterize developmental progression (Garrett et al. 2010). As emphasized in Potential Relation of AVA Patterns to Variance of Time Series section, the 2 measures provide complimentary and independent information, with factors that affect variance (e.g., amplitude modulation) not affecting AVA. Thus, AVA focuses on a specific feature of time series dynamics and is insensitive to the absolute magnitude of signal fluctuations. When studying features of long BOLD resting-state series, combining AVA and raw variance measures may be useful, as variance alone may at least partially reflect power in low-frequency bands. Studying differences between BOLD dynamics in different-aged populations or between clinical and nonclinical groups using variance measures could be sensitive to confounds resulting from vasculature effects on BOLD, which may indicate apparent differences in variance due to non-neural factors. BOLD quantifications using variance will also be sensitive to comparisons across MR scanners with different field strengths, as stronger fields can measure stronger fluctuations. In such cases, using AVA can yield useful information due to its independence from the absolute magnitude of fluctuations.

Limitations and Directions for Future Work

Given the initial treatment of variance asymmetry in the current study, our work is limited in several aspects that should be clarified and tended to in future studies. These have to do with specific issues of the current study, as well as more general concerns.

With respect to the child population, we note that the cohort of nonclinical children had a mean IQ of 119, which is more than 1 standard deviation above the population mean. This could bias the interpretation of the result since the children cohort studied is not representative of the random population and may reflect a particularly intelligent subgroup.

Some analyses decisions could also have impacted the results. From the UNITN dataset, we chose to remove independently acquired autonomic nervous system (ANS) cardiac and respiration data, and from the other datasets we used proxy information from CSF time series. As we noted in prior work (Iacovella and Hasson 2011), while correlates of ANS indices are typically considered physiological noise, BOLD time series in certain regions may correlate with physiological fluctuations due to the fact that those regions are involved in information processing that is tightly linked to immediate fluctuations in ANS activity rather than via physiological effects on the BOLD signal. In that case, removing BOLD variance explained by ANS fluctuations could remove meaningful activity components from these regions. Future work could attempt to independently identify regions where BOLD fluctuations correlate with ANS indices, and evaluate their AVA metrics before and after the removal of ANS-related BOLD variance.

An additional limitation of the current work pertains to the analysis of data from different sites. Because we wanted to average AVA maps within each participant, we obtained datasets where more than 1 resting scan was available for each participant. Consequently, the children and adult data were collected on different scanners using different protocols. Ideally, both datasets should have been optimally collected at the same site and future work should collect data across age groups at the same site. Scanning different populations raises the susceptibility of signal quality differences or differential effect of motion-induced artifacts due to different during-scan movement patterns (Power et al. 2012). For adults, AVA clusters partitioned into 3 different networks. For children, AVA was found in medial regions relatively insensitive to motion, as well as in fronto-parietal regions. These patterns do not rule out an effect of motion (which, although corrected for, surely exists in some form), but do argue against a decisive role for motion as the major cause of these patterns.

When evaluating regions where AVA differed for adults and children, we found that the differences were driven by differences in σ2(peaks) or σ2(pits). Examining features of the peaks and pits separately could prove a fruitful direction for future work. Temporal features of between-peak intervals as well as the dynamics of the amplitude changes themselves may carry important information and could be analyzed via techniques such as phase-space analysis.

Finally, we consider limitations and constraints of the method itself that result from its formal features or its reliance on the BOLD signal in the current work. First, while AVA is intended to identify a type of nonrandom type of fMRI signal fluctuation, it identifies a particular subset of nonrandom fluctuation. It does not identify brain areas exhibiting symmetric bi-stable states during rest (Freyer et al. 2009; Freyer et al. 2011); states in which transitions occur symmetrically between 2 equally strong attractors. Similarly, regions that tend to stay at a modal state, with occasional spikes of activation and deactivation of similar magnitude, will not be identified as the variance of maxima and minima would be similar. Cases where departures are just positive or just negative will also not be identified if departures are consistently of similar magnitude. Such infrequent events have been documented in the DMN (Petridou et al. forthcoming). Well-behaved sine-like oscillations similarly do no show AVA. Finally, if a regions' RSA profile is nonstationary and oscillates equally between floor-mode and ceiling-mode patterns over time, aggregation over the entire time-series would obfuscate this pattern from an AVA analysis. Computing AVA along the time series in appropriate sliding windows would solve this particular issue. To summarize, there exist several profiles of spontaneous RSA that would not be identified with AVA. AVA should be considered complementary to methods such as those mentioned above that can identify these profiles.

AVA analysis of the BOLD signal as performed in the current study may suffer from false negatives (misses) due to the nonlinearity of the BOLD response. Numerous studies examining responses to exogenous stimuli show that BOLD scales nonlinearly with stimulus duration (e.g., Miller et al. 2001). Assuming that AVA reflects neural activity bursts of different lengths, nonlinearity could lead to reduced manifestation of AVA in the resulting BOLD signal or potential suppression of observable AVA altogether. This is because nonlinearity can make the BOLD amplitude of a strong neural response appear similar to that of a weaker response. Combining the amplitude parameter with a temporal signature (e.g., the relative duration in the “up” or “down” states) may allow more sophisticated differentiation between states than afforded by the identification of peaks and pits alone. A related question is the extent to which BOLD transients reflect neuronal, vascular, or hemodynamic effects. BOLD transients may reflect a hemodynamic, non-neural origin, e.g., due to temporal offsets between CBF and CBV changes (e.g., Obata et al. 2004). Whether different CBF:CBV offsets can lead to BOLD AVA, and if so, what is the nature of the relationship between the amplitude of physiologically induced transients and neural activity is an issue for future work. Analyses of perfusion data could potentially address this issue by directly studying cerebral blood flow patterns.

In relation to this, it may be that differences in the nonlinearity of the BOLD response itself may account for differences in AVA seen between adults and children. This is a general issue underlying any fMRI analysis comparing different population groups that differ by age, clinical status, or any other parameter that could be related to flow dynamics and vascular compliance. For instance, if the nonlinearity profile of the BOLD response differs across population, saturates more quickly, or has a more limited dynamic range, then this could result in more uniform activity for one group than another, thus inducing differences in AVA. Indeed, the impact of non-neural factors is a generic concern when comparing populations: these factors can drive differences in functional connectivity, determine differences in BOLD percent signal change, or affect the magnitude of BOLD resting-state variance. The use of nonhemodynamic measures is therefore an important direction for future work.

Summary

We presented novel principles underlying RSA in adults in children, derived from a method for analysis of RSA that is focused on nonoscillatory, rapid aspects of BOLD dynamics. The current work shows that in adults, sensory regions show spontaneous, floor-mode activity patterns with greater variance of local-maxima than local-minima in the signal, and that these spontaneous patterns are driven by 3 separate processes rather than a single one. A comparative analysis with children's data demonstrated qualitative developmental differences, but indicated that even at a young age, signs of adult-like activity patterns can be found, and that these are correlated with IQ. The relation of such spontaneous activity patterns to behavioral measures and other clinical states are interesting avenues for future work.

Supplementary Material

Supplementary material can be found at: http://www.cercor.oxfordjournals.org/.

Funding

This research has received funding from the European Research Council under the 7th framework starting grant program (ERC-STG #263318 to U.H) and from Savings Bank Foundation (Fondazione Cassa di Risparmio) of Trento and Rovereto.

Notes

Conflict of Interest: None declared.

References

Achard
S
Salvador
R
Whitcher
B
Suckling
J
Bullmore
E
A resilient, low-frequency, small-world human brain functional network with highly connected association cortical hubs
J Neurosci
 , 
2006
, vol. 
26
 (pg. 
63
-
72
)
Argall
BD
Saad
ZS
Beauchamp
MS
Simplified intersubject averaging on the cortical surface using SUMA
Hum Brain Mapp
 , 
2005
, vol. 
27
 (pg. 
14
-
27
)
Becker
R
Reinacher
M
Freyer
F
Villringer
A
Ritter
P
How ongoing neuronal oscillations account for evoked fMRI variability
J Neurosci
 , 
2011
, vol. 
31
 (pg. 
11016
-
11027
)
Birn
RM
Diamond
JB
Smith
MA
Bandettini
PA
Separating respiratory-variation-related fluctuations from neuronal-activity-related fluctuations in fMRI
NeuroImage
 , 
2006
, vol. 
31
 (pg. 
1536
-
1548
)
Biswal
B
Yetkin
F
Haughton
V
Hyde
J
Functional connectivity in the motor cortex of resting human brain using echo-planar MRI
Magn Reson Med
 , 
1995
, vol. 
34
 (pg. 
537
-
541
)
Britz
J
van de Ville
D
Michel
CM
BOLD correlates of EEG topography reveal rapid resting-state network dynamics
NeuroImage
 , 
2010
, vol. 
52
 (pg. 
1162
-
1170
)
Buckner
RL
Andrews-Hanna
JR
Schacter
DL
The brain's default network: anatomy, function, and relevance to disease
Ann N Y Acad Sci
 , 
2008
, vol. 
1124
 (pg. 
1
-
38
)
Bunzeck
N
Wuestenberg
T
Lutz
K
Heinze
HJ
Jancke
L
Scanning silence: mental imagery of complex sounds
NeuroImage
 , 
2005
, vol. 
26
 (pg. 
1119
-
1127
)
Caballero-Gaudes
C
Petridou
N
Dryden
IL
Bai
L
Francis
ST
Gowland
PA
Detection and characterization of single-trial fMRI BOLD responses: paradigm free mapping
Hum Brain Mapp
 , 
2011
, vol. 
32
 (pg. 
1400
-
1418
)
Caballero-Gaudes
C
Petridou
N
Francis
ST
Dryden
IL
Gowland
PA
Paradigm free mapping with sparse regression automatically detects single-trial functional magnetic resonance imaging blood oxygenation level dependent responses
Hum Brain Mapp
 , 
forthcoming
Chang
C
Glover
G
Time-frequency dynamics of resting-state brain connectivity measured with fMRI
NeuroImage
 , 
2010
, vol. 
50
 (pg. 
81
-
98
)
Chialvo
DR
Emergent complex neural dynamics
Nat Phys
 , 
2010
, vol. 
6
 (pg. 
744
-
750
)
Chung
MK
Robbins
SM
Dalton
KM
Davidson
RJ
Alexander
AL
Evans
AC
Cortical thickness analysis in autism with heat kernel smoothing
NeuroImage
 , 
2005
, vol. 
25
 (pg. 
1256
-
1265
)
Cordes
D
Haughton
VM
Arfanakis
K
Wendt
GJ
Turski
PA
Moritz
CH
Quigley
MA
Meyerand
ME
Mapping functionally related regions of brain with functional connectivity MR imaging
Am J Neuroradiol
 , 
2000
, vol. 
21
 (pg. 
1636
-
1644
)
Costa
M
Goldberger
AL
Peng
CK
Multiscale entropy analysis of complex physiologic time series
Phys Rev Lett
 , 
2002
, vol. 
89
 pg. 
068102
 
Cox
RW
AFNI: software for analysis and visualization of functional magnetic resonance neuroimages
Comput Biomed Res
 , 
1996
, vol. 
29
 (pg. 
162
-
173
)
de Araujo
DB
Tedeschi
W
Santos
AC
Elias
J
Jr
Neves
UPC
Baffa
O
Shannon entropy applied to the analysis of event-related fMRI time series
NeuroImage
 , 
2003
, vol. 
20
 (pg. 
311
-
317
)
Deco
G
Jirsa
VK
McIntosh
AR
Emerging concepts for the dynamical organization of resting-state activity in the brain
Nat Rev Neurosci
 , 
2011
, vol. 
12
 (pg. 
43
-
56
)
Deco
G
Jirsa
V
McIntosh
AR
Sporns
O
Kötter
R
Key role of coupling, delay, and noise in resting brain fluctuations
Proc Natl Acad Sci USA
 , 
2009
, vol. 
106
 (pg. 
10302
-
10307
)
Deco
G
Rolls
ET
Albantakis
L
Romo
R
Brain mechanisms for perceptual and reward-related decision-making
Progress Neurobiol
 , 
forthcoming
De Luca
M
Beckmann
CF
De Stefano
N
Matthews
PM
Smith
SM
fMRI resting state networks define distinct modes of long-distance interactions in the human brain
NeuroImage
 , 
2006
, vol. 
29
 (pg. 
1359
-
1367
)
Desai
R
Liebenthal
E
Possing
ET
Waldron
E
Binder
JR
Volumetric vs. surface-based alignment for localization of auditory cortex activation
NeuroImage
 , 
2005
, vol. 
26
 (pg. 
1019
-
1029
)
Deshpande
G
Laconte
S
Peltier
S
Hu
XP
Tissue specificity of nonlinear dynamics in baseline fMRI
Magn Reson Med
 , 
2006
, vol. 
55
 (pg. 
626
-
632
)
Dice
LR
Measures of the amount of ecologic association between species
Ecology
 , 
1945
, vol. 
26
 (pg. 
297
-
302
)
Dumontheil
I
Gilbert
SJ
Frith
CD
Burgess
PW
Recruitment of lateral rostral prefrontal cortex in spontaneous and task-related thoughts
Q J Exp Psychol
 , 
2010
, vol. 
63
 (pg. 
1740
-
1756
)
Durstewitz
D
Deco
G
Computational significance of transient dynamics in cortical networks
Eur J Neurosci
 , 
2008
, vol. 
27
 (pg. 
217
-
227
)
Fair
D
Dosenbach
N
Church
J
Cohen
A
Brahmbhatt
S
Miezin
F
Barch
D
Raichle
M
Petersen
S
Schlaggar
B
Development of distinct control networks through segregation and integration
Proc Natl Acad Sci USA
 , 
2007
, vol. 
104
 (pg. 
13507
-
13512
)
Fair
DA
Bathula
D
Mills
KL
Costa Dias
TG
Blythe
MS
Zhang
D
Snyder
AZ
Raichle
ME
Stevens
AA
Nigg
JT
, et al.  . 
Maturing thalamocortical functional connectivity across development
Front Neurosci
 , 
2010
, vol. 
4
 pg. 
10
 
Fair
DA
Cohen
AL
Dosenbach
NUF
Church
JA
Miezin
FM
Barch
DM
Raichle
ME
Petersen
SE
Schlaggar
BL
The maturing architecture of the brain's default network
Proc Natl Acad Sci USA
 , 
2008
, vol. 
105
 (pg. 
4028
-
4032
)
Fair
DA
Posner
J
Nagel
BJ
Bathula
D
Costa Dias
TG
Mills
KL
Blythe
MS
Giwa
A
Schmitt
CF
Nigg
JT
Atypical default network connectivity in youth with ADHD
Biol Psychiatr
 , 
2010
, vol. 
68
 (pg. 
1084
-
1091
)
Fair
DA
Schlaggar
BL
Cohen
AL
Miezin
FM
Dosenbach
NUF
Wenger
KK
Fox
MD
Snyder
A
Raichle
ME
Petersen
SE
A method for using blocked and event-related fMRI data to study “resting state” functional connectivity
NeuroImage
 , 
2007
, vol. 
35
 (pg. 
396
-
405
)
Fischl
B
Sereno
MI
Dale
AM
Cortical surface-based analysis II: inflation, flattening, and a surface-based coordinate system
NeuroImage
 , 
1999
, vol. 
9
 (pg. 
195
-
207
)
Fischl
B
Sereno
MI
Tootell
RB
Dale
AM
High-resolution intersubject averaging and a coordinate system for the cortical surface
Hum Brain Mapp
 , 
1999
, vol. 
8
 (pg. 
272
-
284
)
Forman
SD
Cohen
JD
Fitzgerald
M
Eddy
WF
Mintun
MA
Noll
DC
Improved assessment of significant activation in functional magnetic resonance imaging. fMRI: use of a cluster-size threshold
Magn Reson Med
 , 
1995
, vol. 
33
 (pg. 
636
-
647
)
Fox
J
Applied regression analysis, linear models, and related methods
 , 
1997
Thousand Oaks, CA: Sage Publications, Inc
Fox
MD
Corbetta
M
Snyder
AZ
Vincent
JL
Raichle
ME
Spontaneous neuronal activity distinguishes human dorsal and ventral attention systems
Proc Natl Acad Sci USA
 , 
2006
, vol. 
103
 (pg. 
10046
-
10051
)
Freeman
WJ
How brains make up their minds
 , 
2001
New York, NY
Columbia University Press
Freyer
F
Aquino
K
Robinson
PA
Ritter
P
Breakspear
M
Bistability and non-Gaussian fluctuations in spontaneous cortical activity
J Neurosci
 , 
2009
, vol. 
29
 (pg. 
8512
-
8524
)
Freyer
F
Roberts
J
Becker
R
Robinson
P
Ritter
P
Breakspear
M
Biophysical mechanisms of multistability in resting-state cortical rhythms
J Neurosci
 , 
2011
, vol. 
31
 (pg. 
6353
-
6361
)
Garrett
DD
Kovacevic
N
McIntosh
AR
Grady
CL
Blood oxygen level-dependent signal variability is more than just noise
J Neurosci
 , 
2010
, vol. 
30
 (pg. 
4914
-
4921
)
Garrett
DD
Kovacevic
N
McIntosh
AR
Grady
CL
The importance of being variable
J Neurosci
 , 
2011
, vol. 
31
 (pg. 
4496
-
503
)
Garrett
DD
Kovacevic
N
McIntosh
AR
Grady
CL
The modulation of BOLD variability between cognitive states varies by age and processing speed
Cereb Cortex
 , 
forthcoming
Ghosh
A
Rho
Y
McIntosh
AR
Kötter
R
Jirsa
VK
Noise during rest enables the exploration of the brain's dynamic repertoire
PLoS Comput Biol
 , 
2008
, vol. 
4
 pg. 
e1000196
 
Giraud
AL
Kleinschmidt
A
Poeppel
D
Lund
TE
Frackowiak
RSJ
Laufs
H
Endogenous cortical rhythms determine cerebral specialization for speech perception and production
Neuron
 , 
2007
, vol. 
56
 (pg. 
1127
-
1134
)
Glahn
DC
Winkler
aM
Kochunov
P
Almasy
L
Duggirala
R
Carless
Ma
Curran
JC
Olvera
RL
Laird
aR
Smith
SM
, et al.  . 
Genetic control over the resting brain
Proc Natl Acad Sci USA
 , 
2010
, vol. 
107
 (pg. 
1223
-
1228
)
Glover
GH
Li
TQ
Ress
D
Image-based method for retrospective correction of physiological motion effects in fMRI: RETROICOR
Magn Reson Med
 , 
2000
, vol. 
44
 (pg. 
162
-
167
)
Golland
Y
Bentin
S
Gelbard
H
Benjamini
Y
Heller
R
Nir
Y
Hasson
U
Malach
R
Extrinsic and intrinsic systems in the posterior cortex of the human brain revealed during natural sensory stimulation
Cereb Cortex
 , 
2007
, vol. 
17
 (pg. 
766
-
777
)
Greve
D
Mueller
B
Brown
G
Liu
T
Glover
GH
FBIRN
Processing methods to reduce intersite variability in fMRI
Proceedings of the 16th Annual Meeting Organization for Human Brain Mapping
Barcelona, Spain
 
June, 2010 (abstract 1318)
Hampson
M
Driesen
NR
Skudlarski
P
Gore
JC
Constable
RT
Brain connectivity related to working memory performance
J Neurosci
 , 
2006
, vol. 
26
 (pg. 
13338
-
13343
)
Hasson
U
Nusbaum
HC
Small
SL
Task-dependent organization of brain regions active during rest
Proc Natl Acad Sci USA
 , 
2009
, vol. 
106
 (pg. 
10841
-
10846
)
Hornik
K
A CLUE for CLUster ensembles
J Statist Software
 , 
2005
, vol. 
14
 (pg. 
1
-
25
)
Hunter
MD
Eickhoff
SB
Miller
TWR
Farrow
TFD
Wilkinson
ID
Woodruff
PWR
Neural activity in speech-sensitive auditory cortex during silence
Proc Natl Acad Sci USA
 , 
2006
, vol. 
103
 (pg. 
189
-
194
)
Iacovella
V
Hasson
U
The relationship between BOLD signal and autonomic nervous system functions: implications for processing of “physiological noise”
Magn Reson Imaging
 , 
2011
, vol. 
29
 (pg. 
1338
-
1345
)
Ibanez
F
Grosjean
P
Etienne
M
pastecs
 , 
2009
1.3–11 Edition
Ito
J
Nikolaev
AR
van Leeuwen
C
Dynamics of spontaneous transitions between global brain states
Hum Brain Mapp
 , 
2007
, vol. 
28
 (pg. 
904
-
913
)
Jenkinson
M
Smith
SM
A global optimisation method for robust affine registration of brain images
Medical Image Analysis
 , 
2001
, vol. 
5
 (pg. 
143
-
156
)
Kilner
JM
Mattout
J
Henson
R
Friston
KJ
Hemodynamic correlates of EEG: a heuristic
NeuroImage
 , 
2005
, vol. 
28
 (pg. 
280
-
286
)
Kraemer
DJM
Macrae
CN
Green
AE
Kelley
WM
Musical imagery: sound of silence activates auditory cortex
Nature
 , 
2005
, vol. 
434
 pg. 
158
 
Laufs
H
Holt
JL
Elfont
R
Krams
M
Paul
JS
Krakow
K
Kleinschmidt
A
Where the BOLD signal goes when alpha EEG leaves
NeuroImage
 , 
2006
, vol. 
31
 (pg. 
1408
-
1418
)
Laufs
H
Krakow
K
Sterzer
P
Eger
E
Beyerle
A
Salek-Haddadi
A
Kleinschmidt
A
Electroencephalographic signatures of attentional and cognitive default modes in spontaneous brain activity fluctuations at rest
Proc Natl Acad Sci USA
 , 
2003
, vol. 
100
 (pg. 
11053
-
11058
)
Lin
P
Hasson
U
Jovicich
J
Robinson
S
A neuronal basis for task-negative responses in the human brain
Cereb Cortex
 , 
2011
, vol. 
21
 (pg. 
821
-
830
)
Lu
H
Zuo
Y
Gu
H
Waltz
Ja
Zhan
W
Scholl
Ca
Rea
W
Yang
Y
Stein
Ea
Synchronized delta oscillations correlate with the resting-state functional MRI signal
Proc Natl Acad Sci USA
 , 
2007
, vol. 
104
 (pg. 
18265
-
18269
)
Mantini
D
Perrucci
MG
Del Gratta
C
Romani
GL
Corbetta
M
Electrophysiological signatures of resting state networks in the human brain
Proc Natl Acad Sci USA
 , 
2007
, vol. 
104
 (pg. 
13170
-
13175
)
Mazaheri
A
Jensen
O
Asymmetric amplitude modulations of brain oscillations generate slow evoked responses
J Neurosci
 , 
2008
, vol. 
28
 (pg. 
7781
-
7787
)
McDonnell
MD
Ward
LM
The benefits of noise in neural systems: bridging theory and experiment
Nat Rev Neurosci
 , 
2011
, vol. 
12
 (pg. 
415
-
426
)
McIntosh
AR
Kovacevic
N
Lippe
S
Garrett
D
Grady
C
Jirsa
V
The development of a noisy brain
Arch Ital Biol
 , 
2010
, vol. 
148
 (pg. 
323
-
37
)
Michels
L
Bucher
K
Brem
S
Halder
P
Lüchinger
R
Liechti
M
Martin
E
Jeanmonod
D
Kröll
J
Brandeis
D
Does greater low frequency EEG activity in normal immaturity and in children with epilepsy arise in the same neuronal network?
Brain Topogr
 , 
2011
, vol. 
24
 (pg. 
78
-
89
)
Miller
KL
Luh
WM
Liu
TT
Martinez
A
Obata
T
Wong
EC
Frank
LR
Buxton
RB
Nonlinear temporal dynamics of the cerebral blood flow response
Hum Brain Mapp
 , 
2001
, vol. 
13
 (pg. 
1
-
12
)
Morcom
AM
Fletcher
PC
Does the brain have a baseline? Why we should be resisting a rest
NeuroImage
 , 
2007
, vol. 
37
 (pg. 
1073
-
1082
)
Morgan
VL
Gore
JC
Detection of irregular, transient fMRI activity in normal controls using 2dTCA: comparison to event-related analysis using known timing
Hum Brain Mapp
 , 
2009
, vol. 
30
 (pg. 
3393
-
3405
)
Morgan
VL
Gore
JC
Szaflarski
JP
Temporal clustering analysis: what does it tell us about the resting state of the brain?
Med Sci Monitor
 , 
2008
, vol. 
14
 (pg. 
345
-
352
)
Murphy
K
Birn
RM
Handwerker
DA
Jones
TB
Bandettini
PA
The impact of global signal regression on resting state correlations: are anti-correlated networks introduced
NeuroImage
 , 
2009
, vol. 
44
 (pg. 
893
-
905
)
Musso
F
Brinkmeyer
J
Mobascher
A
Warbrick
T
Winterer
G
Spontaneous brain activity and EEG microstates. A novel EEG/fMRI analysis approach to explore resting-state networks
NeuroImage
 , 
2010
, vol. 
52
 (pg. 
1149
-
1161
)
Obata
T
Liu
TT
Miller
KL
Luh
W
Wong
EC
Frank
LR
Buxton
RB
Discrepancies between BOLD and flow dynamics in primary and supplementary motor areas: application of the balloon model to the interpretation of BOLD transients
NeuroImage
 , 
2004
, vol. 
21
 (pg. 
144
-
153
)
Petridou
N
Gaudes
CC
Dryden
IL
Francis
ST
Gowland
P
Periods of rest in fMRI contain individual spontaneous events which are related to slowly fluctuating spontaneous activity
Hum Brain Mapp
 , 
forthcoming
Power
JD
Barnes
KA
Snyder
AZ
Schlaggar
BL
Petersen
SE
Spurious but systematic correlations in functional connectivity MRI networks arise from subject motion
NeuroImage
 , 
2012
, vol. 
59
 (pg. 
2142
-
2154
)
Raichle
ME
Two views of brain function
Trends Cogn Sci
 , 
2010
, vol. 
14
 (pg. 
180
-
190
)
Ritter
P
Moosmann
M
Villringer
A
Rolandic alpha and beta EEG rhythms’ strengths are inversely related to fMRI BOLD signal in primary somatosensory and motor cortex
Human Brain Mapping
 , 
2009
, vol. 
30
 (pg. 
1168
-
1187
)
Robinson
S
Basso
G
Soldati
N
Sailer
U
Jovicich
J
Bruzzone
L
Kryspin-Exner
I
Bauer
H
Moser
E
A resting state network in the motor control circuit of the basal ganglia
BMC neuroscience
 , 
2009
, vol. 
10
 pg. 
137
 
Rosa
MJ
Kilner
J
Blankenburg
F
Josephs
O
Penny
W
Estimating the transfer function from neuronal activity to BOLD using simultaneous EEG-fMRI
NeuroImage
 , 
2010
, vol. 
49
 (pg. 
1496
-
1509
)
R_Team
R: A language and environment for statistical computing
2008
 
R Foundation for Statistical Computing, Vienna Austria
Saad
Z
Glen
D
Chen
G
Beauchamp
M
Desai
R
Cox
R
A new method for improving functional-to-structural alignment using local Pearson correlation
NeuroImage
 , 
2009
, vol. 
44
 (pg. 
839
-
848
)
Saad
ZS
Reynolds
RC
Argall
B
Japee
S
Cox
RW
SUMA: an interface for surface-based intra-and inter-subject analysis with AFNI
2004
IEEE International Symposium on Biomedical Imaging
Arlgton, VA
(pg. 
1510
-
1513
)
Sadaghiani
S
Scheeringa
R
Lehongre
K
Morillon
B
Giraud
AL
Kleinschmidt
A
Intrinsic connectivity networks, alpha oscillations, and tonic alertness: a simultaneous electroencephalography/functional magnetic resonance imaging study
J Neurosci
 , 
2010
, vol. 
30
 (pg. 
10243
-
10250
)
Scheeringa
R
Bastiaansen
MCM
Petersson
KM
Oostenveld
R
Norris
DG
Hagoort
P
Frontal theta EEG activity correlates negatively with the default mode network in resting state
Int J Psychophysiol
 , 
2008
, vol. 
67
 (pg. 
242
-
251
)
Scheeringa
R
Fries
P
Petersson
K-M
Oostenveld
R
Grothe
I
Norris
DG
Hagoort
P
Bastiaansen
MCM
Neuronal dynamics underlying high- and low-frequency EEG oscillations contribute independently to the human BOLD signal
Neuron
 , 
2011
, vol. 
69
 (pg. 
572
-
583
)
Scheeringa
R
Magnus
K
Oostenveld
R
Norris
DG
Hagoort
P
Bastiaansen
MCM
NeuroImage trial-by-trial coupling between EEG and BOLD identifies networks related to alpha and theta EEG power increases during working memory maintenance
NeuroImage
 , 
2009
, vol. 
44
 (pg. 
1224
-
1238
)
Scholvinck
ML
Maier
A
Ye
FQ
Duyn
JH
Leopold
DA
Neural basis of global resting-state fMRI activity
Proc Natl Acad Sci USA
 , 
2010
, vol. 
107
 (pg. 
10238
-
10243
)
Shehzad
Z
Kelly
A
Reiss
P
Gee
D
Gotimer
K
Uddin
L
Lee
S
Margulies
D
Roy
A
Biswal
B
, et al.  . 
The resting brain: unconstrained yet reliable
Cereb Cortex
 , 
2009
, vol. 
19
 (pg. 
2209
-
2229
)
Skipper
JI
Goldin-Meadow
S
Nusbaum
HC
Small
SL
Gestures orchestrate brain networks for language understanding
Curr Biol
 , 
2009
, vol. 
19
 (pg. 
661
-
667
)
Smith
M
Fast robust automated brain extraction
Hum Brain Mapp
 , 
2002
, vol. 
17
 (pg. 
143
-
155
)
Spivey
MJ
Anderson
SE
Dale
R
The phase transition in human cognition
New Math Nat Computation
 , 
2009
, vol. 
5
 (pg. 
197
-
220
)
Sporns
O
The non- random brain: efficiency, economy, and complex dynamics
Front Comput Neurosci
 , 
2011
, vol. 
5
 pg. 
5
 
Stam
CJ
Pijn
J
Pritchard
W
Reliable detection of nonlinearity in experimental time series with strong periodic components
Phys D
 , 
1998
, vol. 
112
 (pg. 
361
-
380
)
Stam
CJ
Pijn
JP
Suffczynski
P
Silva
FHLd
Dynamics of the human alpha rhythm: evidence for non-linearity?
Clinical Neurophysiology
 , 
1999
, vol. 
110
 (pg. 
1801
-
1813
)
Stark
DE
Margulies
DS
Shehzad
ZE
Reiss
P
Kelly
AM
Uddin
LQ
Gee
DG
Roy
AK
Banich
MT
Castellanos
FX
Regional variation in interhemispheric coordination of intrinsic hemodynamic fluctuations
J Neurosci
 , 
2008
, vol. 
28
 (pg. 
13754
-
13764
)
Tagliazucchi
E
Balenzuela
P
Fraiman
D
Chialvo
DR
Criticality in large-scale brain fMRI dynamics unveiled by a novel point process analysis
Front Fractal Physiol
 , 
2012
, vol. 
3
 pg. 
15
 
Tagliazucchi
E
Balenzuela
P
Fraiman
D
Montoya
P
Chialvo
DR
Spontaneous BOLD event triggered averages for estimating functional connectivity at resting state
Neurosci Lett
 , 
2010
, vol. 
482
 (pg. 
158
-
163
)
Van De Ville
D
Britz
J
Michel
CM
EEG microstate sequences in healthy humans at rest reveal scale-free dynamics
Proc Natl Acad Sci USA
 , 
2010
, vol. 
107
 (pg. 
18179
-
18184
)
Voisin
J
Bidet-Caulet
A
Bertrand
O
Fonlupt
P
Listening in silence activates auditory areas: a functional magnetic resonance imaging study
J Neurosci
 , 
2006
, vol. 
26
 (pg. 
273
-
278
)
Wang
K
Jiang
T
Yu
C
Tian
L
Li
J
Liu
Y
Zhou
Y
Xu
L
Song
M
Li
K
Spontaneous activity associated with primary visual cortex: a resting-state FMRI study
Cereb Cortex
 , 
2008
, vol. 
18
 (pg. 
697
-
704
)
Whitford
TJ
Rennie
CJ
Grieve
SM
Clark
CR
Gordon
E
Williams
LM
Brain maturation in adolescence: concurrent changes in neuroanatomy and neurophysiology
Hum Brain Mapp
 , 
2007
, vol. 
28
 (pg. 
228
-
237
)
Yuan
H
Liu
T
Szarkowski
R
Rios
C
Ashe
J
He
B
Negative covariation between task-related responses in alpha/beta-band activity and BOLD in human sensorimotor cortex: An EEG and fMRI study of motor imagery and movements
NeuroImage
 , 
2010
, vol. 
49
 (pg. 
2596
-
2606
)
Zaitsev
M
Hennig
J
Speck
O
Point spread function mapping with parallel imaging techniques and high acceleration factors: fast, robust, and flexible method for echo-planar imaging distortion correction
Magn Reson Med
 , 
2004
, vol. 
52
 (pg. 
1156
-
1166
)
Zhang
D
Raichle
ME
Disease and the brain's dark energy
Nat Rev Neurol
 , 
2010
, vol. 
6
 (pg. 
15
-
28
)
Zhang
Y
Brady
M
Smith
S
Segmentation of brain MR images through a hidden Markov random field model and the expectation maximization algorithm
IEEE Trans Med Imaging
 , 
2001
, vol. 
20
 (pg. 
45
-
57
)