Infraslow State Fluctuations Govern Spontaneous fMRI Network Dynamics.

Spontaneous brain activity as assessed with resting-state fMRI exhibits rich spatiotemporal structure. However, the principles by which brain-wide patterns of spontaneous fMRI activity reconfigure and interact with each other remain unclear. We used a framewise clustering approach to map spatiotemporal dynamics of spontaneous fMRI activity with voxel resolution in the resting mouse brain. We show that brain-wide patterns of fMRI co-activation can be reliably mapped at the group and subject level, defining a restricted set of recurring brain states characterized by rich network structure. Importantly, we document that the identified fMRI states exhibit contrasting patterns of functional activity and coupled infraslow network dynamics, with each network state occurring at specific phases of global fMRI signal fluctuations. Finally, we show that autism-associated genetic alterations entail the engagement of atypical functional states and altered infraslow network dynamics. Our results reveal a novel set of fundamental principles guiding the spatiotemporal organization of resting-state fMRI activity and its disruption in brain disorders.

Introduction

Spontaneous neural activity is ubiquitously present in the mammalian brain and persists across physiological states [1, 2]. In humans, whole-brain patterns of intrinsic brain activity are typically mapped by measuring spontaneous fMRI signal in the resting brain [3], an approach termed resting-state fMRI (rsfMRI). A large body of experimental work has shown that low-frequency fluctuations in the fMRI signal are temporally synchronous across multiple functional systems, delineating a set of reproducible topographies known as resting-state networks [3], which can be reliably identified also in lower species, such as primates [4] and rodents [5]. These observations have prompted a widespread use of interregional correlation between rsfMRI signals as an index of coupling or “functional connectivity” between regions. Importantly, further investigations have shown that rsfMRI activity is characterized by a rich temporal structure, involving dynamic reconfiguration into transient network states occurring on the timescale of seconds [6]. This research has promoted a view of the resting brain as an inherently dynamic system, in which highly evolving patterns of instantaneous activity interact over time in complex and transient communication patterns [7].

Recent animal studies have linked hemodynamic-based measures of intrinsic brain activity to low-frequency fluctuations in neural activity, as measured with calcium imaging [8, 9, 10]. Specifically, optical imaging studies in mice have implicated global waves of neural activity and transient neural co-activations as key contributors of hemodynamic-based measurements of functional connectivity [8, 11]. At the same time, it has also been shown that human rsfMRI network dynamics can be decomposed into recurring instances of regional peak fMRI activity [12], a feature that, in animals, appears to be related to transient variation in calcium co-activation patterns [11]. Taken together, these findings suggest that reconfiguration of spontaneous network activity may be driven by state transitions reflecting recurring patterns of slow-wave activity. This hypothesis is consistent with the initial recognition of putative temporal sequences of propagated fMRI activity [13, 14] and the observation of recurring metastable states in human rsfMRI time series [15, 16]. However, despite our progress in mapping and modeling whole-brain rsfMRI dynamics [17], the fundamental principles by which brain-wide patterns of spontaneous fMRI activity reconfigure and interact with each other remain undetermined.

To obtain a precise characterization of how spontaneous rsfMRI network reconfiguration occurs, we used a framewise clustering strategy to identify recurrent patterns of spontaneous fMRI activity in the mouse brain and analyzed how these spatial patterns are dynamically coupled. By using this approach, we mapped brain-wide patterns of spontaneous rsfMRI activity with voxel resolution in the resting mouse brain at the group and subject level, demonstrating high reproducibility across different rsfMRI datasets [5, 18, 19]. We report that whole-brain mouse rsfMRI dynamics are governed by a limited number of recurring whole-brain co-activation states and show that the dynamics of their transitions are characterized by infraslow phase coupling. Importantly, we further show that autism-related mutations alter the dynamics and topography of these fluctuating states, giving rise to a set of atypical network configurations that can be related to patterns of aberrant rsfMRI connectivity. Our findings suggest that resting-state brain activity in the mammalian brain is governed by recurring transitions between fluctuating brain-wide states and point at altered infraslow network dynamics as a key contributor to the aberrant rsfMRI network activity observed in brain connectopathies.

Results

Selective fMRI Frame Averaging Recapitulates Networks of Correlated Activity

It has been previously shown that selective averaging of fMRI frames exhibiting regional peaks of blood-oxygen-level-dependent (BOLD) activity spatially recapitulates rsfMRI networks obtained via seed-based correlation analysis [12, 20]. To introduce our analytical framework in the context of existing literature and as a first step toward the implementation of framewise rsfMRI clustering in the mouse, we probed whether a similar relationship holds true also in this species. To this aim, we used predefined anatomical regions as rsfMRI correlation seeds in our main dataset (500 time points) and spatially averaged individual fMRI frames as a function of their peak intensity signal to produce seed-based mean co-activation patterns (SB-CAPs). We next compared the obtained SB-CAPs with canonical rsfMRI correlation maps (Figure S1). Consistent with previous observations [12, 20], we found that previously described rsfMRI systems, including the hippocampal, latero-cortical, auditory-temporal, and default-mode networks (DMNs) [19], can be spatially reproduced by averaging a limited number (e.g., 15%; Figure S1A) of fMRI frames exhibiting peak BOLD activity in corresponding anatomical seed locations.

Whole-Brain fMRI Frame Clustering Reveals a Set of Recurring Functional Brain States

The observation that regional peaks of fMRI activity drive spatially structured network has been previously interpreted as evidence that a limited number of fMRI frames may provide key information about ongoing brain states [12, 20]. Departing from this interpretative framework, we carried out a k-means clustering analysis of all the individual rsfMRI frames, without any a priori temporal, anatomical, or intensity-based restrictions (Figure 1). In contrast to seed-guided selective fMRI frame averaging, this approach identifies, as prototypes of each of k spatial activity clusters, the different types of simultaneous single time frame co-activation patterns (shortened hereafter as CAPs or “functional states”) that most frequently (or most stereotypically) recur in the data [16]. This strategy provides a regionally unbiased characterization of recurring rsfMRI states independent of the local polarity and intensity of BOLD activity and can identify composite functional states characterized by patterns of simultaneous co-activation (above fMRI signal baseline) and co-deactivation (below fMRI signal baseline; Figure 1B). Importantly, because this approach classifies activity in all time frames without any a priori anatomical selection, it has the advantage of allowing the detection of manifold recurring brain-wide patterns of spontaneous brain activity and not only those that occur at either local peaks of BOLD activity in pre-determined regions of interest [16].

Figure 1

Identification of Recurring Brain States via Whole-Brain fMRI Frame Clustering

(A) Glass brain representation of the seed-based resting-state networks described in Figure S1 (DMN, default mode network; HCN, hippocampal network; LCN, latero-cortical network; PLN, postero-lateral network).

(B) Illustrative fMRI BOLD time course (SD units) in the anterior cingulate (ACg) (red), somatosensory cortex (SSs) (blue) as well as fMRI global signal (black) in a representative subject (brain slice illustrated in A). Note the presence of peaks of concordant or diverging BOLD activity in cingulate and somatosensory areas across time, suggestive of time-varying network reconfiguration; cingulate and somatosensory regions are concurrently co-activated and co-deactivated in F1 and F4, but they exhibit opposing BOLD activity in F2 and F3.

