Evolutionarily conserved fMRI network dynamics in the mouse, macaque, and human brain

Dominant coactivation modes parsimoniously describe fMRI network dynamics in humans, macaques, and mice

To compare the dynamic organization of fMRI network dynamics within an evolutionary perspective, we used a frame-wise approach based on the identification of coactivation patterns (CAPs)^15,17. Data for this study consisted of awake rsfMRI datasets from two human cohorts (The Hangzhou Normal University - HNU: 30 subjects, 10 test-retest sessions²² and The Midnight Scan Club - MSC: 10 subjects, 5 test-retest sessions²³); the Newcastle (NC) macaque cohort^24,25 (8 animals with two test-retest sessions); and 44 mice under head-fixed conditions. Human data was sex balanced in both independent datasets, while all results for macaques and mice were obtained using male animals only. Using k-means, we clustered fMRI frames in the concatenated timeseries given their spatial similarity. We averaged the frames in each cluster to produce group-level CAP maps^15–17, then T-scored them at the voxel level. We then obtained single-subject CAP maps by averaging the corresponding clustered frames from each subject (Fig. S1). To identify reproducible functional CAPs that are representative of the dynamic structure of fMRI networks in each species, we used a multi-criteria approach (see methods section) aimed at maximizing CAP reproducibility across individuals, sessions, or datasets. This approach yielded k = 8 CAPs in humans and macaques, and k = 6 in mice as optimal clustering solutions (see methods section and Figs. S2, S3). Previous work¹⁶ has shown that CAPs embody rich fMRI topographies that can be reliably matched into mirrored coactive and anti-coactive pairs characterized by opposite patterns of fMRI coactivation (Fig. 1). As we will demonstrate below, and similar to what has been observed in quasi-periodic patterns²⁶, each matched CAP and anti-CAP pairs describe a cyclical fluctuation of a single fMRI state. As predicted, the employed procedure identified in all species mirrored CAP pairs characterized by opposite fMRI coactivation (spatial correlation, r < −0.65, all species, all CAP pairs, Fig. 1A, B). Importantly, the identified CAPs explained in all species a large proportion of variance in fMRI timeseries (R² > 0.58 mice, >0.68 macaques, and >0.75 humans). Moreover, they also guaranteed robustness within and between dataset spatial and occurrence rate reproducibility (Figs. S2, S3, permutation tests with random CAP-identity shuffling). These results suggest that a few dominant dynamic patterns can parsimoniously describe the dynamic organization of fMRI activity in multiple mammalian species.

To further reduce dimensionality and facilitate cross-species comparisons, we coalesced highly anticorrelated CAP pairs into a single Coactivation Mode (C-mode). CAP and anti-CAP are thus the peak and trough of the same fluctuating C-mode. C-modes are thus computed by reversing the sign of the anti-CAP, and spatially averaging it with its paired CAP such to spatially depict the corresponding coactivation axis (Fig. 1A). The topography of the individual CAPs pairs that constitute each C-mode is reported in Figure S4. To assess the potential confounding contribution of head motion to C-mode topography mapping, we repeated the clustering procedure on time-series in which we did not scrub fMRI frames exhibiting high motion-related. We next compared the spatial topography of C-modes obtained with and without frame censoring. This comparison revealed that the C-mode obtained using frame-censored timeseries were very similar to those we mapped using the entire timeseries (r > 0.95, all C-modes, all species, Fig. S5A). Further corroborating a negligible contribution of head motion to our findings, we also found that no C-mode was preferentially enriched with high-motion frames in all species (one-way ANOVA, p > 0.3, F < 1.14 for all comparisons, Fig. S5B).

fMRI C-modes exhibit evolutionarily conserved functional organization

Having identified dominant fMRI C-modes in all three species, we next asked whether their functional organization would show, on top of foreseeable species-specific features, recognizable evolutionarily conserved anatomic features. For this purpose, we matched C-mode topographies based on the similarity of mean fMRI coactivation profile across a set of evolutionarily conserved resting-state networks (RSNs)^7,27–30 (Fig. 1C). The chosen networks include the default mode (DMN), visual (VIS), somatomotor (SMN), limbic (LIMB), ventral (VAT), and dorsal (DAT) attention, frontoparietal (FPN), as well as key subcortical nuclei of the thalamus (TH) and hippocampus (HCP) and were selected based on the notion that, across these three species, they encompass partly-conserved neuroanatomical substrates^7,10. A key exception of note is the lack of established phylogenetic precursors of the VAT, DAT, and FPN in the rodent brain¹⁰. For this reason, these networks were not included in the coactivation profile in mouse data. For matching purposes, we ordered C-modes for each species in decreasing joint occurrence rates, and labeled them H1-H4 for humans, M1-M4 for macaques, and R1-R3 for rodents (mice). To allow spatial comparisons between species, C-mode maps were z-scored spatially, and the mean of the normalized activity of voxels within each RSN mask was computed to build its corresponding profile vector (Fig. 1D, C-modes numbered by decreasing occurrence rate). Vectors were first matched using the Hungarian Algorithm³¹ from human C-modes (organized in decreasing occurrence rate) to macaques, then from macaques to mice.