(C) These dynamic transitions can be captured and classified into recurring brain states by clustering fMRI frames into spatially congruent patterns (CAPs), using the k-means algorithm.

To identify a set of robust and prominent functional topographies representative of rsfMRI dynamics in the mouse brain, we first computed, for increasing cluster number k (2 < k < 40), how much variance is explained by the clustering algorithm (Figure S2A). The obtained explained variance curve, computed for the main dataset, revealed a clear elbow region encompassing the range k = 4 − 8, with the first six states alone explaining 59.6% of total variance. Notably, all the combined n = 34 additional clusters explained only 3.7% additional variance (totaling 63.3%, combined), with each additional cluster alone accounting for less than 1% of total explained variance (Figure S2B). These results suggest that a few recurring dominant network states explain the vast majority of rsfMRI temporal dynamics in the mouse brain.

We further refined our cluster number selection by selecting the value k as the highest value within the elbow region that ensured maximal between-dataset CAP reproducibility with respect to two additional independently collected rsfMRI datasets (n = 41 and n = 23, respectively; Figure S2C). With this procedure, we identified k = 6 states that are robustly conserved across the three rsfMRI datasets and that together account for the largest portion of mouse rsfMRI dynamics (Figures 2 and S2). Further credibility for the robustness of these six states as sets of genuinely co-active fMRI voxels is given by the identification of the same CAPs at the single subject level (described below) and the fact that these six patterns were conserved even when partitioning the main dataset into a higher number of states (k = 20; Figure S3; r > 0.81; all CAPs). It should be emphasized that the employed cluster selection strategy, although optimally suited to identify a set of states that dominate rsfMRI dynamics, does not rule out the occurrence and contribution of additional, albeit less prominent and less conserved across datasets, rsfMRI topographies in the mouse brain. In keeping with this notion, a high-order partitioning of our dataset into k = 20 CAPs revealed a set of additional anatomically plausible network configurations (k = 20; Figure S3A). The obtained CAPs, however, included also some spatially similar configurations (e.g., CAP 1, CAP 9, and CAP 20 and CAP 2, CAP 8, and CAP 10), indicative of possibly redundant intermediate state transitions. In light of these observations, we restricted our subsequent analyses to the first six clusters, as they account for the largest portion of mouse rsfMRI dynamics, they are reproducible across datasets, and this appears to be the finest partition exhibiting the highest cross-dataset reproducibility.

Figure 2

Recurring Functional States of the Mouse Brain

(A) Whole-brain representation of the functional brain states (CAPs) we identified at the group level. Red-yellow indicates co-activation (i.e., high fMRI BOLD signal), and blue indicates co-deactivation (i.e., low fMRI BOLD signal; p < 0.01; Bonferroni corrected).

(B) CAPs have been ordered based on their spatial properties by numbering consecutively states characterized by opposing BOLD co-activation patterns (i.e., 1-2, 3-4, and 5-6), as denoted by the negative correlations.

(C and D) CAP occurrence rate (C) and mean duration (D; mean ± SEM).

ACg, anterior cingulate cortex; AUD, auditory cortex; dHC, dorsal hippocampus; vHC, ventral hippocampus; HT, hypothalamus; ILA, infralimbic area; LAN, lateral amygdalar nucleus; MOp, primary motor cortex; Mos, secondary motor cortex; ORB, orbitofrontal cortex; PIR, piriform area; PL, pallidum; Rs, retrosplenial cortex; SSp, primary somatosensory cortex; SSs, secondary somatosensory cortex; ST, striatum; TeA, temporal association cortex; TH, thalamus; VIS, visual cortex. See also Videos S1, S2, S3, S4, S5, and S6 and Figures S2–S5.

fMRI States Encompass Known Connectivity Networks of the Mouse Brain and Can Be Identified at the Single Subject Level

An illustration of the six identified states is reported in Figure 2A. One defining characteristic of all the six CAPs is their configuration as a composite assembly of regional substrates encompassing previously characterized large-scale rsfMRI networks of the mouse brain. For example, CAP 1 shows a clear co-activation of primary and secondary motor-sensory areas belonging to the mouse latero-cortical network (LCN) together with deactivation of cortico-limbic regions and peri-hippocampal constituents of the mouse DMN. Similarly, CAP 5 encompasses co-deactivation of the DMN and co-activation of the hippocampal network. The presence of spatially prominent contributions of known rsfMRI connectivity networks in the identified states, plus their duration in the order of a few seconds (Figure 2D), implicates the observed patterns as putative time-varying constituents of spontaneous rsfMRI network dynamics as assessed with correlational techniques.

To investigate whether the selected six states are representative of spontaneous brain dynamics identifiable at the single subject level, we repeated our state detection using k = 6 on a subset of nineteen mice, for which we acquired rsfMRI images over a 30-min window. We next spatially matched each subject-level state with the corresponding CAP obtained at the group-level analyses (Figure 2A) to obtain voxelwise CAP incidence maps. These analyses revealed that all six states can be reliably identified at the single-subject level, with foci of very high cross-subject incidence in key network locations (Figure 3A). Importantly, the spatial distribution of the observed CAPs at the subject level (Figure 3B) closely recapitulates the features that we observed with group-level clustering (Figure 3B). These correspondences suggest that the six identified CAPs correspond to genuine brain states, representative of spontaneous functional reconfigurations occurring at the single-subject level.

Figure 3

Brain States Can Be Detected at the Single-Subject Level

(A) Incidence map of each CAP at the subject level (k = 6, 30-min rsfMRI acquisitions; p < 0.05; FDR corrected). Each voxel represents the proportion of subjects with significant co-activation (p < 0.05; FDR corrected) as its corresponding group-level CAP template.

(B) Spatial distribution of CAPs detected in a representative subject (p < 0.05; FDR corrected). Yellow indicates regional co-activation, blue regional co-deactivation.

See also Figure S4.

We next leveraged the identification of reliable states at the subject level to generate a standard co-activation map for each CAP across all the control subjects of the three datasets (n = 104; Figure S4) with the aim to identify core spatial features that could guide future CAP identifications. The obtained map highlighted opposing patterns of DMN and LCN co-activation and fronto-hippocampal reciprocity as high-confidence defining spatial features of CAP pairs 1-2 and 5-6, respectively. A widespread pattern of cortico-striatal co-deactivation was the most prominent spatial feature in CAPs 3 and 4.

Further analyses also ruled out a spurious contribution of signal-to-noise (SNR) fluctuations or motion to the emergence and spatial topography of CAPs (Figure S5A). Specifically, we found that framewise SNR was uniformly distributed across CAPs (one-way ANOVA; p = 0.98). Similarly, we found that strict censoring of putative “motion-affected” frames, leading to a rejection of 24% and 10% putative “high-motion” frames (framewise displacement thresholds; 75 or 100 μm), resulted in CAPs spatially undistinguishable from those observed by using uncensored frames (r > 0.98; Figure S5B). Moreover, putative “motion-affected” frames appeared to be uniformly distributed across CAPs (one-way ANOVA; p = 0.51 and 0.91, respectively; Figures S5C and S5D). These results argue against a significant contribution of motion artifacts or SNR fluctuations to our imaging results.

Functional States Exhibit Infraslow Network Dynamics

A notable feature of the identified CAPs is their configuration into state and anti-state pairs characterized by opposing patterns of functional co-activation (Figure 2B). This characteristic was especially prominent in CAPs 3 and 4 (r = −0.96) but also apparent in CAPs 1 and 2 (r = −0.80) and CAPs 5 and 6 (r = −0.66), the first pair being characterized by a clear anti-correlation between the DMN and LCN and the latter between hippocampal areas and the DMN (Figures 2A and 2B). The observation of opposite spatial configurations was not necessarily expected or resulting from the employed clustering procedure. This attribute suggests that the networks constituting these states continue to be similarly correlated (or anti-correlated) not only when these elements are co-activated but also when they are co-deactivated. This could happen if the state anti-state pairs undergo periodic state fluctuations, with opposite states reflecting peaks and troughs of this network dynamics. Based on the above results, we tested the hypothesis that the recurring transitions between CAPs could be described in terms of infraslow fluctuations.

To this aim, we computed at each instant of time the spatial correlation between each CAP and the BOLD fMRI signal in that time frame (Figure 4A). The obtained time series, hereafter referred to as “CAP time courses,” describe the spatial match between spontaneous brain activity and the specified CAP in the considered time frame.

Figure 4

fMRI States Exhibit Infraslow Network Dynamics

(A) Illustrative CAP and GS time course from a representative subject.

(B) Mean power spectral density of CAPs and the fMRI GS (mean ± SEM). Blue dashed vertical lines delimit the 0.01- to 0.03-Hz frequency band employed in subsequent phase analyses; horizontal red lines show the spectrum of randomly shuffled surrogate time courses.

When we computed the power spectra of each CAP time course (Figure 4B), we observed a peak of predominant power in the infraslow 0.01- to 0.03-Hz frequency band. Importantly, the power of the fluctuations within this band was significantly higher (p < 0.001; false discovery rate [FDR] corrected) than the corresponding spectra generated using sequence-locked surrogate time series of each CAP’s time course (Figure 4). Together, these results indicate that brain-wide spontaneous brain activity undergoes transitions between a set of fluctuating infraslow network configurations described by CAPs. Within this framework, the CAPs described here can be alternatively conceptualized as up and down states of three fluctuating network systems encompassing distributed patterns of spontaneous activity.

To visualize the spatiotemporal dynamics of CAP assembly and disassembly, we next computed the average of the whole-brain BOLD frames time locked around each CAP time course’s local maxima. These results, which are best visualized in the form of videos (Videos S1, S2, S3, S4, S5, and S6), clearly show how each state builds up and disassembles with spatial patterns that closely recapitulate wave-like propagating activity recently observed with calcium imaging in the mouse dorsal cortex [11].

Video S1. Gradual State Assembly and Disassembly of CAP1, Related to Figure 2

Time-locked averaged frames corresponding to the 2% highest peaks of CAP1 time course. Each selected frame was set as a time t = 0 reference event, and sampled frames within a −30 ≤ t ≤ 30 repetition interval.

Video S2. Gradual State Assembly and Disassembly of CAP2, Related to Figure 2

Time-locked averaged frames corresponding to the 2% highest peaks of CAP2 time course. Each selected frame was set as a time t = 0 reference event, and sampled frames within a −30 ≤ t ≤ 30 repetition interval.

Video S3. Gradual State Assembly and Disassembly of CAP3, Related to Figure 2

Time-locked averaged frames corresponding to the 2% highest peaks of CAP3 time course. Each selected frame was set as a time t = 0 reference event, and sampled frames within a −30 ≤ t ≤ 30 repetition interval.

Video S4. Gradual State Assembly and Disassembly of CAP4, Related to Figure 2

Time-locked averaged frames corresponding to the 2% highest peaks of CAP4 time course. Each selected frame was set as a time t = 0 reference event, and sampled frames within a −30 ≤ t ≤ 30 repetition interval.

Video S5. Gradual State Assembly and Disassembly of CAP5, Related to Figure 2

Time-locked averaged frames corresponding to the 2% highest peaks of CAP5 time course. Each selected frame was set as a time t = 0 reference event, and sampled frames within a −30 ≤ t ≤ 30 repetition interval.

Video S6. Gradual State Assembly and Disassembly of CAP6, Related to Figure 2

Time-locked averaged frames corresponding to the 2% highest peaks of CAP6 time course. Each selected frame was set as a time t = 0 reference event, and sampled frames within a −30 ≤ t ≤ 30 repetition interval.

Functional State Transitions Occur at Specific Phases of Global fMRI Signal

Interestingly, the power spectrum of the fMRI global signal (GS) also exhibits a power spectrum peak in the same infraslow band that characterizes CAP fluctuations (Figure 4B). This observation raises the question of how CAP infraslow dynamics may relate to that of the GS. We hypothesized that intrinsic fluctuations in GS may not reflect spatially undifferentiated ups and downs of whole-brain activity, but on the contrary, each phase of the GS may encompass the specific activation of selected subsets of possible brain states, each characterized by a characteristic profile of brain activity. To test this hypothesis, we measured the frequency of CAP occurrence at different phases of the fMRI GS fluctuations by filtering the GS signal in the infraslow band and computing the phase at each instant [21]. With the phase convention we used, phase values of 0 and correspond to peaks and troughs of the GS, respectively (Figure 5B). To understand whether different sections of GS fluctuation cycles correspond to different states, we next computed the circular distribution of GS phases at which each CAP occurred. Notably, we found that the occurrences of all six CAPs were not distributed uniformly across the GS infraslow cycle (Figures 5B and 5C) but were concentrated at specific ranges of the GS phase cycle, with all distributions exhibiting a significant deviation from circular uniformity (Raleigh test; p < 0.05; Bonferroni corrected). Specifically, CAP pairs 3 and 4 as well as 1 and 2 tended to occur around the trough and the peak of GS fluctuations, respectively (Figures 5B and 5C), albeit with different spread of phases, and occurrences of CAPs 5 and 6 were concentrated around time points in which the GS phase was and , respectively, corresponding to intermediate, off-peak GS levels. This latter observation lends strong support to our hypothesis that different phases of the GS fluctuations do not reflect just spatially unstructured ups and downs of global brain activity. The qualitative differences highlighted above were confirmed by quantitative tests, with evidence of significant differences between GS mean phases of all CAPs (Watson-Williams test; p < 0.05; Bonferroni corrected), with the only exception of CAPs 2 and 4. Additional testing using temporally permuted GS filtered time courses showed that all original GS phase distributions exhibit significant deviation from circular uniformity with significantly high concentration parameters (p < 0.001; permutation test).

Figure 5

fMRI States Occur at Specific Phases of Global Signal Fluctuations

(A) Instantaneous phase of the GS in a representative subject. Colored dots mark the occurrence of each CAP.