The corresponding results are depicted as the topological correlation between network coactivation profiles of human and macaque datasets (HNU and MSC, Fig. 1E, F, respectively), as well as between macaques and mice (Fig. 1G). Human-to-macaque matching gave consistent results for both HNU and MSC datasets, with the spatial topography of human C-modes H1-H4 being best aligned to corresponding macaque C-modes M3, M4, M2, and M1, respectively (see also Fig. 2A). Because the human HNU dataset included more subjects and was performed at daylight hours as animal scans, we describe our results hereafter for the HNU as main dataset and present a summary of the obtained MSC results as Supplementary Fig. (Figs. S6). Macaque to mouse matching linked C-modes M1 to R1, M2 to R2, and M4 to R3, respectively (Fig. 1G,  2A). The resulting overall cross-species matchings yielded an organization linking C-modes H1-M3, H2-M4-R3, H3-M2-R2, and H4-M1-R1. For brevity, we henceforth refer to these as C-mode 1, C-mode 2, C-mode 3, and C-mode 4, respectively. As human and macaque C-mode 1 was not preferentially matched to any mouse C-modes, we refer to mouse C-modes as 2, 3, and 4 for consistency with those mapped in higher species throughout the manuscript.

Collectively, cross-species C-mode matching revealed four topographically related C-modes in humans and macaques, three of which (C-modes 2–4) were also represented in the rodent brain (Fig. 2A). In all species, C-modes exhibited rich spatial organization encompassing positive and negative peaks of fMRI activity that together delineated a set of stereotypic network coactivation profiles (Fig. 2A, right panels). Specifically, C-mode 1 encompassed peaks of fMRI activity in DMN, accompanied by below-baseline activity in somatomotor and ventral attention areas in both humans and macaques. C-mode 2 exhibited coactivation of DMN areas in anticorrelation with somatomotor, visual, and limbic networks in all three species. Interestingly, while in humans, the hippocampus had below-baseline fMRI activity, macaques and mice showed above-baseline fMRI activity in this region. C-mode 3 activity peaked in the somatomotor areas and concomitantly engaged most cortical regions of the human brain, albeit with considerably weaker or negative fMRI activity in visual, thalamic, and basal forebrain areas, reminiscent of global fMRI signal (GS) fluctuations. Finally, C-mode 4 was characterized in humans and macaques by positive coactivation in visual and posterior cortical regions, and negative coactivation in frontoparietal cortical regions. This topography in mice and macaques, but not humans, was associated with co-deactivation of the anterior cingulate and prefrontal regions of the DMN. Taken together, these results point to the presence of notable topographic correspondences in the functional organization of dominant C-modes in human, macaque, and mouse brains.

We next computed for each C-mode its occurrence rate, defined here as the proportion of fMRI frames assigned to each C-mode for each subject/animal. Interestingly, while the spatial organization of C-modes exhibited species-invariant topographic features, their occurrence rate, showed variation across species (Fig. 2B). Specifically, in humans, we observed a dominant occurrence of C-mode 1 and 2 (Kruskal–Wallis test, p < 0.001) while mice occurrence of sensory-oriented C-mode 3 and 4 was observed instead (Kruskal–Wallis test, p < 0.001), with macaques showing equiprobable (Kruskal–Wallis test, p = 0.25) C-mode occurrence (although a slight trend for a mouse-like profile was apparent). Thus, the temporal structure of C-modes in humans is a preferred occurrence of configurations led by associative cortical networks such as the DMN (C-modes 1 and 2), with mice and possibly macaques exhibiting instead greater occurrence of C-modes encompassing coactivation of sensory unimodal areas (C-modes 3 and 4). To further corroborate these results, we additionally computed the occurrence rates of individual CAPs (Fig. S7) and compared them with their corresponding anti-CAP (both of which compose a C-mode), finding no significant differences between these pairs of distributions in any of the analyzed datasets (p > 0.2, signed-rank test). This result supports the notion that C-modes reliably encapsulate key dynamic properties of their constituting CAPs.

fMRI C-modes exhibit infraslow fluctuations in humans, macaques, and mice

Our formalism allowed us to compare C-Modes with single fMRI frames and thus uncover the temporal structure of the spatial configurations that underlie spontaneous fMRI dynamics. We leveraged this property to describe, in all species, the temporal evolution of each C-mode for each subject/animal by computing the instantaneous spatial correlation between each C-mode map and each fMRI frame in all timeseries. The power spectra of the resulting C-mode timeseries revealed that C-modes undergo infraslow fluctuations in all species, with most of the power peaking within the 0.01–0.03 Hz range (Fig. 3). Power within the infraslow band was significantly greater than the mean power of surrogate timeseries generated from shuffling the joint C-mode timeseries 1000 times (Fig. 3, black traces). Peaks of infraslow activity were distinct and sharp in humans, and slightly less prominent (yet clearly recognizable) in mice and macaques.

We next investigated the assembly of each C-mode by normalizing (z-score) its time-course and by time locking peaks of C-mode time-courses (i.e., event-wavelet) at the group level. We found that C-modes in awake humans, macaques, and mice assemble and disassemble in a slow and gradual fashion (Fig. 3, red insets), reminiscent of damped oscillations in the dominant frequency band. These results suggest that fMRI C-modes exhibit comparable infraslow dynamic cycling in all the examined mammalian species.

fMRI C-modes occur at specific phases of fMRI global signal cycles in humans, macaques, and mice

Prompted by previous investigations in anesthetized mice^16,29, we next probed if C-modes in higher species have a preferred occurrence within fMRI global signal (GS) cycles. We thus first ascertained that also in awake conditions, the GS would dominantly fluctuate within the infraslow range in all species (Fig. 4A). A spectral analysis showed that the power spectrum of GS sharply peaks within a 0.01–0.03 Hz band in awake humans, with analogous (albeit less pronounced) peaks of activity in the same frequency range in both awake macaques and mice. We then built a circular distribution of the phases of the filtered (0.01–0.03 Hz) GS at which each C-mode occurred, by sampling only occurrences in which the normalized C-mode timeseries surpassed threshold values of 1 SD. C-mode occurrence was significantly phase-locked within GS-cycles in humans, macaques, and mice (Fig. 4B, Rayleigh test, p < 0.05, FDR corrected). Interestingly, GS-phase distributions in macaques and humans presented key similarities, with C-modes 1, 3 and 4 (but not 2) exhibiting remarkably conserved cross-species phase alignment. Moreover, C-modes 3 and 4 showed broadly similar circular means across all the three species examined here.