(B) CAPs occur at specific phases of global signal fluctuations. Phase dispersion is plotted as circular variance around the circular mean (top horizontal bars). CAP encoding is described according to color scheme in (A).

(C) Circular distribution of GS phases at each CAP’s occurrence within a GS cycle. For each distribution, the resulting vector (magnitude and phase) is shown as a black radial line.

See also Figures S3 and S6.

Importantly, these dynamic features appear to be independent of the number of selected clusters, because a high order portioning of our dataset (k = 20; Figure S3B) did not affect CAP 1-6 phase relationship. Moreover, a computation of circular-linear correlation between spatial similarity of each CAP pair, and their GS phase difference, did not reveal any statistically significant association (p > 0.2; all CAP pairs; FDR corrected), ruling out an artifactual collation at similar phases of spatially related CAPs, because of possible spatial collinear effects.

Functional States Are Phase Coupled

Although the above analysis suggests that CAPs act as infraslow fluctuating networks, it does not fully clarify whether these CAPs fluctuate independently or whether the phases at which they fluctuate within a single cycle of the global fMRI signal are coupled to each other. If the latter is the case, phase relationships between CAPs must be observed not only when pooling all phase occurrences together (as in Figure 5C) but also when considering phase relationships between occurrences of different CAPs within the same (or immediately adjacent) GS cycle. We thus computed the GS phase angular difference between CAP occurrences in the same GS cycle or across cycles that are immediately adjacent in time (i.e., next or previous GS cycle; Figure S6). We found that GS phase difference between occurrences of the same CAP (diagonal panels, Figure S6) were concentrated around zero, suggesting that a given state appears in general for a short range of adjacent phases during a given cycle. Moreover, reciprocal CAPs appeared with a phase difference of in the same cycle, suggesting that each GS cycle reflects at least in part the alternation between peaks and troughs of a specific spatially structured network. Finally, we found significantly locked distribution of phase differences between non-reciprocal CAP pairs (e.g., between CAPs 4 and 5 or CAPs 3 and 6). Interestingly, when non-reciprocal CAPs co-occur, they do so mostly within the same half of the cycle (i.e., all or most instances occur with a phase difference of less than ). Thus, these CAPs can co-occur within the same cycle with a relatively broad phase difference. The fraction of cycles in which non-reciprocal CAPs co-occur was high, between 75% and 87%, depending on the non-reciprocal CAP pair considered. Collectively, these results show that the identified functional states are phase coupled and suggest that GS fluctuations, rather than being the result of unstructured changes of activity, reflect coupled dynamics of interacting structured networks along the entire fluctuation cycle.

Altered Patterns of rsfMRI Connectivity Entail Non-canonical Functional State Dynamics

The relationship between patterns of spontaneous fMRI activity and correlational rsfMRI network topographies (Figure S1) suggests that aberrant spatiotemporal structure of instantaneous fMRI states may underlie alterations in steady-state rsfMRI activity and “functional connectivity” [22, 23]. To probe this hypothesis, we mapped functional states in mice haploinsufficient for the chromatin-remodeling gene Chd8. This mutation recapitulates a major genetic risk factor for autism spectrum disorders (ASDs) and produces rsfMRI over-connectivity in hippocampal and motor cortical areas [24]. As described above, brain state identification using k = 6 in control animals (Chd8+/+) recapitulated the CAPs observed in our two larger control datasets (Figure S2C). Notably, CAPs in Chd8 mutants (Chd8+/−) showed robust region-specific alterations in fMRI signal intensity, entailing an atypical engagement of specific anatomical substrates with respect to the corresponding states in control littermates (Figure 6B; T(40) > 2.8; cluster corrected; p < 0.05). Specifically, CAPs 1 and 2 in Chd8+/− mice exhibited an aberrant reciprocal involvement of thalamic and dorsal hippocampal regions. Similarly, cingulate and mid-thalamic recruitment appeared to be defective in CAPs 3 and 4 of Chd8+/− mice. Foci of aberrant motor sensory and prefrontal co-activation were found in CAPs 5 and 6. These results provide evidence of non-canonical brain states in mice harboring a key ASD-risk mutation, entailing robust alterations in regional patterns of spontaneous BOLD fMRI activity. We did not observe any significant between-group difference in CAP occurrence rates (p > 0.16) or average CAP duration (p > 0.7).

Figure 6

Altered Brain States in a Genetic Mouse Model of Autism

(A) Representative coronal and sagittal cuts showing the functional organization of CAPs 1-6 in control (Chd8+/+, left) and mutant mice (Chd8+/−, right; p < 0.001; familywise error [FWE] cluster corrected).

(B) Anatomical location of the regions exhibiting significant between-group differences in BOLD co-activation intensity for each of the six identified CAPs (p < 0.001; FWE cluster corrected). Red or blue indicate areas exhibiting regionally increased or decreased BOLD co-activation, respectively, in mutant mice with respect to control littermates.

(C) Delayed CAP occurrence within GS cycles with respect to control mice in Chd8+/− mutants (p < 0.01; William-Watson test for circular mean homogeneity; Bonferroni corrected; all CAPs, except for CAP3).

(D) Seed-based rsfMRI correlation differences in Chd8 mutants. Inter-group differences (p < 0.05; FWE cluster corrected) are depicted with respect to a seed-pair in the SS to replicate the rsfMRI findings reported in [24]. The location of hippocampal (HCP) seed used for correlational profiling is also illustrated (HCP1, 2, and 3).

(E) Seed-based correlation profiling upon regression of individual CAP time courses or the GS in both groups (∗∗∗p < 0.001; two-way repeated-measures ANOVA; mean ± SEM).

Cg, cingulate cortex; dHCP, dorsal hippocampus; HCP, hippocampus; vHCP, ventral hippocampus; SS, somatosensory cortex.

We next investigated whether these alterations also affected CAP infraslow network dynamics (Figure 6C) by sampling the GS phase of the occurrence of each state in control and Chd8 mutants. We found CAPs in Chd8+/− mutants to occur at the same phase cycle and in the same order as the corresponding state of control mice, corroborating between-group CAP matching. However, all the functional states of Chd8+/− mice exhibited a significantly delayed occurrence (−30.9° < Δφ < −46.5°) with respect to GS fluctuation phases (Figure 6C; p < 0.01; William-Watson test), with the exception of CAP 3, in which no genotype-dependent phase delay was observed. The distribution of phases happening within the same or adjacent GS cycles was found to significantly deviate from circular uniformity in both mutant and control mice (Raleigh test; p < 0.05; Bonferroni corrected), showing that CAPs are phase coupled in both cohorts. However, there was a significant genotype-dependent difference in the circular mean of the phase distributions of functional states occurring within same or adjacent cycles, with the only exception of phase relationships between CAPs 1-4 and 2-4 (Williams-Watson test; p < 0.05; Bonferroni corrected), with major GS phase differences occurring between CAP 3 and CAPs 4-6. Duration and occurrence of CAPs were instead similar across cohorts (p > 0.16; all CAPs; FDR corrected). Collectively, these findings suggest that ASD risk mutations can alter both the regional topography and dynamics of recurring rsfMRI states.