C-modes were computed by coalescing a CAP with an anti-CAP of nearly identical spatial shape but opposite fMRI coactivation. This description of brain dynamics works well under the assumption that CAP and anti-CAPS represent a single fluctuating brain sub-state whose pattern of activity changes sign cyclically. To corroborate the cyclic nature of C-mode fluctuations within GS-cycles, we repeated the sampling of GS-phases but with inverted (negative) C-mode time-courses so to capture the full temporal evolution of C-mode. This approach revealed clearly opposite distributions of sampling of the GS-phases from positive and negative C-mode occurrences (Fig. 4B, red insets). These results suggest that C-modes describe cycling spatiotemporal sub-state that fluctuates in magnitude and sign according to the infraslow structure of the global fMRI signal.

To further provide evidence of the infraslow oscillatory nature of C-mode dynamics, we investigated whether different C-modes were phase-coupled within individual GS-cycles. We computed the GS-phase difference between the occurrences of a C-mode in a GS-cycle, and the following occurrences of another C-mode either within the same GS-cycle, or in an immediately subsequent cycle (Fig. S8). These analyses showed highly consistent positioning of a C-mode within each GS-cycle (Fig. S8, diagonals) and phase coupling between some C-mode pairs (see C-mode 2 and 3 and 3 and 4 in humans, and C-modes 3 and 4 in mice). These results suggest that in mammalian species, intrinsic fMRI signal fluctuations do not reflect spatially undifferentiated peaks and troughs of fMRI activity, but instead encompass infralow cycling between dominant patterns of fMRI activity. This general principle can be extended to entail the evolutionary conservation of the phase relationship with fMRI GS-cycles for most (albeit not all, see C-mode 2) of the explored C-modes. In sum, these results show that C-modes are phase-locked to intrinsic GS fluctuations and imply that different C-modes can be conceptualized as networks of coupled oscillators in multiple mammalian species.

Coupled oscillatory activity explains C-mode transition dynamics

We next considered whether the temporal structure of C-mode instantaneous transitions could similarly be underpinned by species-invariant principles. For each species, we modeled the system as a Markov process from sequences of concatenated C-mode occurrences and computed the transition probability into a different C-mode, as well as C-mode self-transitions (also termed persistence probability, Fig. 5A). We found that the most recurring C-modes (1–2 in humans, 3–4 in macaques, and 3–4 in mice) were those with the highest probability of being transitioned to (p < 0.01, black crosses in Fig. 5A). Thus, they can be conceptualized as sinks of preferred directional transitions.

Furthermore, we investigated, for each species, the accessibility of C-modes from one another by computing the corresponding Entropy of Markov Trajectories (HMT)³² from the transition probability matrices. This parameter measures the complexity of a transition: low entropy values imply an almost deterministic direct path or high accessibility. On the contrary, high entropy values suggest high uncertainty, requiring random steps through different C-modes before reaching the destination, i.e., low accessibility. In keeping with the C-mode occurrence rates we describe in Fig. 2B, the most recurring C-modes in all species were also those that were most accessible, i.e., they were the ones characterized by the lowest entropy values (Fig. 5B, C). Specifically, in humans, C-modes 1 and 2 were the most accessible ones (p < 0.0001 against C-modes 3–4). In macaques, we did not find any C-mode to have preferred accessibility (p = 0.25) Conversely, in mice, C-modes 3 and specially 4 were the most accessible ones (p < 0.0001). Importantly, in all species, the ensuing accessibility profile (Fig. 5C) recapitulated the C-mode occurrence rates we described in Fig. 2B, with the most accessible C-modes being also the most occurring ones (Fig. 2B). This was true also in macaques, where the entropy of Markov trajectories required to reach C-modes 3 and 4 followed the tendency of higher occurrence rates of these C-modes.

In the sections above, we characterized intrinsic fMRI dynamics both in terms of coupled quasi-periodic fluctuations between coactive networks captured by different C-modes, and also in terms of transitions between C-Modes. If coupled cyclic activity is key to explaining state transition dynamics, we expect that trajectories between C-modes with lower entropy (i.e., small HMT, corresponding to more direct transitions) would occur on average with shorter infraslow phase differences. Conversely, C-mode pairs with larger entropy (i.e., higher HMT corresponding to less direct transitions) would occur on average with longer infraslow phase differences. To test this hypothesis, we computed the circular-linear correlation between the mean GS-phase differences between C-mode occurrences (depicted in Fig. S8) and the HMT for their trajectories. We found clear and significant correlation values in all species (r = 0.66, 0.66, 0.87 with p = 0.029, 0.031, and 0.032 for humans, macaques, and mice, respectively). This implies that the transition structure of C-modes is, in part, described by infraslow-coupled cyclic dynamics across all species, embedding fast transition phenomena within GS-cycles in the dominant infraslow band.

fMRI C-mode occurrence predicts the ranking of connectivity gradients in humans, macaques, and mice

Previous investigations have shown that a high portion of the variance in static fMRI connectivity is explained by a limited fraction (5–15%) of fMRI frames exhibiting exceedingly high co-fluctuation amplitude^16,33,34. We thus investigated whether the dynamic occurrence of dominant C-modes alone could similarly be sufficient to reconstitute key organizational features of the static fMRI connectome. To this aim, we first examined whether the occurrence-weighted average of C-modes would reproduce the static architecture of fMRI connectivity. To this purpose, we calculated for each C-modes its co-fluctuation matrix, i.e., the cross-multiplication of each C-mode map with itself ²⁹. We found that in all species, the weighted average by occurrence rate of C-mode co-fluctuation matrices exhibited a high correlation with the corresponding group-mean static fMRI connectivity matrix (r > 0.57, all species, with r = 0.75 in humans, Fig. S9). These findings corroborate the notion that C-modes dynamics account for high co-fluctuation events critical for the topographic organization of the static functional connectome.