In an attempt to relate the observed dysfunctional state dynamics to aberrant rsfMRI “connectivity,” we anatomically profiled interregional rsfMRI correlation between somatosensory and hippocampal areas previously shown to be hyperconnected in Chd8 mutants [24], before and after regressing each CAP’s time course independently (Figure 6E). Independent regression of the time course of either CAP 1 or 2 removed inter-group differences in cortico-hippocampal rsfMRI connectivity observed with seed-based correlation analyses (p > 0.6; two-way ANOVA; Figure 6E). Regression of all the other CAPs or the GS did not substantially affect somatosensory-hippocampal rsfMRI coupling (p < 0.0005; two-way ANOVA; Figure 6E). These results suggest that altered interregional synchronization in rsfMRI time series may involve aberrant co-activation topographies.

Discussion

Dynamic mapping of rsfMRI datasets has revealed that spontaneous brain activity is a non-stationary phenomenon, involving reconfiguration into recurring, possibly propagating [13, 14, 15, 16], functional states [17, 25]. Here, we expand this view by reporting that time-varying patterns of spontaneous fMRI activity can be described in terms of simple dynamical rules. Specifically, we show that the spatiotemporal dynamics of spontaneous rsfMRI activity is governed by a limited number of infraslow network transitions occurring at different phases of global fMRI signal fluctuations. Importantly, we found the dynamics and spatial topography of these states to be altered in a genetic model of developmental disorders.

These findings advance our understanding of the principles underlying spontaneous fMRI network dynamics in three directions. First, our results point at multi-scale infraslow fluctuations as a superstructure of states that critically regulate integration and decoupling of concurrently active neuronal communities, as previously postulated [26]. Within this framework, our results provide key, yet so far elusive, experimental evidence to the influential theory that spontaneous brain-wide network activity undergoes infraslow assembly and disassembly [27, 28, 29]. Second, our work also provides a novel framework for the description and interpretation of dysfunctional rsfMRI “connectivity” in autism [30], an endophenotype replicated in mouse lines harboring autism risk gene mutations [22, 31, 32, 33, 34]. The observation that ASD-relevant rsfMRI connectivity aberrancies are associated with topographically altered fMRI states challenges the predominant view of aberrant rsfMRI connectivity as a simple proxy for circuit-specific decoupling or miswiring and suggests an interpretational model in which stationary rsfMRI patterns of disconnectivity represent only a partial representation of a possibly larger and more widespread dysfunction in regional coupling and neuronal dynamics. Given the established contribution of synchronized rhythms in mediating long-range neuronal communication [35], alterations in brain-wide infraslow network dynamics may be indicative of abnormal mechanisms of integration and decoupling of distributed neural inputs, a hypothesis consistent with the prominent sensory processing and cognitive impairments observed in ASD patients. Finally, departing from correlation-based methods to describe rsfMRI dynamics, the use of fMRI signal amplitude as the basis for inter-group differences in state topographies provides an easily interpretable readout that can be used to translate patterns of disconnectivity into physiologically relevant patterns of spontaneous neural activity.

The identified brain states exhibit a composite network structure, providing novel insights into the functional and evolutionary structure of rsfMRI network dynamics in the mammalian brain. The opposing DMN and LCN co-activations in CAPs 1-2 and 5-6 recapitulate a cardinal feature of human DMN organization [36] and corroborate the presence of a tight inverse coupling between these two systems [5, 37]. A similarly competing relationship between DMN regions and latero-cortical constituents of the task-positive network (TPN) has been identified in humans using autoregressive pattern-finding algorithms [13]. An evolutionary back projection of these network features would implicate the LCN as a putative precursor of the TPN in rodents, as previously postulated [18]. Notably, the same CAPs also provide a spatial delineation of the rodent DMN and LCN in non-correlative terms, supporting initial descriptions of this network to comprise thalamo-frontal and peri-hippocampal cortical components in rodents [5, 18]. Finally, the contrasting polymodal and unimodal organization of the DMN and LCN areas, together with the segregable co-activation of visual and motor-sensory cortical areas along a principal cortical axis (CAPs 3 and 4), recapitulate the hierarchical topography of cortical gradients recently mapped in humans with rsfMRI and structural imaging [38].

Recent investigations have linked fluctuations of the fMRI GS to vigilance [39], glucose metabolism [40], and arousal mediated by ascending nuclei [41, 42]. Our results are consistent with the hypothesis that the fMRI GS may encode for key neuronal-relevant information and support a view in which each GS cycle is the sum of different, partially overlapping network configurations. Although the intrinsic drivers of these reconfigurations remain elusive, the observation of foci of co-deactivation in basal forebrain areas in CAP 4 (Figure 2A) recapitulates similarly contrasting patterns of global activity in cortical and basal forebrain regions observed in humans [41] and is consistent with a role of ascending modulatory activity in driving these infraslow fluctuations [42].

Recent optical imaging studies in mice provide key mechanistic clues as to the neural determinants of the state fluctuations we describe in this work. Calcium imaging in mice has revealed that spontaneous neural activity entails infraslow neural waves as well as transient co-activations spanning the entire dorsal cortex [2, 8]. Concurrent hemodynamic and neural measurements have convincingly linked the two phenomena [8, 10]. Recently, this relationship has been expanded to relate cortical co-activation patterns of calcium activity with spatially structured hemodynamic fluctuations [11]. Importantly, the spatial extension of these fluctuations appears to recapitulate the contrasting involvement of DMN and LCN areas, similar to what we describe here. These spatial correspondences argue for a neural origin of the identified fMRI states and support a view in which brain-wide fMRI fluctuations are guided by intrinsic infraslow cycling between slowly propagating neural activity. In keeping with this framework, the aberrant state topographies observed in the autism model could be a proxy for chronic alterations in mesoscale neuronal co-activity (i.e., neuronal ensembles), a dysfunction recently described in other mouse models of developmental disorders [43].

Our results suggest that a few recurring network states explain the most rsfMRI temporal dynamics. It is interesting to compare this number with the considerably higher proportion (ca. 15) of spatially independent sources identified using ICA in mouse rsfMRI datasets [18]. In keeping with our findings, multiple spatially independent components have been typically identified in humans, (e.g., 15–40), but the dynamic states arising from these spatial sources is typically very small (e.g., 4–7) [44]. Our observations suggest that, under resting conditions, most of these spatially independent sources may be concomitantly engaged in stereotyped and coordinated patterns of activity. Further work entailing perturbational manipulations of some of these sources [45] may help clarify whether and how different spatial sources combine to produce the rich, yet highly stereotyped, structure of resting state dynamics.

A methodological aspect that deserves further consideration is how to best compare CAP topographies across different datasets or experimental conditions. Leveraging the reproducibility of our CAPs and the presence of a genetically identified reference population in our mouse studies, here, we compared CAPs identified separately in Chd8 control and mutant littermates. However, alternative methods, including the identification of CAPs from combined datasets, appear to be equally applicable to the detection of between-group functional differences. Understanding how to best perform such comparisons and design the analytical approach accordingly is an important goal of future research.

In conclusion, we document that fMRI network activity is governed by infraslow transitions between fluctuating brain-wide states and show that the spatiotemporal structure of these reconfigurations can be described in terms of simple dynamical rules. Our findings reveal a novel fundamental principle guiding the spatiotemporal organization of resting-state fMRI activity and point at altered infraslow network dynamics as a key contributor to the aberrant rsfMRI network activity observed in brain connectopathies.

STAR★Methods

Key Resources Table

Lead Contact and Materials Availability

Further information and requests for resources and reagents should be directed to and will be fulfilled by the Lead Contact, Alessandro Gozzi (alessandro.gozzi@iit.it).

Experimental Model and Subject Details

All in vivo experiments were conducted in accordance with the Italian law (DL 26/214, EU 63/2010, Ministero della Sanità, Roma) and the recommendations in the Guide for the Care and Use of Laboratory Animals of the NIH. Animal research protocols were reviewed and consented by the animal care committee of the Italian Institute of Technology, and Italian Ministry of Health. Three independent datasets were used in this study: datasets 1 and 2 consisted of cohorts of n = 40 and n = 41 (12-18 week-old) male C57BI6 mice, respectively, purchased from the Jackson Laboratory (Bar Harbor, USA). A third dataset consisted of n = 23, 15-18 week-old Chd8+/+ (9 males, 14 females), and n = 19 Chd8+/− (6 males, 13 females) mice. These animals have been generated on a C57Bl6/J background by Albert Basson, as described in [24]. The animals were housed in a 12:12 hr light-dark cycle in individually ventilated cages, with access to food and water ad libitum.

Method Details

Resting state fMRI data acquisition

Experiments were performed on a set of n = 40 adult (12-18 week-old) male C57Bl6/J mice. This set of animals was imaged at the IIT laboratory in Rovereto (Italy) in 2016, and is the dataset that we have used throughout our article to describe rsfMRI dynamics, unless otherwise stated. The animal preparation protocol for experimental measurements has been described in great detail [30, 44]. Briefly, mice were anaesthetized with isoflurane (5% induction), intubated and artificially ventilated (2%, surgery). The left femoral artery was cannulated for continuous blood pressure monitoring and terminal arterial blood sampling. At the end of surgery, isoflurane was discontinued and substituted with halothane (0.75%). Functional data acquisition commenced 45 minutes after isoflurane cessation. Mean arterial blood pressure was recorded throughout imaging sessions. Arterial blood gases (paCO2 and paO2) were measured at the end of the functional time series to exclude non-physiological conditions.

All rsfMRI data were acquired with a 7.0 Tesla MRI scanner (Bruker Biospin, Ettlingen) equipped with BGA-9 gradient set (380 mT/m, max. linear slew rate 3,420 T/m/s), using a 72 mm birdcage transmit coil, and a four-channel solenoid coil for signal reception as previously described [19]. For each session, high-resolution anatomical images were acquired with a fast spin echo sequence (repetition time (TR)/echo time (TE) 1200/15 ms, matrix 192 × 192, field of view 2 × 2 cm2, 18 coronal slices, slice thickness 0.60 mm). Co-centered single-shot blood-oxygen level dependent (BOLD) EPI time series were acquired using an echo planar imaging sequence with the following parameters: TR/TE 1200/15 ms, flip angle 30°, matrix 100 × 100, field of view 2 × 2 cm2, 18 coronal slices, slice thickness 0.50 mm, 500 (n = 21) or 1500 (n = 19) volumes and a total rsfMRI acquisition time of 10 or 30 minutes, respectively. All the group analyses were carried out on the first 500 time points (10 minutes). The single subject CAP analysis was limited to the n = 19 subjects in which we acquired 1500 time points.

To corroborate the reproducibility of our findings across independent datasets, we applied the whole analytical pipeline to an additional dataset composed of n = 41 male C57Bl6/J mice in which we acquired rsfMRI timeseries (n = 300, 6 minutes) using the same sedation protocol and image parameters employed in the present study. This set of animals was scanned at the IIT laboratory in Pisa (Italy) in 2014. A characterization of functional network organization in this set of animals has been previously described [19]. We refer to this animal cohort as to dataset 2. Finally, to assess the ability of our analytical pipeline to detect aberrant states in mouse models of brain pathology, we applied our analytical framework to a cohort of Chd8 haploinsufficient mice, a relevant subtype of autism spectrum disorder, which has been previously described to present aberrant rsfMRI network activity [24]. Briefly, rsfMRI imaging was performed on 15-18 week old mice (n = 23 Chd8+/+, 9 males n = 19 Chd8+/−, 6 males both lines have a C57Bl6/J background), each with 500 volumes (10 minutes), using the same animal preparation protocol and rsfMRI acquisition parameters of dataset 1. The control group in this study was employed as a third independent dataset for the validation of our clustering procedure (dataset 3). We merged data from both males and females from these cohorts as we previously reported no major sex-related differences in behavior and connectivity in these mice [24]. This set of animals was scanned at the IIT laboratory In Rovereto, in 2017.

Data preprocessing

Data preprocessing was carried out as recently described [29]. Briefly, fMRI time series were despiked, motion corrected, and spatially normalized to an in-house mouse brain template [18] yielding a final normalized spatial resolution of 0.1 × 0.1 × 0.5 mm3 (192 × 192 × 24 matrix). Head motion traces and the mean ventricular signal (average fMRI time series within a manually-drawn ventricle mask from the template) were regressed out. The resulting images were band-pass filtered using a 0.01 – 0.1 Hz band, spatially smoothed using a Gaussian kernel of 0.5 mm FWHM, and z-scored voxel-wise.

Quantification and Statistical Analysis

Seed based CAPs and correlation analysis

To probe the relationship between networks inferred from conventional rsfMRI seed-based correlation analyses, and those described by high regional fMRI activity at only a few critical time points [12, 20], we first extracted individual rsfMRI volumes (here referred to as “frames”) at time points with suprathreshold signal in a set of a priori regions of interest. We next spatially averaged the extracted rsfMRI volumes to provide CAPs, which we refer to here as seed-based mean CAPs.

Seed location was chosen based on prior rsfMRI mapping in the mouse [18, 19]. For each of the probed regions, we also computed a canonical group-level correlation map using the corresponding regional rsfMRI signal as seed. We next computed the spatial correlation between the seed-based CAPs and their corresponding correlation maps by retaining rsfMRI frames exceeding a predefined intensity threshold, covering the whole 0-99th percentile range as previously described [12]. For illustrative purposes, we generated representative seed-based CAPs by averaging all the fMRI frames with the highest 15% BOLD signal intensity across all subjects (68 ± 5 out of 500 frames, mean ± SD, all subjects), and using a T threshold of 7, corresponding to p < 0.01, Bonferroni corrected [12, 45]. For each region we also generated a group-level canonical seed-based correlation map, which we thresholded at T = 7, corresponding to p < 0.01, Bonferroni corrected).

Whole-brain CAPs analysis