To further investigate the relationship between C-mode dynamics and static fMRI connectivity, we next inquired whether C-mode occurrence could also be related to the organization of fMRI connectivity gradients^20,35. Gradient approaches were favored here because of their ability to chart the high-dimensional static functional connectome into low-dimensional manifold representations. In humans, the principal gradient of resting-state fMRI delineates a hierarchical organizational axis that reconstitutes the unimodal-transmodal spectrum in cortical regions. Crucially, this arrangement spatially recapitulates key myeloarchitectural, cytoarchitectural, and transcriptional features of the human and animal neocortex^36,37. Higher-order gradients have similarly been linked to axes of differentiation between sensory cortical modalities. Here, we examined how the dynamic organization of fMRI activity relates to the spatial organization of gradients. We posited that C-mode occurrence rate may be linked to the ranking of functional connectivity gradients, with dominant (i.e., most occurring) C-modes aligning with the gradients that explain the most variation in fMRI connectivity. To test this hypothesis, we first computed for each species the top five fMRI connectivity gradients, which we next ranked by decreasing variance explained (Fig. S10). We then compared, for each species, the spatial similarity of the obtained gradients to that of each C-mode, matching them according to their highest absolute spatial correlation. Supporting our hypothesis, plausible spatial correspondences between C-mode and gradient topographies were observed in all species (Fig. 6A). Specifically, in humans, the most occurring C-modes 1 and 2 were matched with dominant gradients 1 and 2, while in macaques and mice, dominant gradients 1 and 2 were matched with most occurring C-modes 4 and 3. Moreover, in all species, C-mode occurrence was linearly related to the variance explained by each gradient (R² > 0.78 for all species, Fig. 6B). This analysis shows how recurring network interaction represented by dominant C-modes, and their relative occurrence rate, shape the organization of the static fMRI connectome and its principal axis of variance in all the probed species.

Data and preprocessing