We used spatial clustering of individual fMRI frames to identify whole-brain patterns of simultaneous co-activation of brain activity with voxel-resolution. Following prior studies [12, 15, 16], we carried out frame-wise clustering of brain-wide mouse rsfMRI images by retaining the voxels that were in the top 10% or in the bottom 5% of all BOLD signal in a single time frame. Pilot studies showed that this step expedites the clustering procedure without qualitatively affecting the outcome of CAP clustering. Specifically, we found that frame-wise clustering without the use of intensity thresholding led to a marginal change (< 1.64%) in the time-frames assigned to a specific CAP, and that the resulting CAPs are indistinguishable from those obtained with intensity thresholding (Pearson’s correlation, > 0.98). All the subsequent frame-averaging procedures, leading to the definition of mean CAP maps, were carried out by averaging un-thresholded fMRI frames. Similarly, all our dynamic analyses were carried out using native (intensity un-thresholded) CAPs, and/or rsfMRI timeseries.

The preprocessed fMRI frames were formatted into N-dimensional vectors , with being the number of frames, and the amount of voxels. Such frames were clustered using the k-means++ algorithm [46], which partitions the vector-set into k clusters , such that the sum of within-cluster distances D in the following equation is minimized:being the mean of the fMRI frames in each cluster: , and the distance between the fMRI frame at time j and the cluster mean, measured as one minus the spatial Pearson’s correlation coefficient. The k-means++ algorithm provides an optimized initialization of each centroid (seeding), such that distant centroids have a higher chance of being chosen as starting seeds for the clustering procedure. The addition of this step has been shown to outperform conventional k-means with random seeding both in terms of clustering performance and speed [46].

k-means clustering was carried out on the whole rsfMRI dataset 1 (40 × 500 frames), with k ranging from k = 2 up to 40, using Pearson correlation between fMRI data in different time frames as clustering distance measure, and using 15 replicates with 500 iterations each. Departing from previous investigations [12, 15, 16], we carried out k number selection generalizing a more quantitative procedure which was previously employed to map robust, recurrent states in electrophysiological recordings [47]. To this aim, we computed, for increasing k, how much variance is explained by the clustering algorithm, defined as the ratio between the between-cluster variance and the total variance (within-cluster + between-cluster variance). Within-cluster variance was computed as the averaged (over clusters) sum of square distances between elements in a cluster and its centroid. Between-cluster variance was computed as the averaged square distance between a cluster centroid and the centroid of all clusters or centroid of all data [47].

To assess reproducibility of the corresponding CAPs across the three independently collected datasets used in this study we computed CAPs progressively increasing k in each control dataset (n = 40, n = 41 and n = 23, respectively). Between-dataset CAP matching was carried out using the Hungarian algorithm [48, 49], followed by a canonical measurements of spatial correlation (Pearson’s r), assuming anatomical correspondence with a CAP correlation > 0.45, corresponding to a p < 10e-5, permutation test. This procedure identified k = 6 CAPs conserved across all the three rsfMRI datasets.

To compute the centroid of each cluster (which we took as a CAP), the fMRI frames assigned to each cluster were averaged voxel-wise, and normalized to T-scores (p < 0.01, Bonferroni corrected), permitting to visualize the mean voxel-wise distribution of fMRI BOLD signal for each of the identified CAP. For visualization purposes, the obtained maps were thresholded to T-scores > 7, corresponding to a p < 0.01, Bonferroni corrected. For each CAP we next computed its occurrence rate (i.e., the proportion of frames assigned to each CAP) and duration (i.e., the average number of consecutive frames belonging to the same CAP at each occurrence) for each of the n = 40 subjects. We also computed between-CAP spatial similarity, defined as pairwise Pearson’s spatial correlation between all the identified CAPs.

To assess a possible confounding effect of fluctuation in signal-to-noise ratio to the detection of CAPs, we computed SNR as the ratio between the mean BOLD signal in a somatosensory cortical area (e.g., a region highly affected by GS fluctuations) and the standard deviation of the background signal outside the brain for each subject’s timeseries, frame by frame. We next computed for each fMRI frame the distribution of SNR at the occurrence of each CAP.

To assess the ability of our approach to detect CAPs also at the subject level, we repeated the clustering analysis in a subset of 19 subjects of the main dataset for which we acquired extended (30 minutes, 1500 frames) time series, using the group-level CAPs as initial centroids. The CAPs found in each subject were next matched to the group-level templates using the Hungarian Algorithm as previously described [48, 49], and their spatial correlation was computed using Pearson’s spatial correlation between maps (r). To illustrate CAP incidence across individual subjects, we mapped, for each voxel in a CAP, the proportion of subjects which had a significant co-activation (t test, p < 0.05, FDR corrected) and polarity (i.e., positive or negative BOLD signal) of the corresponding group-level CAP template.

Pilot analyses in which we undersampled each subject’s timeseries from 500 down to 72 frames, showed that the resulting subject-level CAPs exhibit a stable spatial correlation with the corresponding group-level CAPs up to a minimal length of approximately 125 time-points (not shown).

Inverse spatial configuration between CAP anti-CAP pairs at the subject level was quantified by building CAP maps for each subject (n = 40, dataset #1), and cross-comparing their spatial similarity. We found CAPs 1 and 2 to exhibit a mean anticorrelation of - 0.60 C.I. [-0.63, −0.57], p = 2-31, CAPs 2 and 4, −0.68, C.I.[-0.72, −0.7 = 65], p = 1.4e-32, and CAPs 5 and 6 −0.44, C.I.[-0.50, −0.38], p = 3.4e-17.

To obtain a standard co-activation models of all CAPs, we first computed subject-level CAPs (k = 6) for all the n = 104 control subjects of datasets #1, #2 and # 3. We next generated a cross-subject incidence map based on voxel polarity observed in single-subject CAP maps and computed for each CAP the corresponding p values from a one-tailed directional binomial test against chance (null, q = 0.5).

CAP dynamics

To investigate the dynamics of CAP evolution, we generated a CAP-to-frame correlation time course (which we term “CAP time course”). We computed these spatial correlations in each time frame as the Pearson correlation between the centroid of the considered cluster and the fMRI activity in the considered time frame of a given subject as previously described [42]. To describe the assembly and disassembly dynamics of CAPs, we used the method devised by Liang et al. [42]. Briefly, for each CAP we selected the fMRI frames corresponding to the local maxima within the CAP’s time course, limiting the selection of peak events to 2% of all the concatenated frames, and adjusting this number for each CAP’s occurrence rate. For example, for CAP 1 (occurrence rate = 0.18) we sampled 0.18 × 400 frames corresponding to the CAP’s time course local maxima. Each selected frame was set as a time t = 0 reference event, and we next sampled frames within a −30 ≤ t ≤ 30 repetition interval. All events were time-lock averaged, leading to a dynamic portrayal of mean temporal evolution of CAP assembly and disassembly in the form of a concatenated frames (Videos S1, S2, S3, S4, S5, and S6).

The temporal structure of the CAP time courses for each subject was also assessed by computing its power spectrum. Pilot analyses showed that the power spectra of CAPs 1-6 are not affected by the number of clusters k, a finding consistent with the power spectral profile of the fMRI GS being very similar to the one of CAPs. To assess the significance of the observed infra-slow fluctuations (0.01 - 0.03 Hz), we calculated the power spectrum of surrogate rsfMRI timeseries obtained by jointly shuffling all the six CAP time courses using the same random sequence of time points. This strategy preserves the signal amplitude distribution as well as the correlation structure of the original. Statistical assessments were next carried out upon 1000 iterations, using a one-tailed two-sample t test for each frequency in the band (FDR corrected). Further analyses also demonstrated that the spectral density of CAPs obtained with k = 20 are characterized by peaks of infra-slow activity closely resembling those observed with k = 6 (data not shown).

To compute the instantaneous phase of the GS in the infraslow range, we first band-pass filtered the GS time-courses between 0.01-0.03 Hz. We then used the Hilbert Transform [21] to decompose the GS signal into an analytical signal with a characteristic instantaneous phase and amplitude. We divided each subject’s instantaneous GS phase signal into cycles in the range [0, 2π], and within each cycle, collected the GS phase values at each CAP’s occurrence. Using the same filter design described above, we filtered and again normalized the CAP time courses, and sampled the GS phase at each CAP occurrence only when the CAP time course at that instant was above 1 SD, in order to ensure that a specific frame pertained to a specific CAP. Using the MATLAB CircStats toolbox [50], we computed circular statistics of the obtained distribution of GS phases at each CAP, and represented their dispersion in a cosine cycle representing a GS fluctuation, using the circular variance to quantify dispersion around the circular mean. We tested the significance of the obtained distributions by repeating the above-mentioned procedure using randomly permuted GS time courses (1000 iterations) and comparing, for each CAP, the concentration parameters of the original GS-phase distributions with the ones obtained with permuted GS time courses. To probe the presence of phase-coupling between CAP occurrence and GS infraslow dynamics, we computed the angular differences (phase differences) of the GS between occurrences of a given CAP within a GS cycle, and occurrences of another CAP within the previous, current, and subsequent cycle. Again, GS phase samples at each CAP were only considered if their filtered values at that instant were above 1 SD at the corresponding instances.

To rule out a contribution of collinearity to phase analyses related to the use of our projection method, we computed circular-linear correlation [50] between spatial similarity (r) of each pairs of CAPs, and the GS phase difference between their occurrence. This analysis did not reveal any significant association for any of the CAP pairs examined (p > 0.2, all CAP pairs, FDR corrected), suggesting that the projection method in our case does not artificially collate together clusters with similar spatial profiles at similar phases.

CAP dynamics in a genetic model of autism

The same analytical pipeline described above was applied to rsfMRI time series recorded in n = 23 Chd8+/+ and n = 19 Chd8+/− littermates [24] after standard preprocessing step. CAPs were computed independently in each group using k = 6, as this value ensures the highest cross-dataset reproducibility of the identified clusters (Figure S1). Between-group CAP matching was obtained via the intersection of three convergent criteria. We first analyzed spatial correlation between the CAPs of each group, to independently identify three CAP and anti-CAP pairs. We next computed the spatial similarity of the identified CAP pairs across the two groups, resulting in a spatial correlation between CAP pairs 1-2, 3-4 and 5-6 ranging from 0.45 to 0.72. A visual inspection of the obtained correlation matrix confirmed that any other reordering of the Chd8+/−CAPs would lead to a much lower correlation match between genotypes. Finally, the obtained matching was corroborated by the observation that the phase occurrences of the matched CAPs were consistent across genotypes (Figure 6C).

Inter-group differences in CAP anatomy were mapped using a two-sample t test and family-wise error (FWE) cluster correction (p < 0.05, cluster-defining threshold T(40) = 2.8). We then computed CAP features and CAP time-courses for all subjects in each strain, and sampled the filtered GS-phase distributions at each CAP occurrence for each group of mice independently. Differences between preferred GS-phase at CAP occurrences were computed using a William-Watson test for circular mean homogeneity, Bonferroni corrected for six comparisons. We also sampled the angular differences (phase differences) of the GS between occurrences of a given CAP within a GS cycle, and occurrences of another CAP within the previous, current, and subsequent cycle. Again, GS phase samples at each CAP were only considered if their filtered values at that instant were above 1 SD at the corresponding instances.

Finally, to assess the involvement of individual CAPs in the seed-based rsfMRI correlation differences previously described between hippocampal and motor-sensory areas [24], we computed differences in whole-brain correlation maps between a 3 × 3 X 1 voxel seed placed bilaterally in the somato-sensory cortex. Inter-group differences were assessed using a two-sample t test and family-wise error (FWE) cluster correction (p < 0.05, cluster-defining threshold T(40) = 2.8). We then recomputed Pearson correlation between this seed and selected foci of over-connectivity in the hippocampus. The rsfMRI correlation profiles were computed before and after regressing each CAP’s time course independently in each subject. Differences in rsfMRI correlations were assessed by means of a repeated-measures two-way ANOVA.

Data and Code Availability

rsfMRI datasets #1 and #2 are available for download at the following links (https://doi.org/10.17632/7y6xr753g4.1, https://doi.org/10.17632/r2w865c959.1 and https://doi.org/10.17632/thpszcwcgx.2). rsfMRI dataset #3 and the code employed to carry out our analyses can be obtained upon reasonable request to the lead contact, Alessandro Gozzi.

REAGENT or RESOURCE	SOURCE	IDENTIFIER
Deposited Data
rsfMRI dataset #1, part 1 (20 mice)	This paper	Mendeley Data: https://doi.org/10.17632/7y6xr753g4.1
rsfMRI dataset #1, part 2 (20 mice)	This paper	Mendeley Data: https://doi.org/10.17632/r2w865c959.1
rsfMRI dataset #2 (41 mice)	This paper	Mendeley Data: https://doi.org/10.17632/thpszcwcgx.2
Experimental Models: Organisms/Strains
Mouse: C57BL/6J	Charles River	C57BL/6J; RRID: IMSR_JAX:000664
Mouse: Chd8+/− transgenic mouse line	Albert Basson. See [24]	https://doi.org/10.1093/cercor/bhy058
Software and Algorithms
MATLAB 2017b	MathWorks	https://www.mathworks.com; RRID: SCR_001622
Image clustering	MATLAB Statistics and Machine Learning Toolbox	https://www.mathworks.com/help/stats/index.html; RRID: SCR_001622
Signal processing	MATLAB Signal Processing Toolbox	https://www.mathworks.com/help/signal/index.html; RRID: SCR_001622
Circular Statistics	MATLAB Circular Statistics Toolbox	https://www.mathworks.com/matlabcentral/fileexchange/10676-circular-statistics-toolbox-directional-statistics; RRID: SCR_016651
SPM	Wellcome Department of Cognitive Neurology, London, UK	https://www.fil.ion.ucl.ac.uk/spm/; RRID: SCR_002823
FSL	FMRIB Software Library	https://fsl.fmrib.ox.ac.uk/fsl/fslwiki/; RRID: SCR_002823
AFNI	SSCC at NIMH	https://afni.nimh.nih.gov/; RRID: SCR_005927