Resting-state fMRI (rsfMRI) datasets from awake, freely breathing humans, macaques, and mice were used in this study. Human and macaque datasets were previously published (see references below). Preprocessing included most steps suggested by the guidelines of the Human Connectome Project⁷⁴, using a combination of fMRI-dedicated software AFNI⁷⁵, FSL⁷⁶, FreeSurfer⁷⁷, and SPM12 (http://fil.ion.ucl.ac.uk/spm/).

Human

The main dataset, the Hangzhou Normal University of the Consortium for Reliability and Reproducibility (CoRR-HNU, or HNU)²² includes 10 sessions of 10-min scans over the course of a month from n = 30 young healthy adults with no history of neurological or psychiatric disorders, head injuries, nor substance abuse (balanced sexes, age = 24 ± 2.41 years). Before scanning, participants were asked to relax and remain still with their eyes opened, avoiding falling asleep. During scanning, a black crosshair was shown in the middle of a gray background. The study was approved by the ethics committee of the Center for Cognition and Brain Disorders at Hangzhou Normal University, and all participants signed written consent before data collection. A GE MR750 3T scanner (GE Medical Systems, Waukesha, WI, USA) was used to acquire MRI data. Functional scans were acquired with an echo-planar imaging sequence - EPI: TR = 2 s, TE = 30 ms, flip angle = 90°, FOV = 220 × 220 mm, matrix = 64 × 64, voxel-size = 3.4 mm isotropic, 43 slices. Data were downloaded from the International Neuroimaging Data-Sharing Initiative (INDI - http://fcon_1000.projects.nitrc.org/indi/CoRR/).

The Midnight Scan Club (MSC) dataset²³ was used as a secondary, replication dataset. 5 out of 10 randomly selected sessions were used on separate days, and included 30-min scans from n = 10 healthy young adults (balanced sexes, age = 29.1 ± 3.3 years). Participants were asked to visually fix on a white crosshair against a black background. The study was approved by the Washington University School of Medicine Human Studies Committee and Institutional Review Board, and all participants signed written consent before scanning. Functional scans were acquired with a Siemens TRIO 3T MRI scanner (Erlangen, Germany) using a gradient-echo EPI sequence: TR = 2.2 s, TE = 27 ms, flip angle = 90°, voxel-size = 4 mm isotropic, 36 slices. Data were downloaded from OpenNeuro (10.18112/openneuro.ds000224.v1.0.3).

Preprocessing. The first 5 fMRI volumes were removed from each subject’s raw data, then despiking (AFNI 3dDespike) and slice-timing correction (AFNI 3dTshift) was performed. Data subsequently underwent motion-correction (AFNI 3dvolreg); skull-stripping (FSL fast and bet⁷⁸); co-registration to the MNI 3 mm isotropic template (FSL flirt); regression of nuisance parameters (white matter, cerebrospinal fluid, and 24 motion parameters (6 parameters, 6 derivatives, and their respective squared time-series) (AFNI 3dDeconvolve); band-pass filtering between 0.01–0.1 Hz (AFNI 3dBandpass); spatial smoothing with a 6 mm FWHM kernel (AFNI 3dBlurInMask); and voxel time-series were finally normalized to z-scores (zero-mean, and standard deviation units).

Macaque

The Newcastle dataset (NC) includes n = 14 rhesus macaque monkeys (Macaca Mulatta) scanned with no contrast agents, from which N = 10 animals (two females, eight males), in which two independent fMRI sessions were available, were used for our analyses (age = 2.28 ± 2.33, weight = 11.76 ±  3.38). Two additional (female) animals were discarded from our analyses due to excessive head motion exceeding, in either session, over 30% of fMRI volumes with Frame-wise Displacement above a 0.3 mm threshold. Animal procedures, head-fixation, and protocols were approved by the UK Home Office and comply with the Animal Scientific Procedures Act (1986) on the care and use of animals in research and with the European Directive on the protection of animals used in research (2010/63/EU) (see ref. ²⁴ for protocol specifics on animal preparation for awake imaging²⁴). A Vertical Bruker 4.7T primate dedicated scanner was used, and rsfMRI experiments were performed in awake, head-fixed animals for two separate sessions with TR = 2 s; TE = 16 ms, voxel-size = 1.2 mm isotropic. Data were downloaded from NHP data-sharing consortium PRIME-DE (http://fcon_1000.projects.nitrc.org/indi/indiPRIME.html)²⁵.

Preprocessing. The first 5 rsfMRI volumes were removed from each animal’s raw data, then despiking (AFNI 3dDespike) and slice-timing correction (AFNI 3dTshift) was performed. Data subsequently underwent motion-correction (AFNI 3dvolreg); skull-stripping (FSL fast and bet⁷⁸); co-registration to the Yerkes19 2 mm isotropic template⁷⁹ (FSL flirt); regression of nuisance parameters (white matter, cerebrospinal fluid, and 24 motion parameters (six parameters, six derivatives, and their respective squared time-series) (AFNI 3dDeconvolve); band-pass filtering between 0.01–0.1 Hz (AFNI 3dBandpass); spatial smoothing with a 3 mm FWHM kernel (AFNI 3dBlurInMask); and voxel time-series were finally normalized to z-scores (zero-mean, and standard deviation units).

Mouse

C57BL/6J mouse data were obtained from n = 44, head-fixed awake male mice undergoing a 12-min rsfMRI scan²⁹. In vivo experiments were conducted in accordance with the Italian law (DL 26/214, EU 63/2010, Ministero della Sanita, Roma) and with the National Institute of Health recommendations for the care and use of laboratory animals. The animal research protocols for this study were reviewed and approved by the Italian Ministry of Health and the animal care committee of the University of Trento, and Istituto Italiano di Tecnologia (IIT). All surgeries were performed under anesthesia. Young adult (<12 months old) male C57BL/6 J mice were used. RsfMRI scans, both retrieved and newly acquired, were acquired at the IIT laboratory in Rovereto (Italy) using a Bruker 7T scanner (Bruker Biospin, Ettlingen) with a BGA-9 gradient set, 72 mm birdcage transmit coil, and a four-channel solenoid receiver coil: TR = 1 s, TE = 15 ms, flip angle = 60°, matrix = 100 × 100, FOV = 2.3 × 2.3 cm, 18 coronal slices 0.6 mm thick, 12 min total acquisition time.

Preprocessing. The first 120 rsfMRI volumes (2 min) were removed from each animal’s raw data to account for thermal gradient equilibration, then despiking (AFNI 3dDespike) was performed. Due to the short TR we did not perform slice-timing correction in these images. Data subsequently underwent motion-correction (FSL mcflirt); skull-stripping (FSL fast and bet⁷⁸); co-registration to an in-house mouse brain template of 0.23 × 0.23 × 0.6 mm³ (ANTS registration suite⁸⁰); regression of nuisance parameters (white matter, cerebrospinal fluid, and 24 motion parameters (6 parameters, 6 derivatives, and their respective squared time-series) (AFNI 3dDeconvolve); band-pass filtering between 0.01–0.1 Hz (AFNI 3dBandpass); spatial smoothing with a 0.5 mm FWHM kernel (AFNI 3dBlurInMask); and voxel time-series were finally normalized to z-scores (zero-mean, and standard deviation units).

Whole brain CAP detection and cluster-number selection

To identify recurrent rsfMRI whole-brain states, we used the whole-brain coactivation patterns (CAPs) approach^15–17 in which fMRI frames are cluster based on their spatial similarity and then averaged to define recurrent patterns of fMRI coactivation. Specifically, for each species, we first performed censoring of motion-contaminated frames (frame-wise displacement: FD >0.3, 0.3, and 0.075 mm for humans², macaques²⁸, and mice²⁹ respectively), and then concatenated the frames from all subjects or animals. Given that clustering human fMRI data in the 3 mm MNI-space became computationally challenging owing to its large dimensionality (n-voxels = 43,539, compared to 11,402 in macaques, and 8937 in mice), this step was carried out upon reducing data using the coarse 950-ROI Craddock Parcellation⁸¹. The choice of this specific parcellation regards its ability to provide a fair dimensionality reduction, while preserving information present at the voxel scale⁸¹. After these final steps, we ran, for each species, the k-means clustering algorithm^16,82 (spatial correlation as the distance metric, 500 iterations, 5 replications with different random initializations, from k = 2:20, 5 independent runs). CAP maps were obtained at the group level by averaging the fMRI frames belonging to a cluster at the voxel level, then normalizing these values to T-scores from the concatenated datasets (Fig. S1). At the single-subject level, we obtained, for each subject/animal, a CAP map for each cluster by only averaging and converting to T-scores, the frames belonging to a cluster but only within a subject/animal’s data/session. We note that CAP mapping through frame averaging in humans was done at the voxel level, as the parcellated data was only used for clustering purposes. After clustering, we recovered the censored fMRI frames to the cropped datasets and assigned them to the CAP with the highest spatial correlation. This was done in order to have a continuum of frames for subsequent analyses.

Selection of the optimal number of clusters was done following a set of previously proposed empirical rules^16,19, as well as newly proposed metrics. These were dependent on the availability of test-retest sessions within a dataset in macaques and humans, as well as a full independent dataset in humans from a different site. Specifically, for human HNU and MSC datasets independently, we first ran the k-means clustering algorithm with the concatenated dataset from k = 2:20, selecting, for each partition, the solution with the highest variance explained¹⁶ from five replications in five independent runs. Here, the variance of the data explained by each partition is defined as 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^16,83. For each dataset independently, we first computed the variance explained by data partitioned into an increasing number of clusters (2 ≤ k ≤ 20) (Fig. S2A, B). We next assessed the topographic consistency of CAPs at increasing partitions. To this purpose, we computed how consistent CAPs are by assessing, for each k partition, the spatial correlation between a mean CAP map, and its matched map in the previous order (k-1) partition. Matching was done using the Hungarian Algorithm³¹. We found that partitions between k = 6:10 yielded, in both datasets, topographically stable CAPs that could be reliably identified at higher partitions (Fig. S2C, D).

We then assessed the within-dataset repeatability by comparing the spatial correlations between the mean CAP maps of each subject between each independent fMRI sessions (10 for HNU, and 5 for MSC, Fig. S2E, F). Statistical significance of the mean within-subject repeatability of each CAP was assessed by recomputing the spatial correlations between subject-level CAP maps after randomly shuffling the CAP-identity of fMRI frames, preserving occurrence rates. This process was repeated 1000 times, and repeatability values for each subject beneath the highest permutation value were flagged as non-repeatable (asterisks in Fig. S2E, F). The result of these comparisons showed that within an upper limit of k = 13 (HNU) or 15 (MSC), all the mapped CAPs were represented in all fMRI sessions of each subject, with significant spatial correlation across sessions (p < 0.05, surrogate testing with randomly shuffled cluster associations). At higher partitions, one or more CAPs were instead no longer represented in one or more subjects. We further probed the within-session stability of clustering by computing the CAP occurrence rate obtained across imaging sessions. We found that the CAP occurrence rate (i.e., the proportion of fMRI frames associated with a CAP in each subject), for both datasets, was stable across sessions at k = 2, 6, 8, and 9 (Kruskal–Wallis test, ten groups for HNU, five groups for MSC, FDR corrected for k comparisons, Fig. S2G, H).

Finally, to maximize the generalizability of our partitioning, we compared the main dataset’s (HNU) CAP maps with those obtained in the MSC in terms of topography matching and occurrence rates. To this aim, we spatially compared mean group-level CAP maps (matched with the Hungarian Algorithm for each partition) and tested the significance of this comparison by recomputing the values after randomly shuffling the fMRI frames within each dataset, while preserving occurrence rates (Fig. S2I). The mean occurrence rates from each dataset was assessed with a Wilcoxon signed-rank test, p < 0.05, FDR corrected for k comparisons (Fig. S2J). This analysis revealed that clusters from all partitions, except k = 6, 7, 9, 14, 15, and 16, were topographically reproducible across datasets (Fig. S2I). By contrast, a comparison of the mean CAP occurrence rates between datasets (Fig. S2I) showed that k = 8 was the only partition in which this parameter was conserved across datasets (Wilcoxon signed-rank test, p < 0.05, FDR corrected for k comparisons). Based on these analyses, k = 8 was the only partition level meeting all the required within and between subject/dataset reproducibility criteria. We thus based all our subsequent analyses of human fMRI dynamics using k = 8 clusters.

Selection of optimal clusters in macaque fMRI time-series (Newcastle test-retest datasets, n = 8 animals, two sessions) followed the same strategy employed for human time-series. Computation of explained variance in macaque fMRI time-series revealed an elbow region within the k = 6:10 range (Fig. S3A). Macaque CAPs were topographically stable at increasing k partitions until an upper limit at k = 10 (Fig. S3B). Comparing the repeatability of CAP topographies across sessions revealed that above k = 8, some CAPs present a lack of topographical reproducibility (p > 0.05, surrogate testing with randomly shuffled cluster associations), as well as significant differences in CAP occurrence rates (Kruskal–Wallis test, FDR corrected for k comparisons, Fig. S3C, D). This cumulative evidence suggests that k = 8 was the highest partition that guarantees CAP stability, as well as test-retest topographical and frame distribution repeatability amongst clusters.

For mice, the variance explained curve for the awake dataset we used in this work (N = 44) showed an elbow between 6 and 8 clusters (Fig. S3E, F). Topographic CAP stability as a function of increasing partition number revealed that mouse CAPs were topographically stable up to k = 6. This value is in agreement with the results of CAP number selection in prior independent studies^29,39, in which k = 6 was consistently identified as the optimal partition in this species.

Coactivation modes and between-species matching

In previous work we demonstrated that CAPs appear in pairs of spatial configurations with opposite fMRI coactivations¹⁶ (CAP, and anti-CAP, Fig. 1A). As this feature was initially detected in awake and anesthetized mice²⁹, we first showed that these anatomical configurations were present also in all of our datasets. For each species, we thus organized CAPs in descending occurrence rate and computed the Pearson’s correlation between the vectorized mean CAP maps as a topological similarity metric, matching CAPs and anti-CAPs as pairs with the highest anticorrelation coefficient (Fig. 1B). We then reduced pairs into C-modes, by taking each CAP map pair (T-score normalized), subtracting the least recurring one from its counterpart, then dividing the resulting map by two (Fig. 1A). In order to have a normalized basis for between-species comparisons, we further z-scored each resulting C-mode map. Given that C-modes represent unique configurations between known resting-state networks^16,17,19, we extracted the mean fMRI values from voxels within masks of evolutionarily conserved^30,53 cortical RSN in humans²⁷, macaques²⁸ and mice^29,69, as well as hippocampal and thalamic masks^74,84 (Fig. 1C). To best match DMN organization across species, we split this network into an anterior and a posterior component. This choice was made to account for recent primate and mouse evidence of rostrocaudal disassembly of this network in some datasets and species^29,58. In rodents, this is especially apparent in awake conditions, where parietal and peri-hippocampal cortical components that are synchronously connected to the mouse PFC (in anesthesia) tend to form dissociable networks in awake animals²⁹. For consistency with previous studies, we included the caudate-putamen in the mouse DMN (but not in primate and human DMN) as in rodents, this region is proximal and highly synchronous with anterior cingulate hubs of the DMN both under anesthesia and in awake conditions¹². Pilot studies showed that the inclusion (or lack thereof) of the CPu in mouse DMN did not alter monkey-to-human C-Mode matching.

For each C-mode, we computed their occurrence rates as the average occurrence of the CAP and anti-CAP conforming to the C-mode. We then organized the vectorized network coactivation profiles of each C-mode in descending occurrence rate order, and matched C-mode profiles of humans (C-mode H1-4) to macaques (C-mode M1-4), then macaques to mice (C-mode R1-3, rodent) using the Hungarian Algorithm³¹. This matching was conducted using vector-correlations between profiles as distance metrics for the cost function to be minimized (Fig. 1D). We represented the similarities between network coactivation profiles of C-modes in the human HNU dataset and macaques (Fig. 1E), as well as with the MSC dataset and macaques, confirming our matching (Fig. 1F). Finally, we matched macaque C-mode network coactivation profiles with those of mice, without accounting for frontoparietal (FPN), ventral (VAT) and Dorsal Attention (DAT) networks, as these have not been described in the mouse brain (Fig. 1G). We note that since C-mode 1 was not reliably identified in mice, most likely due to the lack of weighting by higher-order VAT, DAT, and FPN, we hereafter exclude C-mode 1 from this species.

Because C-modes fundamentally include an anticorrelated counterpart, in pilot studies, we probed cross-species matching using absolute value correlations as opposed to positive correlation only. We found that only macaque M4 and mouse R3, as well as M1 and the inverse of R1, produced a different match when absolute value correlation was used as distance metrics. Matching between human and macaque C-modes remained unchanged. However, while this alternative matching appeared to be spatially plausible, only matching using positive spatial correlation as a cost function resulted in a plausible cross-species preservation of phase coupling with fMRI global signal cycles in macaques and mouse. For this reason, we retained positive correlation-based matching as a reference dataset for the present manuscript, as we believe it provides a more plausible representation of key evolutionarily relevant dynamic features of C-modes.

After matching and organizing C-modes, we mapped the voxel-wise C-modes as well as the network coactivation profiles within C-modes (within-network voxel mean ± SD, Fig. 2A). Independent group-level CAP maps for each human (HNU and MSC); macaque and mouse dataset are shown in Fig. S4). For each subject/animal in each dataset, we computed a C-mode’s occurrence rate as the proportion of fMRI frames associated with the CAP pairs composing the C-mode. A comparison of the distributions of occurrence rates was then performed (Kruskal–Wallis test, FDR corrected).

We further tested that in-scanner head motion did not selectively affect the configuration of any C-mode. This was done by first independently clustering the non-motion censored rsfMRI datasets with k = 8, 8, and 6 for humans, macaques, and mice, respectively; building C-modes; and then comparing the topographies of matched C-mode maps (Hungarian Algorithm³¹) to the ones obtained after censoring (Pearson’s Correlation) (Fig. S5A). We then computed, for each C-mode in each subject/animal, the proportion of frames in it that were classified as high-head motion (FD >0.3, 0.3, and 0.075 mm for humans, macaques, and mice, respectively), and compared these distributions across C-modes (one-way ANOVA, Fig. S5B).

C-mode infraslow dynamics, formation, and temporal structure

We investigated the temporal evolution of C-modes by generating a time-course of instantaneous C-mode to fMRI frame spatial correlation (hereafter referred to as C-mode time-courses) for each subject/animal, and then computing the power spectrum of these time-courses as well as that of the Global fMRI signal (GS, i.e., the average of all voxels at each frame). We next computed the group mean ± SEM power spectra. C-mode time-courses were normalized to standard deviation units and mean zero. The significance of the observed infraslow portion of C-mode spectra was assessed by calculating the power spectrum of surrogate C-mode timeseries obtained after jointly shuffling the timepoints of the original C-mode time-courses using the same random sequence of time points. This method preserves the signal amplitude distribution and the correlation structure of the original timeseries¹⁶. We next carried out, for each frequency, a one-tailed two-sample t-test (FDR corrected for comparisons within the 0.01–0.1 Hz band) comparing the power of the original C-mode PSD with the mean of 1000 iterations of the surrogates. The evolution of the formation of a C-mode was computed by time-lock averaging the C-mode to frame correlation values at the vicinity of a peak surpassing 1 SD (±30 s), then averaging these event-kernels at the group level (Fig. 3 insets).

Given that the GS also presented a strong peak in the power spectrum in the [0.01–0.03 Hz] infraslow band (Fig. 4A), we computed whether the occurrence of C-modes was phase-locked to infraslow GS fluctuations^16,29. We extracted the instantaneous infraslow phase of the GS (filtered between 0.01–0.03 Hz) using the Hilbert Transform, then divided the trace into cycles of minimum 30 s and maximum 100 s. For each cycle, if a C-mode was present, we sampled the GS-phase at which the C-mode occurred, retaining only samples from normalized C-mode time-courses exceeding 1 SD, thus guaranteeing the selection of frames that are reasonably well assigned to a specific C-mode¹⁶. We then built the distribution of GS-phases from the group-concatenated sampling of C-mode occurrences, and using MATLAB’s CircStats toolbox⁸⁵, we performed a Rayleigh test (p < 0.001, FDR corrected) for deviations from circular uniformity (Fig. 4B). To confirm the cyclical features of C-mode fluctuations within GS-cycles, we repeated this analysis but sampling GS-phases from occurrences of inverted (negative) C-mode time-courses (Fig. 4B, red insets). As in our previous work¹⁶, we also demonstrated that C-modes are phase-locked within GS-cycles by sampling the GS-phase difference between the occurrence of C-modes within a GS-cycle, and the following occurrences of another C-mode within the same cycle or the immediately subsequent one, limiting sampling to instances in which C-mode time-courses exceeded 1 standard deviation (Fig. S8).

To further explore the temporal sequencing of C-mode occurrences and transitions (Fig. 5), we defined, for each dataset, a concatenated sequence of C-mode occurrences across subjects/animals, and calculated the transition probability matrix from a C-mode i at time t to another C-mode j at time $t+1$ as the proportion of transitions $i\to j$ and all other transitions from i^19,29. Only transitions within the same subject were included, and we first considered sequences in which we counted the auto-transitions ($i=j$), namely persistence probabilities. Off-diagonal elements ($i\ne j$), named transition probabilities, were then computed after building sequences in which we removed the repeating elements in order to control for auto-correlations given the C-mode’s dwell time^19,86. We quantified the directional prevalence of a transition ($P_{ij}>P_{ji}$) by taking their difference. To measure the relative facilitation to reach a C-mode from another, we computed the Entropy of Markov Trajectories (HMT)^29,32 from the off-diagonal matrix elements (i.e. transition probabilities). This method calculates the descriptive complexity of the paths between C-mode (in bits), where a lower complexity refers to less information required to access a destination C-mode j from a source C-mode i, hence being more accessible as the path travels through less path before reaching its destination. Specifically, for a Markov Chain (MC) defined by transition probability matrix P, we define the Entropy Rate per step $t\to t+1$ as:

Where $\mu$ is the stationary distribution solving ${\mu }_j={\sum }_i{\mu }_i{P}_{ij}.$ Also define the Entropy of a Markov Trajectory $T_{ij}$ from C-mode i to C-mode j as:

Where $\mathbf{K}={\left(\mathbf{I}-\mathbf{P}-\mathbf{M}\right)}^{-\mathbf{1}}\left({\mathbf{H}}^{*}-{\mathbf{H}}_{\mathbf{\Delta }}\right)$, M is a matrix of stationary probabilities $M_{ij}={\mu }_{ij}$; I the identity matrix; ${\mathbf{H}}^{*}$ is the matrix of single-step entropies $H_{ij}^{*}=H\left(P_i\right)={\sum }_k{\mu }_k{P}_{ik}$ (from C-mode i to any C-mode k); and ${\mathbf{H}}_{\mathbf{\Delta }}$ is a diagonal entropy matrix with trajectories from a C-mode to itself ${\left(H_{\mathrm{\Delta }}\right)}_{ii}=H\left(\mathrm{MC}\right)/{\mu }_i$, with zeros if $i\ne j$³².

Statistical significance of persistence probabilities was tested from random sequences by generating 1000 permutations of C-mode occurrence sequences at the subject/animal level, concatenating the sequences, and counting how many times a random iteration exceeded the true value to reach a p value. Transition probabilities and directional transition prevalence were tested against matrices built after randomly permuting the non-repeating sequences 1000 times at the subject/animal level, concatenating them and removing samples if ${C\mathrm{mode}}_i\left(t+1\right)={C\mathrm{mode}}_i\left(t\right)$. These surrogates were also used to test the significance of Markov trajectory entropies. To further quantify and assess the relative complexity of transitions to a particular destination C-mode, we quantified, after repeating the analyses with subject level sequences, the sum of the HMT matrix columns, and compared their distributions across C-mode destinations (one-way ANOVA and Tukey test for multiple comparisons, Fig. 5C). Finally, to assess the relationship of the mean GS-phase difference between occurrences of different C-modes and the HMT of their transitions, we computed the circular-linear correlations of their mean values (Fig. 5B and S8, respectively) using the circ_corrcl.m Matlab function⁸⁵.

C-mode influence on static fMRI connectome and functional connectivity gradients

We first generated C-mode co-fluctuation matrices by cross-multiplying the mean C-mode map with itself⁴², which yielded a voxel-wise representation of the co-activating (or co-deactivating) peaks of fMRI activity. We then computed the group-level mean Functional Connectivity matrices by computing the Pearson’s Correlation between the concatenated time-series of voxels in all subjects/animals in each group⁸⁷, and computed the correlation between the vectorized upper triangular part of the FC matrix to the weighted (by occurrence rate) average of the C-mode Co-fluctuation matrices²⁹ (Fig. S9A, B). Voxels were ordered according to the RSN they belong to (see Fig. 1C).

We further explored if the hierarchies that dominate FC are in accordance with the dynamic structure of C-modes. We first computed the gradients of FC matrices from each species using the Diffusion Mapping method (BrainSpace⁸⁸ - diffusion_mapping.m Matlab function: FC matrix sparse to top 10-percentile per node, 5 components, anisotropic diffusion parameter = 0.5)^69,89. We mapped the gradients (Figure S10) and computed their spatial correlation to the obtained C-modes, matching each C-mode to the gradient with the highest absolute Pearson’s correlation (Fig. 6A). Finally, we plotted the variance explained by each gradient (lambda) to the group mean occurrence rate of their matched C-mode, fitting a linear model to this relationship (Fig. 6B).

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.

Peer review information

Nature Communications thanks Joana Cabral and Sung-Ho Lee for their contribution to the peer review of this work. A peer review file is available.