1. Neuroscience
Download icon

Topographic gradients of intrinsic dynamics across neocortex

  1. Golia Shafiei  Is a corresponding author
  2. Ross D Markello
  3. Reinder Vos de Wael
  4. Boris C Bernhardt
  5. Ben D Fulcher
  6. Bratislav Misic  Is a corresponding author
  1. McConnell Brain Imaging Centre, Montréal Neurological Institute, McGill University, Canada
  2. School of Physics, The University of Sydney, Australia
Research Article
  • Cited 7
  • Views 1,950
  • Annotations
Cite this article as: eLife 2020;9:e62116 doi: 10.7554/eLife.62116

Abstract

The intrinsic dynamics of neuronal populations are shaped by both microscale attributes and macroscale connectome architecture. Here we comprehensively characterize the rich temporal patterns of neural activity throughout the human brain. Applying massive temporal feature extraction to regional haemodynamic activity, we systematically estimate over 6000 statistical properties of individual brain regions’ time-series across the neocortex. We identify two robust spatial gradients of intrinsic dynamics, one spanning a ventromedial-dorsolateral axis and dominated by measures of signal autocorrelation, and the other spanning a unimodal-transmodal axis and dominated by measures of dynamic range. These gradients reflect spatial patterns of gene expression, intracortical myelin and cortical thickness, as well as structural and functional network embedding. Importantly, these gradients are correlated with patterns of meta-analytic functional activation, differentiating cognitive versus affective processing and sensory versus higher-order cognitive processing. Altogether, these findings demonstrate a link between microscale and macroscale architecture, intrinsic dynamics, and cognition.

Introduction

The brain is a complex network of anatomically connected and perpetually interacting neuronal populations (Sporns et al., 2005). Inter-regional connectivity promotes signaling via electrical impulses, generating patterned electrophysiological and haemodynamic activity (Avena-Koenigsberger et al., 2017; Suárez et al., 2020). Neuronal populations are organized into a hierarchy of increasingly polyfunctional neural circuits (Jones and Powell, 1970; Mesulam, 1998; Hilgetag and Goulas, 2020; Bazinet et al., 2020), manifesting as topographic gradients of molecular and cellular properties that smoothly vary between unimodal and transmodal cortices (Huntenburg et al., 2018). Recent studies have demonstrated cortical gradients of gene transcription (Fulcher et al., 2019; Burt et al., 2018), intracortical myelin (Huntenburg et al., 2017), cortical thickness (Wagstyl et al., 2015) and laminar profiles (Paquola et al., 2019).

The topological and physical embedding of neural circuits in macroscale networks and microscale gradients influence their dynamics (Kiebel et al., 2008; Gollo et al., 2015; Wang, 2020). For a neuronal population, the confluence of local micro-architectural properties and global connectivity shapes both the generation of local rhythms, as well as its propensity to communicate with other populations. Specifically, cell type composition, their morphology and their configuration in local circuits determine how signals are generated, transmitted and integrated (Payeur et al., 2019). These micro-architectural properties – increasingly measured directly from histology or inferred from other measurements, such as microarray gene expression – provide a unique opportunity to relate circuit architecture to temporal dynamics and computation. Indeed, multiple studies have focused on how intrinsic timescales vary in relation to microscale and macroscale attributes (Murray et al., 2014; Mahjoory et al., 2019; Shine et al., 2019; Gao et al., 2020; Ito et al., 2020; Raut et al., 2020). The primary functional consequence of this hierarchy of timescales is thought to be a hierarchy of temporal receptive windows: time windows in which a newly arriving stimulus will modify processing of previously presented (i.e. contextual) information (Hasson et al., 2008; Honey et al., 2012; Baldassano et al., 2017; Huntenburg et al., 2018; Chaudhuri et al., 2015; Chien and Honey, 2020). Thus, areas at the bottom of the hierarchy preferentially respond to immediate changes in the sensory environment, while responses in areas at the top of the hierarchy are modulated by prior context. Altogether, previous work highlights a hierarchy of a small number of manually selected time-series features, but it is possible that different types of local computations manifest as different organizational gradients.

The relationship between structure and dynamics is also observed at the network level (Suárez et al., 2020). Intrinsic or ‘resting state’ networks possess unique spectral fingerprints (Keitel and Gross, 2016). The signal variability of brain areas, measured in terms of standard deviations or temporal entropy, is closely related to their structural and functional connectivity profiles (i.e. network embedding) (Mišić et al., 2011; Burzynska et al., 2013; Garrett et al., 2017; Shafiei et al., 2019). More generally, the autocorrelation of blood oxygenation level-dependent (BOLD) signal is correlated with topological characteristics of structural brain networks, such that areas with greater connectivity generate signals with greater autocorrelation (Sethi et al., 2017; Fallon et al., 2020). Finally, in computational models of structurally coupled neuronal populations (neural mass and neural field models [Breakspear, 2017]), highly interconnected hubs exhibit slower dynamic fluctuations, while sensory areas exhibit fast fluctuating neural activity (Gollo et al., 2015). Indeed, these models offer better fits to empirical functional connectivity if they assume heterogeneous local dynamics (Cocchi et al., 2016; Demirtaş et al., 2019; Wang et al., 2019; Deco et al., 2020).

Altogether, multiple lines of evidence suggest that local computations may reflect systematic variation in microscale properties and macroscale network embedding, manifesting as diverse time-series features of regional neural activity. How molecular, cellular and connectomic architecture precisely shapes temporal dynamics, and ultimately, cortical patterns of functional specialization, is poorly understood. A significant limitation is that conventional computational analysis is based on specific, manually selected time-series features, such as the decay of the autocorrelation function, bands of the Fourier power spectrum, or signal variance. Yet the time-series analysis literature is vast and interdisciplinary; how do other metrics of temporal structure vary across the brain and what can they tell us about cortical organization?

Here we comprehensively chart summary features of spontaneous BOLD signals across the cerebral cortex (hereafter referred to as ‘intrinsic dynamics’), mapping temporal organization to structural organization. We apply massive temporal feature extraction to resting state BOLD signals to derive a near-exhaustive time-series profile for each brain region. We then systematically investigate the relationship between local time-series features and gene expression, microstructure, morphology, structural connectivity and functional connectivity. Finally, we map time-series features to a meta-analytic atlas of cognitive ontologies to investigate how temporal dynamics shape regional functional specialization. We show that intrinsic dynamics reflect molecular and cytoarchitectonic gradients, as well as patterns of structural and functional connectivity. These spatial variations in intrinsic dynamics ultimately manifest as patterns of distinct psychological functions.

Results

All analyses were performed on four resting state fMRI runs from the Human Connectome Project (Van Essen et al., 2013). The data were pseudorandomly divided into two samples of unrelated participants to form Discovery and Validation samples with n=201 and n=127, respectively (Vos de Wael et al., 2018). External replication was then performed using data from the Midnight Scan Club (Gordon et al., 2017). Massive temporal feature extraction was performed using highly comparative time-series analysis, hctsa (Fulcher et al., 2013; Fulcher and Jones, 2017), yielding 6441 features per regional time-series, including measures of frequency composition, variance, autocorrelation, fractal scaling and entropy (Figure 1). The results are organized as follows. We first investigate whether regions that are structurally and functionally connected display similar intrinsic dynamics. We then characterize the topographic organization of time-series features in relation to microstructural attributes and cognitive ontologies.

Temporal phenotyping of regional dynamics.

The highly comparative time-series analysis toolbox, hctsa (Fulcher et al., 2013; Fulcher and Jones, 2017), was used to extract 6441 time-series features of the parcellated time-series for each brain region and participant, including measures of autocorrelation, variance, spectral power, entropy, etc. Regional time-series profiles were then entered into two types of analyses. In the first analysis, pairs of regional time-series feature vectors were correlated to generate a region × region temporal profile similarity network. In the second analysis, principal component analysis (PCA) was performed to identify orthogonal linear combinations of time-series features that vary maximally across the cortex.

Inter-regional temporal profile similarity reflects network geometry and topology

We first assessed the extent to which intrinsic dynamics depend on inter-regional physical distance, anatomical connectivity and functional connectivity. We estimated similarity between inter-regional dynamics by computing Pearson correlation coefficients between regional time-series feature vectors (Figure 1). Two regional time-series are judged to be similar if they have similar temporal profiles, estimated across a comprehensive and diverse set of time-series features (e.g. similar entropy, stationarity, linear correlation properties) (Fulcher, 2018). This measure of similarity identifies pairs of regions that have similar dynamical features, but not necessarily coherent or synchronous dynamics (Figure 2a). We refer to correlations between regional time-series feature profiles as ‘temporal profile similarity’.

Figure 2 with 2 supplements see all
Inter-regional temporal profile similarity reflects network geometry and topology.

(a) Temporal profile similarity networks are constructed by correlating pairs of regional time-series feature vectors. Brain regions are ordered based on their intrinsic functional network assignments (Yeo et al., 2011; Schaefer et al., 2018). (b) Temporal profile similarity between regions significantly decreases as a function of Euclidean distance between them. The black line represents an exponential fit as y=0.37e-0.03x+0.01, where y is temporal profile similarity and x is Euclidean distance. (c, d) Regional time-series features are compared between pairs of cortical areas using their structural and functional connectivity profiles. Pairwise temporal profile similarity is significantly higher among structurally-connected areas (c), and among regions that belong to the same intrinsic functional networks (d). Asterisks denote a statistically significant difference of the means (two-tailed t-test; p0). For structural networks, statistical significance of the difference of the mean temporal profile similarity of connected and unconnected node pairs is also assessed against a null distribution of differences generated from a population degree- and edge length-preserving rewired networks (Betzel and Bassett, 2018) (c, right-most panel). For functional networks, statistical significance of the difference of the mean temporal profile similarity of within and between intrinsic networks is also assessed against a null distribution of differences generated by spatial autocorrelation-preserving label permutation (‘spin tests’; Alexander-Bloch et al., 2018) (d, right-most panel). (e) Temporal profile similarity is positively correlated with functional connectivity. This relationship remains after partialling out Euclidean distance between regions from both measures using exponential trends. rs denotes the Spearman rank correlation coefficient; linear regression lines are added to the scatter plots for visualization purposes only. Connections are color-coded based on the intrinsic network assignments (Yeo et al., 2011; Schaefer et al., 2018). VIS = visual, SM = somatomotor, DA = dorsal attention, VA = ventral attention, LIM = limbic, FP = fronto parietal, DMN = default mode.

Figure 2b shows a negative exponential relationship between spatial proximity and temporal profile similarity, meaning that regions that are spatially close exhibit similar intrinsic dynamics. Interestingly, regions that share an anatomical projection have greater temporal profile similarity than those that do not (Figure 2c; two-tailed t-test; t(79,798)=40.234, p0). To test whether this anatomically-mediated similarity of time-series features is not due to spatial proximity, we performed two additional comparisons. First, we regressed out the exponential trend identified above from the temporal profile similarity matrix, and repeated the analysis on the residuals, yielding a significant difference in temporal profile similarity between connected and non-connected regions (two-tailed t-test; t(79,798)=9.916, p0). Second, we generated an ensemble of 10,000 degree- and edge length-preserving surrogate networks (Betzel and Bassett, 2018), and compared the difference of the means between connected and non-connected pairs in the empirical and surrogate networks. Again, we observe a significant difference in temporal profile similarity between connected and non-connected regions (two-tailed; prewired=0.0001; Figure 2c).

Likewise, regions belonging to the same intrinsic functional network have greater temporal profile similarity compared to regions in different networks (Figure 2d; two-tailed t-test; t(79,798)=61.093, p0). To confirm this finding is not driven by spatial proximity, we repeated the analysis with distance-residualized values (Mišić et al., 2014), finding a significant difference (two-tailed t-test; t(79,798)=47.112, p0). We also repeated the analysis using a nonparametric label-permutation null model with preserved spatial autocorrelation (10,000 repetitions) (Alexander-Bloch et al., 2018; Markello and Misic, 2020), again finding significantly greater within- compared to between-network temporal profile similarity (two-tailed; pspin=0.006; Figure 2d). These results are consistent when applying the 17 network partition of intrinsic networks (Yeo et al., 2011; Schaefer et al., 2018; Figure 2—figure supplement 1).

More generally, we find a weak positive correlation between temporal profile similarity and functional connectivity (original: Spearman rank rs=0.23, p0; distance-corrected: rs=0.18, p0; Figure 2e), suggesting that areas with similar time-series features exhibit coherent spontaneous fluctuations, but that the two are only weakly correlated. Figure 2e shows the correlation between temporal profile similarity and functional connectivity; points represent node pairs and are colored by their membership in intrinsic networks (Yeo et al., 2011; Schaefer et al., 2018). The results are consistent when functional connectivity is estimated using partial correlations (Figure 2—figure supplement 2). In other words, two regions could display similar time-series features, but they do not necessarily fluctuate coherently. Thus, representing time-series using sets of features provides a fundamentally different perspective compared to representing them as the raw set of ordered BOLD measurements.

As a final step, we sought to assess the distinct contributions of Euclidean distance, structural connectivity and functional connectivity to temporal profile similarity. Dominance analysis revealed the relative importance of each predictor (collective R2 = 0.28; distance = 56%, structural connectivity = 20.4%, functional connectivity = 23.6%; Supplementary file 1), suggesting that distance contributes the most to temporal profile similarity, while structural and functional connectivity make distinct but approximately even contributions (Budescu, 1993; Azen and Budescu, 2003) (https://github.com/dominance-analysis/dominance-analysis). Altogether, we find that the organization of intrinsic dynamics is closely related to both the geometric and topological embedding of brain regions in macroscale networks.

Two distinct spatial gradients of intrinsic dynamics

We next investigate the topographic organization of time-series features. The hctsa library generates 6 441 time-series features, with the aim of being comprehensive in coverage across scientific time-series analysis algorithms and, as a result, contains groups of features with correlated outputs (Fulcher et al., 2013). We therefore sought to identify groups of correlated features that explain maximal variance and that span different conceptual types of time-series properties. Applying principal component analysis (PCA; Scikit-learn [Pedregosa et al., 2011]) to the region × feature matrix yielded mutually orthogonal patterns of intrinsic dynamics (Figure 1), with the top two components collectively accounting for more than 70% of the variance in time-series features (Figure 3a). Figure 3a shows the spatial distribution of the top two components. The first component (PC1) mainly captures differential intrinsic dynamics along a ventromedial-dorsolateral gradient, separating occipital-parietal cortex and anterior temporal cortex. The second component (PC2) captures a unimodal-transmodal gradient, reminiscent of recently reported miscrostructural and functional gradients (Huntenburg et al., 2018). Both components show considerable hemispheric symmetry. In the following sections, we focus on these two components because of their (a) effect size (percent variance accounted for), (b) close resemblance to previously reported topographic gradients, and (c) reproducibility (only the first two components were reproducible in both the HCP and MSC datasets; see Sensitivity and replication analyses below). Note that neither spatial maps were significantly correlated with temporal signal-to-noise ratio map, computed as the ratio of the time-series mean to standard deviation (tSNR; PC1: rs=0.28, pspin=0.19; PC2: rs=0.21, pspin=0.16).

Figure 3 with 1 supplement see all
Topographic gradients of intrinsic dynamics.

(a) PCA analysis identified linear combinations of hctsa time-series features with maximum variance across the cortex. Collectively, the first two components (PC1 and PC2) account for 75% of the total variance in time-series features of BOLD dynamics. To estimate the extent to which cortical regions display the patterns of intrinsic dynamics captured by each component, hctsa matrices were projected back onto the PC weights (eigenvectors), yielding spatial maps of brain scores for each component. Spatial maps are depicted based on the standard deviation σ of their respective brain score distributions. (b) To understand the feature composition of the intrinsic dynamic patterns captured by PC1 and PC2, feature loadings were computed by correlating individual hctsa feature vectors with the PC score maps. PC loadings thus estimate the shared spatial variance between an individual time-series feature and the composite intrinsic dynamic map captured by a PC. time-series features are ordered by their individual loadings. Grey indicates non-significance based on 10,000 spatial permutation tests (FDR correction, α=0.001). Features corresponding to the top and bottom 5% of PC1 and PC2 are visualized using word clouds. The complete list of features (ranked by loading), their definitions, correlations and p-values for both components is presented in machine-readable format in Supplementary files 3 and 4. Feature nomenclature in hctsa is organized such that the term prefix indicates the broad class of measures (e.g. AC = autocorrelation, DN = distribution) and the term suffix indicates the specific measure (for a complete list, see https://hctsa-users.gitbook.io/hctsa-manual/list-of-included-code-files). (c) The spatial distributions of two high-loading representative time-series features are depicted for each component, including lag-1 linear autocorrelation (AC_1) and lag-[0,2,3] nonlinear autocorrelation (AC_nl_023, estimated as average <xt2xt-2xt-3> across time-series x) for PC1; and kurtosis (DN_Moments_4) and entropy (EN_DistributionEntropy_ks__02) of the time-series points distribution for PC2. To build intuition about what each component reflects about regional signals, three regional time-series from one participant are selected based on their lag-1 autocorrelation and kurtosis (circles on the brain surface: pink = 5th percentile, green = 50th percentile, purple = 95th percentile). Going from low-ranked to high-ranked regions results in a slowing down of BOLD fluctuations for the former and increasingly heavier symmetric tails of the signal amplitude distributions for the latter. Note that normalized feature values are shown in the first row, whereas the raw feature values are shown in the second row.

Which time-series features contribute most to these topographic gradients of intrinsic dynamics? To address this question, we systematically assess the feature composition of PC1 and PC2. We compute univariate correlations (i.e. loadings) between individual time-series feature vectors and PC scores (Figure 3b). Each loading is assessed against 10 000 spin tests and the results are corrected for multiple comparisons by controlling the false discovery rate (FDR [Benjamini and Hochberg, 1995]; α=0.001). The top 5% positively and negatively correlated features are shown in word clouds. The complete list of features (ranked by loading), their definitions, loadings and p-values for both components is presented in machine-readable format in Supplementary files 3 and 4. Altogether, we find that PC1 is sensitive to temporal dependencies in BOLD signals, while PC2 is sensitive to the distribution shape of time-series amplitudes. For PC1, in line with previous reports, we observe strong contributions from multiple measures of autocorrelation (e.g. linear autocorrelation; nonlinear autocorrelation; automutual information). Short-lag autocorrelation measures load positively, while long-lag autocorrelation measures load negatively, consistent with the notion that autocorrelation decays with increasing time lag (Murray et al., 2014; Gao et al., 2020; Raut et al., 2020; Figure 3—figure supplement 1). For PC2, we observe strong contributions from measures of distribution shape, captured by measures of distributional entropy (e.g. entropy of kernel-smoothed distribution; kurtosis; distribution balance about the mean). In other words, PC2 captures the spread of time-series amplitudes away from the mean. Interestingly, none of the odd moments (distribution asymmetry) are high in the PC2 loading list, just even moments, suggesting that PC2 captures the shape of the deviations of time-series data points in both directions from the mean. Thus, PC2 indexes the range or diversity of values that a regional time-series can realize. Hereafter, we refer to the time-series profile of PC1 as ‘autocorrelation’ and PC2 as ‘dynamic range’.

To illustrate the spatial organization and time-series attributes of these components, Figure 3c shows the spatial distributions of two high-loading representative time-series features for each component. Ventromedial areas (lower in the PC1 gradient) have lower linear and nonlinear autocorrelation, while doroslateral areas (higher in the PC1 gradient) have greater autocorrelation. Sensory areas (lower in the PC2 gradient) have greater distributional entropy and kurtosis, while transmodal areas (higher in the PC2 gradient) have lower distributional entropy and kurtosis. Finally, to build intuition about what each component reflects about regional signals, we select three regional time-series from one participant based on their lag-1 autocorrelation and kurtosis (Figure 3c; pink = 5th percentile, green = 50th percentile, purple = 95th percentile). For the former, going from low-ranked to high-ranked regions results in a slowing down of BOLD fluctuations. For the latter, going from low-ranked to high-ranked regions results in increasingly heavier symmetric tails of the signal amplitude distributions.

Intrinsic dynamics reflect microscale and macroscale hierarchies

To assess whether the dominant variation in time-series features of BOLD dynamics varies spatially with structural and functional gradients, we next quantify the concordance between PC1/PC2 and multiple microstructural and functional attributes (Figure 4). Specifically, we compare PC1 and PC2 with the following microscale and macroscale features: (1) the first component of microarray gene expression computed from the Allen Institute Human Brain Atlas (Hawrylycz et al., 2012; Burt et al., 2018) using PCA analysis, (2) the principal gradient of functional connectivity estimated using diffusion map embedding (Margulies et al., 2016; Coifman et al., 2005; Langs et al., 2015) (https://github.com/satra/mapalign), (3) T1w/T2w ratio, a putative proxy for intracortical myelin (Huntenburg et al., 2017), (4) cortical thickness (Wagstyl et al., 2015). We use Spearman rank correlations (rs) throughout, as they do not assume a linear relationship among variables. Given the spatially autocorrelated nature of both hctsa features and other imaging features, we assess statistical significance with respect to nonparametric spatial autocorrelation-preserving null models (Alexander-Bloch et al., 2018; Markello and Misic, 2020).

Figure 4 with 2 supplements see all
Hierarchical organization of intrinsic dynamics.

PC1 and PC2 brain score patterns are compared with four molecular, microstructural and functional maps. These maps include the first principal component of microarray gene expression data from the Allen Human Brain Atlas (Hawrylycz et al., 2012; Burt et al., 2018), the first (principal) gradient of functional connectivity estimated using diffusion map embedding (Margulies et al., 2016; Coifman et al., 2005; Langs et al., 2015), group-average T1w/T2w ratio, and group-average cortical thickness. The three latter indices were computed from the HCP dataset (Van Essen et al., 2013). Statistical significance of the reported Spearman rank correlation rs is assessed using 10,000 spatial permutations tests, preserving the spatial autocorrelation in the data (‘spin tests’; Alexander-Bloch et al., 2018). Linear regression lines are added to the scatter plots for visualization purposes only.

PC1 topography is correlated with the first principal component of gene expression (rs=0.57, pspin=0.03), but no other attributes. PC2 topography is significantly correlated with the first principal component of gene expression (rs=-0.45, pspin=0.0008), with the principal gradient of functional connectivity (rs=0.77, pspin=0.0001), with T1w/T2w ratio (rs=-0.57, pspin=0.0001), and with cortical thickness (rs=0.43, pspin=0.008). Altogether, the two topographic gradients of intrinsic dynamics closely mirror molecular and microstructural gradients, suggesting a link between regional structural properties and regional dynamical properties. Figure 4—figure supplement 1 further confirms this intuition, showing the mean score of each component for three well-known cortical partitions, including intrinsic functional networks (Yeo et al., 2011; Schaefer et al., 2018), cytoarchitectonic classes (von Economo and Koskinas, 1925; von Economo et al., 2008; Vértes et al., 2016) and laminar differentiation levels (Mesulam, 2000).

For completeness, we also tested associations with two maps that were previously related to cortical hierarchies: evolutionary expansion (indexing enlargement of cortical areas in the human relative to the macaque) (Hill et al., 2010; Baum et al., 2020) and node-wise functional participation coefficient (indexing the diversity of a node’s links) (Bertolero et al., 2017; Baum et al., 2020). PC2 is significantly correlated with evolutionary expansion (rs=0.52, pspin=0.0002), but neither component is correlated with participation coefficient (Figure 4—figure supplement 2).

Spatial gradients of intrinsic dynamics support distinct functional activations

Given that topographic patterns of intrinsic dynamics run parallel to microstructural and functional gradients, and are marked by specific time-series features, we next asked whether these topographic patterns of intrinsic dynamics are related to patterns of functional activation and psychological processes. To address this question, we used Neurosynth to derive probability maps for multiple psychological terms (Yarkoni et al., 2011). The term set was restricted to those in the intersection of terms reported in Neurosynth and in the Cognitive Atlas (Poldrack et al., 2011), yielding a total of 123 terms (Supplementary file 2). Each term map was correlated with the PC1 and PC2 score maps to identify topographic distributions of psychological terms that most closely correspond to patterns of intrinsic dynamics (Bonferroni corrected, α=0.05; Figure 5). Consistent with the intuition developed from comparisons with intrinsic networks, PC1 intrinsic dynamics mainly defined a cognitive-affective axis (e.g. ‘attention’ versus ‘stress’, ‘fear’, ‘loss’, ‘emotion’; Figure 5a), while PC2 dynamics defined a sensory-cognitive axis (e.g. ‘perception’, ‘multisensory’, ‘facial expression’ versus ‘cognitive control’, ‘memory retrieval’, ‘reasoning’; Figure 5b).

Spatial gradients of intrinsic dynamics support distinct functional activations.

We used Neurosynth to derive probability maps for multiple psychological terms (Yarkoni et al., 2011). The term set was restricted to those in the intersection of terms reported in Neurosynth and in the Cognitive Atlas (Poldrack et al., 2011), yielding a total of 123 terms (Supplementary file 2). Each term map was correlated with the PC1 (a) and PC2 (b) score maps to identify topographic distributions of psychological terms that most closely correspond to patterns of intrinsic dynamics. Grey indicates non-significance based on 10,000 spatial permutation tests (Bonferroni correction, α=0.05). Statistically significant terms are shown on the right.

Sensitivity and replication analyses

As a final step, we sought to assess the extent to which the present findings are replicable under alternative processing choices and in other samples (Figure 6). For all comparisons, we correlated PC1 and PC2 scores and weights obtained in the original analysis and in each new analysis. Significance was assessed using spatial autocorrelation preserving nulls as before. We first replicated the results in individual subjects in the Discovery sample by applying PCA to individual region × feature matrices and aligning PCA results through an iterative process using Procrustes rotations (https://github.com/satra/mapalign [Langs et al., 2015]). The mean individual-level PC scores and weights were then compared to the original findings (Figure 6a). We next replicated the results by repeating the analysis after grey-matter signal regression (similar to global signal regression as the global signal is shown to be a grey-matter specific signal following sICA+FIX) (Glasser et al., 2018; Glasser et al., 2016), with near identical results (Figure 6b). To assess the extent to which results are influenced by choice of parcellation, we repeated the analysis using the 68-region Desikan-Killiany anatomical atlas (Desikan et al., 2006), which were then further divided into 200 approximately equally-sized cortical areas. Again, we find near-identical results (Figure 6c).

Sensitivity and replication analyses.

For all comparisons, we correlated PC1 and PC2 scores and weights obtained in the original analysis and in each new analysis. Significance was assessed using spatial autocorrelation preserving nulls. Specific analyses include: (a) comparing group-level and individual subject-level results, (b) comparing data with and without grey-matter signal regression, (c) comparing functional (Schaefer) and anatomical parcellations (Desikan-Killiany), (d) comparing HCP Discovery and Validation datasets, (e) comparing HCP Discovery and Midnight Scan Club datasets.

In the last two analyses, we focused on out-of-sample validation. We first repeated the analysis on the held-out Validation sample of n=127 unrelated HCP subjects, with similar results (Figure 6d). Finally, we repeated the analysis using data from the independently collected Midnight Scan Club (MSC) dataset, again finding highly consistent results (Figure 6e).

Discussion

In the present report, we comprehensively characterize intrinsic dynamics across the cortex, identifying two robust spatial patterns of time-series features. The patterns, capturing spatial variation in signal autocorrelation and dynamic range, follow microscale gradients and macroscale network architecture. Importantly, the two patterns underlie distinct psychological axes, demonstrating a link between brain architecture, intrinsic dynamics, and cognition. These findings are robust against a wide range of methodological choices and were validated in two held-out samples.

Our results demonstrate that regional haemodynamic activity, often overlooked in favour of electrophysiological measurements with greater temporal resolution, possesses a rich dynamic signature (Garrett et al., 2013; Uddin, 2020; Preti et al., 2017; Lurie et al., 2020; Li et al., 2019; Bolt et al., 2018). While multiple reports have suggested the existence of a timescale or temporal receptive window hierarchy (Kiebel et al., 2008; Murray et al., 2014; Honey et al., 2012; Hasson et al., 2008; Watanabe et al., 2019; Golesorkhi et al., 2020; Ito et al., 2020), these investigations typically involved (a) incomplete spatial coverage, making it difficult to quantitatively assess correspondence with other microscale and macroscale maps, and (b) a priori measures of interest, such as spectral power or temporal autocorrelation, potentially obscuring other important dynamical features. Here we comprehensively benchmark the entire dynamic profile of the brain, by near-exhaustively estimating 6000+ features from the wider time-series literature. We identify a much broader spectrum of time-series features that relate to microstructure, connectivity and behavior. As we discuss below, feature-based time-series phenotyping offers a powerful, fundamentally new and entirely data-driven method to quantify and articulate neural dynamics.

Applying a data-driven feature extraction method to high-resolution BOLD fMRI, we decompose regional signals into two intrinsic modes, with distinct topographic organization and time-series features. One pattern, characterized by variation in signal autocorrelation, follows a ventromedial-dorsolateral gradient, separating the limbic and paralimibic systems from posterior parietal cortex. Another pattern, characterized by dynamic range, follows a unimodal-transmodal gradient, separating primary sensory-motor cortices from association cortex. The first is closely associated with gene expression PC1 (itself closely related to cell type composition, synaptic physiology and cortical cytoarchitecture [Burt et al., 2018]), suggesting a molecular and cellular basis for regional differences in temporal autocorrelation. The second is closely associated with the principal functional gradient, as well as with intracortical myelin and cortical thickness, suggesting that the dynamic range of BOLD signals is related to regional variation in macroscale circuit organization. Taken together, we find evidence that molecular and cellular properties (gene expression PC1) relate to regional autocorrelation, while micro-circuit properties (T1w/T2w, cortical thickness) and macroscale network embedding (principal functional gradient) relate to regional dynamic range.

An emerging literature emphasizes the hierarchical organization of neural systems, whereby systematic variation in laminar architecture across the cortical sheet is mirrored by multiple cytological properties, including neuron density, spine count, branching and neurotransmitter receptor profiles (Mesulam, 1998; Margulies et al., 2016; Hilgetag and Goulas, 2020). These variations manifest as spatially ordered gradients of structural and functional attributes (Huntenburg et al., 2018), including gene expression (Burt et al., 2018; Fulcher et al., 2019; Hansen et al., 2020), cortical thickness (Wagstyl et al., 2015), intracortical myelin (Huntenburg et al., 2017), laminar differentiation (Paquola et al., 2019; Wagstyl et al., 2020) and excitability (Demirtaş et al., 2019; Wang, 2020; Markicevic et al., 2020; Straub et al., 2020). Indeed, we find that the two patterns of intrinsic dynamics are closely related to gene expression, intracortical myelin and cortical thickness. Our results build on this literature, demonstrating that microscale and connectional hierarchies leave an indelible mark on intrinsic dynamics (Lurie and D’Esposito, 2020), perhaps through variation in local excitability (Demirtaş et al., 2019; Wang et al., 2019; Wang, 2020; Deco et al., 2020). How these patterns are related to underlying cell types and subcortical afferent input – in particular, thalamocortical feedback – is an important ongoing question (Abeysuriya et al., 2015; Garrett et al., 2018; Shine et al., 2019; Wang et al., 2019; Muller et al., 2020; Paquola et al., 2020).

Importantly, the two patterns are related to two dominant axes of meta-analytic functional activation. We show that topographic variations in microcircuitry and connectomic embedding yield variations in intrinsic dynamics and may explain regional differences in functional specialization. The ventromedial-dorsolateral autocorrelation pattern differentiates affective versus cognitive activation (mainly visual cognition and visuo-spatial attention), whereas the unimodal-transmodal dynamic range pattern differentiates primary sensory versus higher-order cognitive processing. Collectively, these results provide evidence that local computations reflect systematic variation in multiple anatomical circuit properties, and can be measured as unique temporal signatures in regional activity and patterns of functional specialization.

More generally, the present findings are part of a larger trend in the field to understand structure-function relationships by considering molecular (Richiardi et al., 2015; Fulcher and Fornito, 2016; Anderson et al., 2018; Zheng et al., 2019), cellular (Scholtens et al., 2014; Anderson et al., 2020; Shafiei et al., 2020; Muller et al., 2020) and physiological (Sethi et al., 2017; Fallon et al., 2020) attributes of network nodes, thereby conceptually linking local and global brain organization (Khambhati et al., 2018; Suárez et al., 2020). In such ‘annotated networks’, macroscale network architecture is thought to reflect similarity in local properties, and vice versa, such that areas with similar properties are more likely to be anatomically connected and to functionally interact with one another (Beul et al., 2017; Goulas et al., 2019; Wei et al., 2019; Hilgetag et al., 2019). Indeed, we find that two regions are more likely to display similar intrinsic dynamics if they are anatomically connected and if they are part of the same functional community, suggesting that network organization and local intrinsic dynamics are intertwined (Gollo et al., 2015; Cocchi et al., 2016). A significant corollary of the present work is that functional connectivity – presently conceptualized as coherent fluctuations in neural activity and operationalized as correlated BOLD values over time – misses out on an important set of inter-regional relationships. Namely, two regions may display identical time-series profiles, suggesting common circuit dynamics and function, but unless they also display time-locked activity, current methods would miss out on this potentially biologically meaningful inter-regional relationship.

The present results are consistent with contemporary theories linking brain structure and function, but they must be interpreted with respect to several methodological caveats. First, all analyses were performed on BOLD time-series with lower sampling rate compared to electromagnetic recordings, potentially obscuring more subtle dynamics occurring on faster timescales. To mitigate this concern, all analyses were performed in high-resolution multiband HCP data with multiple runs, and replicated in MSC data, but in principle, these analyses could be repeated and validated in magnetoencephalographic recordings (Watanabe et al., 2019). Second, all analyses were performed on haemodynamic time courses that may not completely reflect the underlying neuronal population dynamics. Despite this caveat, we observe a close correspondence between the isolated patterns of intrinsic dynamics and molecular, structural, functional, and psychological gradients. Third, the pattern of temporal signal-to-noise ratio in the BOLD is known to be non-uniform, but it is not correlated with the intrinsic dynamics patterns observed in the present report. Fourth, the analysis included all features from hctsa, potentially biasing results towards specific properties of BOLD signals. We attempted to mitigate this challenge by applying PCA to directly examine correlation patterns among features, but PCA components may still lend greater weight to over-represented feature classes (Fulcher et al., 2013). This may obscure the contribution of under-represented feature classes, and should be investigated further in future work.

Altogether, the present results point towards highly patterned intrinsic dynamics across the neocortex. These patterns reflect prominent molecular and microstructural gradients, as well as macroscale structural and functional organization. Importantly, spatial variation of intrinsic dynamics parallels spatial variation of meta-analytic cognitive functional activation. These findings demonstrate that structural organization of the brain shapes patterns of intrinsic dynamics, ultimately manifesting as distinct axes of psychological processes.

Materials and methods

Dataset: human connectome project (HCP)

Request a detailed protocol

Following the procedure described in Vos de Wael et al., 2018, we obtained structural and functional magnetic resonance imaging (MRI) data of two sets of healthy young adults (age range 22–35 years) with no familial relationships (neither within nor between sets) as Discovery (n=201) and Validation (n=127) sets from Human Connectome Project (HCP; S900 release [Van Essen et al., 2013]). All four resting state fMRI scans (two scans (R/L and L/R phase encoding directions) on day 1 and two scans (R/L and L/R phase encoding directions) on day 2, each about 15 min long; TR=720 ms), as well as structural MRI and diffusion weighted imaging (DWI) data were available for all participants.

HCP data processing

All the structural and functional MRI data were pre-processed using HCP minimal pre-processing pipelines (Van Essen et al., 2013; Glasser et al., 2013). We provide a brief description of data pre-processing below, while detailed information regarding data acquisition and pre-processing is available elsewhere (Van Essen et al., 2013; Glasser et al., 2013). The procedure was separately repeated for Discovery and Validation sets.

Structural MRI

Request a detailed protocol

T1- and T2- weighted MR images were corrected for gradient nonlinearity, and when available, the images were co-registered and averaged across repeated scans for each individual. The corrected T1w and T2w images were co-registered and cortical surfaces were extracted using FreeSurfer 5.3.0-HCP (Vos de Wael et al., 2018; Dale et al., 1999; Fischl et al., 1999). For each individual, cortical thickness was estimated as the difference between pial and white matter surfaces and T1w/T2w ratio was calculated as a putative proxy for intracortical myelin content. The pre-processed data were parcellated into 400 cortical areas using Schaefer parcellation (Schaefer et al., 2018).

Resting state functional MRI

Request a detailed protocol

All 3T functional MRI time-series were corrected for gradient nonlinearity, head motion using a rigid body transformation, and geometric distortions using scan pairs with opposite phase encoding directions (R/L, L/R) (Vos de Wael et al., 2018). Further pre-processing steps include co-registration of the corrected images to the T1w structural MR images, brain extraction, normalization of whole brain intensity, high-pass filtering (>2000s FWHM; to correct for scanner drifts), and removing additional noise using the ICA-FIX process (Vos de Wael et al., 2018; Salimi-Khorshidi et al., 2014). The pre-processed time-series were then parcellated into 400 areas as described above. The parcellated time-series were used to construct functional connectivity matrices as a Pearson correlation coefficient between pairs of regional time-series for each of the four scans of each participant. A group-average functional connectivity matrix was constructed as the mean functional connectivity across all individuals and scans.

Diffusion weighted imaging (DWI)

Request a detailed protocol

DWI data was pre-processed using the MRtrix3 package (Tournier et al., 2019) (https://www.mrtrix.org/). More specifically, fiber orientation distributions were generated using the multi-shell multi-tissue constrained spherical deconvolution algorithm from MRtrix (Dhollander et al., 2016; Jeurissen et al., 2014). White matter edges were then reconstructed using probabilistic streamline tractography based on the generated fiber orientation distributions (Tournier et al., 2010). The tract weights were then optimized by estimating an appropriate cross-section multiplier for each streamline following the procedure proposed by Smith and colleagues (Smith et al., 2015) and a connectivity matrix was built for each participant using the same parcellation as described above. Finally, we used a consensus approach to construct a binary group-level structural connectivity matrix, preserving the edge length distribution in individual participants (Mišić et al., 2015; Betzel et al., 2019; Shafiei et al., 2020; Mišic et al., 2018).

Replication dataset: Midnight Scan Club (MSC)

Request a detailed protocol

We used resting state fMRI data of n=10 healthy young adults, each with 10 scan sessions of about 30 min long, from Midnight Scan Club (MSC Gordon et al., 2017) dataset as an independent replication dataset. Details about the participants, MRI acquisition, and data pre-processing are provided by Gordon and colleagues elsewhere (Gordon et al., 2017). We obtained the surface-based, pre-processed resting state fMRI time courses in CIFTI format through OpenNeuro (https://openneuro.org/datasets/ds000224/versions/1.0.0). The pre-processing steps include motion correction and global signal regression (Gordon et al., 2017). Following the pre-processing methods suggested by Gordon et al., 2017, we smoothed the surface-level time-series data with geodesic 2D Gaussian kernels (σ=2.55 mm) using the Connectome Workbench (Marcus et al., 2011). Finally, we censored the motion-contaminated frames of time-series for each participant separately, using the temporal masks provided with the dataset. The pre-processed data were parcellated into 400 cortical regions using Schaefer parcellation (Schaefer et al., 2018). One participant (MSC08) was excluded from subsequent analysis due to low data reliability and self-reported sleep as described in Gordon et al., 2017. The parcellated time-series were then subjected to the same analyses that were performed on the HCP Discovery and Validation datasets.

Microarray expression data: Allen Human Brain Atlas (AHBA)

Request a detailed protocol

Regional microarray expression data were obtained from six post-mortem brains provided by the Allen Human Brain Atlas (AHBA; http://human.brain-map.org/) (Hawrylycz et al., 2012). We used the abagen (https://github.com/netneurolab/abagen; Markello et al., 2020) toolbox to process and map the data to 400 parcellated brain regions from Schaefer parcellation (Schaefer et al., 2018).

Briefly, genetic probes were reannotated using information provided by Arnatkeviciute et al., 2019 instead of the default probe information from the AHBA dataset. Using reannotated information discards probes that cannot be reliably matched to genes. The reannotated probes were filtered based on their intensity relative to background noise levels (Quackenbush, 2002); probes with intensity less than background in 50% of samples were discarded. A single probe with the highest differential stability, ΔS(p), was selected to represent each gene (Hawrylycz et al., 2015), where differential stability was calculated as:

(1) ΔS(p)=1(N2)i=1N1j=i+1Nρ[Bi(p),Bj(p)]

Here, ρ is Spearman’s rank correlation of the expression of a single probe p across regions in two donor brains, Bi and Bj, and N is the total number of donor brains. This procedure retained 15,656 probes, each representing a unique gene.

Next, tissue samples were mirrored across left and right hemispheres (Romero-Garcia et al., 2018) and then assigned to brain regions using their corrected MNI coordinates (https://github.com/chrisfilo/alleninf) by finding the nearest region, up to 2 mm away. To reduce the potential for misassignment, sample-to-region matching was constrained by hemisphere and cortical/subcortical divisions (Arnatkeviciute et al., 2019). If a brain region was not assigned any sample based on the above procedure, the sample closest to the centroid of that region was selected in order to ensure that all brain regions were assigned a value. Samples assigned to the same brain region were averaged separately for each donor. Gene expression values were then normalized separately for each donor across regions using a robust sigmoid function and rescaled to the unit interval (Fulcher and Fornito, 2016). Scaled expression profiles were finally averaged across donors, resulting in a single matrix with rows corresponding to brain regions and columns corresponding to the retained 15,656 genes. The expression values of 1906 brain-specific genes were used for further analysis (Burt et al., 2018).

Massive temporal feature extraction using hctsa

Request a detailed protocol

We used the highly comparative time-series analysis toolbox, hctsa (Fulcher et al., 2013; Fulcher and Jones, 2017), to perform a massive feature extraction of the time-series of each brain area for each participant. The hctsa package extracted over 7000 local time-series features using a wide range of operations based on time-series analysis (Fulcher et al., 2013; Fulcher and Jones, 2017). The extracted features include, but are not limited to, distributional properties, entropy and variability, autocorrelation, time-delay embeddings, and nonlinear properties of a given time-series (Fulcher et al., 2013; Fulcher, 2018).

The hctsa feature extraction analysis was performed on the parcellated fMRI time-series of each run and each participant separately (Figure 1). Following the feature extraction procedure, the outputs of the operations that produced errors were removed and the remaining features (6441 features) were normalized across nodes using an outlier-robust sigmoidal transform. We used Pearson correlation coefficients to measure the pairwise similarity between the time-series features of all possible combinations of brain areas. As a result, a temporal profile similarity network was constructed for each individual and each run, representing the strength of the similarity of the local temporal fingerprints of brain areas (Figure 1). The resulting similarity matrices were then compared to the underlying functional and structural brain networks.

Neurosynth

Request a detailed protocol

Functional activation probability maps were obtained for multiple psychological terms using Neurosynth (Yarkoni et al., 2011) (https://github.com/neurosynth/neurosynth). Probability maps were restricted to those for terms present in both Neurosynth and the Cognitive Atlas (Poldrack et al., 2011), yielding a total of 123 maps (Supplementary file 2). We used the volumetric ‘association test’ (i.e. reverse inference) maps, which were projected to the FreeSurfer fsaverage5 mid-grey surface with nearest neighbor interpolation using Freesurfer’s mri_vol2surf function (v6.0.0; http://surfer.nmr.mgh.harvard.edu/). The resulting surface maps were then parcellated to 400 cortical regions using the Schaefer parcellation (Schaefer et al., 2018).

Null model

Request a detailed protocol

A consistent question in the present work is the topographic correlation between time-series features and other features of interest. To make inferences about these links, we implement a null model that systematically disrupts the relationship between two topographic maps but preserves their spatial autocorrelation (Alexander-Bloch et al., 2018; Markello and Misic, 2020) (see also Burt et al., 2018; Burt et al., 2020 for an alternative approach). We first created a surface-based representation of the Cammoun atlas on the FreeSurfer fsaverage surface using the Connectome Mapper toolkit (https://github.com/LTS5/cmp; Daducci et al., 2012). We used the spherical projection of the fsaverage surface to define spatial coordinates for each parcel by selecting the vertex closest to the center-of-mass of each parcel (Vázquez-Rodríguez et al., 2019; Shafiei et al., 2020; Vézquez-Rodríguez et al., 2020). The resulting spatial coordinates were used to generate null models by applying randomly-sampled rotations and reassigning node values based on the closest resulting parcel (10,000 repetitions). The rotation was applied to one hemisphere and then mirrored to the other hemisphere.

Data availability

All data used in this study is publicly available. Detailed information about the datasets is available in the manuscript.

The following previously published data sets were used

References

  1. Conference
    1. Dhollander T
    2. Raffelt D
    3. Connelly A
    (2016)
    Unsupervised 3-tissue response function estimation from single-shell or multi-shell diffusion mr data without a co-registered t1 image
    ISMRM Workshop on Breaking the Barriers of Diffusion MRI.
  2. Conference
    1. Fulcher BD
    (2018)
    Feature-based time-series analysis
    Feature Engineering for Machine Learning and Data Analytics.
    1. Garrett DD
    2. Epp SM
    3. Lindenberger U
    (2017)
    Localtemporal variability reflects functional network integrationin the human brain: on the crucial role of thethalamus
    NeuroImage 183:e19.
    1. Hilgetag CC
    2. Goulas A
    (2020) 'Hierarchy' in the organization of brain networks
    Philosophical Transactions of the Royal Society B: Biological Sciences 375:20190319.
    https://doi.org/10.1098/rstb.2019.0319
  3. Conference
    1. Langs G
    2. Golland P
    3. Ghosh SS
    (2015)
    Predicting activation across individuals with resting-state functional connectivity based multi-atlas label fusion
    International Conference on Medical Image Computing and Computer-Assisted Intervention. pp. 313–320.
  4. Conference
    1. Lurie DJ
    2. D’Esposito M
    (2020)
    Gradients of intrinsic dynamics follow connectomic, anatomical and microstructural hierarchies
    Annual Meeting of the Organization for Human Brain Mapping.
    1. Mesulam M-M
    (2000)
    Behavioral neuroanatomy
    Principles of behavioral and cognitive neurology 2:1–120.
    1. Pedregosa F
    2. Varoquaux G
    3. Gramfort A
    4. Michel V
    5. Thirion B
    6. Grisel O
    7. Blondel M
    8. Prettenhofer P
    9. Weiss R
    10. Dubourg V
    (2011)
    Scikit-learn: machine learning in Python
    Journal of Machine Learning Research 12:2825–2830.
  5. Conference
    1. Tournier JD
    2. Calamante F
    3. Connelly A
    (2010)
    Improved probabilistic streamlines tractography by 2nd order integration over fibre orientation distributions
    Proceedings of the International Society for Magnetic Resonance in Medicine.
  6. Book
    1. von Economo CF
    2. Koskinas GN
    3. Triarhou LC
    (2008)
    Atlas of Cytoarchitectonics of the Adult Human Cerebral Cortex, 10
    Karger Basel.
  7. Book
    1. von Economo CF
    2. Koskinas GN
    (1925)
    Die cytoarchitektonik der hirnrinde des erwachsenen menschen
    J. Springer.

Decision letter

  1. Lucina Q Uddin
    Reviewing Editor; University of Miami, United States
  2. Chris I Baker
    Senior Editor; National Institute of Mental Health, National Institutes of Health, United States
  3. Lucina Q Uddin
    Reviewer; University of Miami, United States
  4. Maxwell Bertolero
    Reviewer; University of Pennsylvania, United States

In the interests of transparency, eLife publishes the most substantive revision requests and the accompanying author responses.

Acceptance summary:

We believe that your work contributes novel insights into functional brain organization.

Decision letter after peer review:

Thank you for submitting your article "Topographic gradients of intrinsic dynamics across neocortex" for consideration by eLife. Your article has been reviewed by three peer reviewers, including Lucina Q Uddin as the Reviewing Editor and Reviewer #1, and the evaluation has been overseen by Chris Baker as the Senior Editor. The following individual involved in review of your submission has agreed to reveal their identity: Maxwell Bertolero (Reviewer #2).

The reviewers have discussed the reviews with one another and the Reviewing Editor has drafted this decision to help you prepare a revised submission.

We would like to draw your attention to changes in our revision policy that we have made in response to COVID-19 (https://elifesciences.org/articles/57162). Specifically, we are asking editors to accept without delay manuscripts, like yours, that they judge can stand as eLife papers without additional data, even if they feel that they would make the manuscript stronger. Thus the revisions requested below only address clarity and presentation.

Summary:

This manuscript describes a study that examines topographic gradients of intrinsic dynamics in the human brain. The authors identify two gradients with distinct temporal compositions and gene expression, myelin, cortical thickness, network embedding, and functional activation properties. The reviewers agreed that the reproducibility analyses (replication on independent datasets, use of two different parcellation schemes), and comparisons of data with and without (grey-matter signal regression) are particularly compelling. The similarities in temporal features between distinct brain regions revealed a ventromedial-dorsolateral (PC1) and a unimodal-transmodal (PC2) topographic gradients. Interestingly, analyses suggest a weak correlation between "temporal profile similarities" and functional connectivity. Broadly, these initial set of results suggest that the organization of spontaneous BOLD signal fluctuations link to both the geometrical and topological embedding of regions in macroscopic networks. The second set of analyses suggest that various measures of signal autocorrelations link to PC1 while measures of distribution shape (dynamic range) link to PC2. Moreover, the detected gradients are correlated with the spatial distribution of gene expression, cortical thickness, T1w/T2w, and functional hierarchy in an opposite fashion. The third set of analyses focus on linking topographic patterns in spontaneous fMRI signals with probability maps of task-based fMRI signals. The output of this analysis reveals that PC1 recapitulates a cognitive-affective axis, whereas the PC2 link to a sensory-cognitive axis. Overall, the reviewers were all in agreement regarding the execution, significance and approach of the study.

Essential revisions:

1) The overall narrative and rationale for the approach is missing. Some additional comments and suggestions regarding conceptual framing are below:

a) The Introduction as it is written includes a great deal of jargon that should be unpacked. In particular, the terms "contextual information", "network embedding", "spectral power", and "nonlinear dynamic models" should be defined clearly as they are introduced.

b) Cognitive ontologies come out of the blue in the Results, and are not mentioned at all in the Introduction. Similarly, the significance of gene expression and its relationship to network structure is not at all discussed in the Introduction. Again in the Discussion, the rationale, significance, and relevance of the gene expression findings are missing.

c) The paper provides numerous "confirmatory" findings, but after reading the paper a few times, I am left wondering what new fundamental knowledge the study provides. We know that fluctuations in spontaneous fMRI signal are patterned across the cortex, and that such patterns link to genetics, microstructural gradients and macroscopic functional and structural patterns. Perhaps the most novel/unexpected result is that biologically meaningful time-series features are not associated with functional connectivity. However, this result is only superficially discussed. The authors should make an extra effort to highlight better how the paper confirms (still a vital endeavour) and extends existing knowledge.

d) I am confused about the use of the term "dynamic". The adopted measures have little to do with dynamic measures of neural synchronization per se; they are more summary measure of these dynamics. The authors should be careful in their terminology and discuss the results for what they are (i.e., summary measures of spontaneous changes in fMRI signal).

2) Points of clarification and queries regarding analyses are below:

a) It seems that temporal features (autocorrelation, variance, spectral power, etc) are all considered together in the hctsa analysis. How do the authors account for the fact that these features are not necessarily independent? In other words, does the analysis consider the inter-relationships between these individual features?

b) Perhaps consider including a comparison evolutionary expansion or the participation coefficient (PC2 is likely anti-correlated with these).

c) The statistics reporting needs work. Please report all degrees of freedom. Also, please either report the results in the figure that shows the data, or in the text. In many places, the statistics are divided up, and that makes it hard to read. Also, there seem to be many points where a t-test was done, but only a p-value, not a t-value, is reported. Just saying p is close to zero is not informative if I don't know the DoF and the t value. For each finding, report the test type, the test value, the p value (with DoF), and then reference the panel. Or just reference the panel and put those stats in it. Or both! Also, don't put a * if you are not going to explain what it means.

d) Several temporal features do not seem orthogonal (e.g., power spectrum and autocorrelations). I may have missed this information but, were redundant features removed to obtain a balanced set of features describing the resting-state fMRI signal? If similar features were not removed, results might have been biased towards aspects of rest fMRI signals that are over-represented. The Results section suggests that this may be the case.

e) Patterns of structural connectivity are known to be linked to spontaneous fluctuations in functional connectivity. I suggest extending these results by providing a quantification of the shared/distinct contributions of functional network topology and anatomical factors to temporal profile similarity.

f) Could the authors consider using an estimate of functional connectivity that accounts for shared signal (e.g., partial correlation, multiple regression)? I think this would shed light on the relationships identified in the paper and how much can be related to likely false positives found in correlation FC.

g) There is little discussion regarding why the variables (e.g., Figure 4) should (or should not) be related to each other (see also reviewer comment a). Given that many of these measures are inter-related, it is not entirely clear to me if they are meaningful. An example of this is the positive/negative relationships observed in Figure 4: Is this bound to happen because of the relationship between PC1 and PC2? Is this an interesting relationship, or just a statistical consequence of PCA?

https://doi.org/10.7554/eLife.62116.sa1

Author response

Essential revisions:

1) The overall narrative and rationale for the approach is missing. Some additional comments and suggestions regarding conceptual framing are below:

a) The Introduction as it is written includes a great deal of jargon that should be unpacked. In particular, the terms "contextual information", "network embedding", "spectral power", and "nonlinear dynamic models" should be defined clearly as they are introduced.

We have now clarified these definitions and have modified the manuscript accordingly. Specific changes are explained below (Introduction):

Re: “contextual information”:

“The primary functional consequence of this hierarchy of timescales is thought to be a hierarchy of temporal receptive windows: time windows in which a newly arriving stimulus will modify processing of previously presented (i.e. contextual) information (Hasson et al., 2008; Honey et al., 2012; Baldassano et al., 2017; Huntenburg et al., 2018; Chaudhuri et al., 2015; Chien and Honey, 2020). […] Altogether, previous work highlights a hierarchy of a small number of manually selected specific time-series features, but it is possible that different types of local computations manifest as different organizational gradients.”

Re: “network embedding”:

“The signal variability of brain areas, measured in terms of standard deviations or temporal entropy, is closely related to their structural and functional connectivity profiles (i.e. network embedding) (Misic et al., 2011; Burzynska et al., 2013; Garrett et al., 2017; Shafiei et al., 2019).”

Re: “Spectral power”:

“A significant limitation is that conventional computational analysis is based on specific, manually selected time-series features, such as the decay of the autocorrelation function, bands of the Fourier power spectrum, or signal variance.”

Re: “nonlinear dynamic models”:

“Finally, in computational models of structurally coupled neuronal populations (neural mass and neural field models; Breakspear, 2017), highly interconnected hubs exhibit slower dynamic fluctuations, while sensory areas exhibit fast fluctuating neural activity (Gollo et al., 2015). Indeed, these models offer better fits to empirical functional connectivity if they assume heterogeneous local dynamics that follow a unimodal—transmodal gradient (Cocchi et al., 2016; Demirtas et al., 2019; Wang et al., 2019; Deco et al., 2020).”

b) Cognitive ontologies come out of the blue in the Results, and are not mentioned at all in the Introduction. Similarly, the significance of gene expression and its relationship to network structure is not at all discussed in the Introduction. Again in the Discussion, the rationale, significance, and relevance of the gene expression findings are missing.

We agree, and have modified the manuscript accordingly. Specific changes are explained below:

Re: cognitive ontologies:

We first introduce the rationale for the analysis; namely, that spatial variation in circuit structure may shape circuit dynamics and computations, manifesting as patterns of regional functional specialization (Introduction):

“Altogether, multiple lines of evidence suggest that local computations may reflect systematic variation in microscale properties and macroscale network embedding, manifesting as diverse time-series features of regional neural activity. How molecular, cellular and connectomic architecture precisely shapes temporal dynamics, and ultimately, cortical patterns of functional specialization, is poorly understood.”

We then explicitly lay out how we will address this question: by mapping time-series features to reverse-inference meta-analytic maps (Introduction):

“Finally, we map time-series features to a meta-analytic atlas of cognitive ontologies to investigate how temporal dynamics shape regional functional specialization. […] These spatial variations in intrinsic dynamics ultimately manifest as patterns of distinct psychological functions.”

Re: gene expression

We have modified the Introduction to first intuitively set the idea that gene expression may convey information about cellular structure and function, and may therefore be related to intrinsic dynamics (Introduction):

“For a neuronal population, the confluence of local micro-architectural properties and global connectivity shapes both the generation of local rhythms, as well as its propensity to communicate with other populations. […] These micro-architectural properties – increasingly measured directly from histology or inferred from other measurements, such as microarray gene expression – provide a unique opportunity to relate circuit architecture to temporal dynamics and computation.”

We have also substantially revised the Discussion to explain the significance of these findings. We discuss this change in more detail in our response to comment #2, point g, but have included the new text here as well (Discussion section):

“Applying a data-driven feature extraction method to high-resolution BOLD fMRI, we decompose regional signals into two intrinsic modes, with distinct topographic organization and time-series features. […] Taken together, we find evidence that molecular and cellular properties (gene expression PC1) relate to regional autocorrelation, while micro-circuit properties (T1w/T2w, cortical thickness) and macroscale network embedding (principal functional gradient) relate to regional dynamic range.”

c) The paper provides numerous "confirmatory" findings, but after reading the paper a few times, I am left wondering what new fundamental knowledge the study provides. We know that fluctuations in spontaneous fMRI signal are patterned across the cortex, and that such patterns link to genetics, microstructural gradients and macroscopic functional and structural patterns. Perhaps the most novel/unexpected result is that biologically meaningful time-series features are not associated with functional connectivity. However, this result is only superficially discussed. The authors should make an extra effort to highlight better how the paper confirms (still a vital endeavour) and extends existing knowledge.

We concur that the original manuscript did not fully highlight the novelty of the findings and conceptual advances made by the study. In the revised manuscript, we emphasize 5 key contributions:

1) Comprehensive time-series phenotyping. The majority of the field focuses on specific univariate time-series features-of-interest, such as measures of intrinsic time scale, variance, or spectral power. Our approach is the first to comprehensively benchmark the entire dynamic profile of the brain, by near-exhaustively estimating 6000+ features from the wider time-series literature. This represents a significant conceptual advance, as it offers a fundamentally new and entirely data-driven method to quantify fMRI dynamics. Focusing on univariate time-series features is a powerful but under-utilized direction to understand neural dynamics.

A salient example is that when multiple time-series features are taken into account, we observe two prominent yet distinct modes of intrinsic dynamics, with different topographies and dynamical signatures. Consider the spatial variation of BOLD autocorrelation; multiple accounts posit that the unimodal-transmodal hierarchy is expressed as regional differences in autocorrelation or intrinsic timescale (Hasson et al., 2008; Murray et al., 2014; Honey et al., 2012; Ito et al., 2020). Yet the present results show that this unimodal-transmodal hierarchy is more closely related to regional differences in dynamic range, whereas regional differences in autocorrelation display a different spatial pattern. Indeed, the dynamic range pattern (PC2) more closely corresponds to multiple measures of microarchitecture than the autocorrelation pattern (PC1). By comprehensively studying multiple time-series features and multiple anatomical features, the present report shows that — what is typically assumed to be a single hierarchy — is actually two hierarchies with distinct structural, dynamic and psychological attributes. To clarify this contribution, we have added the following text to the Discussion (Discussion section):

“Our results demonstrate that regional haemodynamic activity, often overlooked in favour of electrophysiological measurements with greater temporal resolution, possesses a rich dynamic signature (Garrett et al., 2013; Uddin, 2020; Preti et al., 2017; Lurie et al., 2020; Li et al., 2019; Bolt et al., 2018). […] As we discuss below, feature-based time-series phenotyping offers a powerful, fundamentally new and entirely data-driven method to quantify and articulate neural dynamics.”

2) Systematic mapping of time-series features to architectural features. A related point is that most previous literature has sought to relate small manually selected sets of time-series features-of-interest to small manually selected sets of structural features-of-interest. The present work integrates comprehensive dynamic profiles to architectural features across multiple scales, including molecular, cellular, laminar and connectomic features. To our knowledge, there exist no reports or records of these relationships, so the present findings will be a touchstone for future research in this domain. To clarify this contribution, we have added the following paragraph to the Discussion (Discussion section): (see also our response to comment #2, point g)

“Applying a data-driven feature extraction method to high-resolution BOLD fMRI, we decompose regional signals into two intrinsic modes, with distinct topographic organization and time-series features[…] How these patterns are related to underlying cell types and subcortical afferent input – in particular, thalamocortical feedback – is an important ongoing question (Abeysuriya et al., 2015; Garrett et al., 2018; Shine et al., 2019; Wang et al., 2019; Muller et al., 2020; Paquola et al., 2020).”

3) Linking structure, dynamics, and cognition. We show that topographic variations in microcircuitry and connectomic embedding yield variations in intrinsic dynamics and, importantly, may explain regional differences in functional specialization. Specifically, we show that patterns of intrinsic dynamics are concomitant with two dominant axes of functional activation, derived from Neurosynth. Thus, the present findings conceptually link structural, temporal and neurocognitive architectures. To clarify this contribution, we have added the following paragraph to the Discussion (Discussion section):

“Importantly, the two patterns are related to two dominant axes of meta-analytic functional activation. […] Collectively, these results provide evidence that local computations reflect systematic variation in multiple anatomical circuit properties, and can be measured as unique temporal signatures in regional activity and patterns of functional specialization.”

4) Re-conceptualizing functional connectivity. At present, functional relationships are typically conceptualized as bivariate relationships among regional time-series, with the implicit assumption that if two regions display coherent dynamics, they may be engaged in a common function. Here we show that, by deconstructing time-series into 6000+ features, it is possible to operationalize inter-regional relationships as similarity of their dynamic properties. Importantly, we show that traditional, correlation-based FC is only weakly associated with temporal profile similarity. In other words, the traditional conceptualization of functional connectivity misses out on an important set of inter-regional relationships. Namely, two regions may display identical time-series profiles, suggesting common circuit dynamics and function, but unless they also display time-locked activity, current methods would miss out on this potentially biologically meaningful inter-regional relationship. We have clarified this contribution in the revised manuscript (Discussion section):

“Indeed, we find that two regions are more likely to display similar intrinsic dynamics if they are anatomically connected and if they are part of the same functional community, suggesting that network organization and local intrinsic dynamics are intertwined (Gollo et al., 2015; Cocchi et al., 2016). […] Namely, two regions may display identical time-series profiles, suggesting common circuit dynamics and function, but unless they also display time-locked activity, current methods would miss out on this potentially biologically meaningful inter-regional relationship.”

5) Extensive robustness testing. The results are subjected to numerous sensitivity and replication analyses, which demonstrate that the findings are consistent across:

a) Group- and individual-level

b) With and without grey-matter signal regression

c) Parcellation (anatomical and functional) and parcellation resolution

d) Out-of-sample validation in held-out HCP data and MSC data

We have added the following sentence to the Discussion (Discussion section):

“These findings are robust against a wide range of methodological choices and were validated in two held-out samples.”

d) I am confused about the use of the term "dynamic". The adopted measures have little to do with dynamic measures of neural synchronization per se; they are more summary measure of these dynamics. The authors should be careful in their terminology and discuss the results for what they are (i.e., summary measures of spontaneous changes in fMRI signal).

We have added the following sentence in the Introduction (at the first mention of the term “dynamic”) to clarify that the time-series features are summary measures of local neural activity (Introduction):

“Here we comprehensively chart summary features of spontaneous BOLD signals across the cerebral cortex (hereafter referred to as “intrinsic dynamics”), mapping temporal organization to structural organization.”

We have also changed the term “dynamic profile” to “time-series profile” throughout the manuscript.

2) Points of clarification and queries regarding analyses are below:

a) It seems that temporal features (autocorrelation, variance, spectral power, etc) are all considered together in the hctsa analysis. How do the authors account for the fact that these features are not necessarily independent? In other words, does the analysis consider the inter-relationships between these individual features?

This is precisely why we applied PCA to the hctsa output – to identify groups of correlated features. The analysis produces 6,441 time-series features that are nearly-exhaustive but some may potentially be correlated. As a result, we sought to find linear combinations of features that explain maximal variance (principal components) and that span different conceptual types of time-series properties.

To clarify this rationale, we have added the following passage to the revised manuscript (Results section, “Two distinct spatial gradients of intrinsic dynamics” subsection):

“The hctsa library generates 6,441 time-series features, with the aim of being comprehensive in coverage across scientific time-series analysis algorithms and, as a result, contains groups of features with correlated outputs (Fulcher et al., 2013). We therefore sought to identify groups of correlated features that explain maximal variance and that span different conceptual types of time-series properties.”

To provide an intuitive sense of whether top-loading features in each component (full list in Supplementary files 3, 4) display fMRI-specific patterning, we additionally plotted two high-loading, representative features for each component in Figure 3C (top). To build intuition about what each component reflects about regional signals, we selected three regional time-series from one participant based on their lag-1 autocorrelation and kurtosis (Figure 3C; circles on the brain surface: pink = 5th percentile, green = 50th percentile, purple = 95th percentile). For the former, going from low-ranked to high-ranked regions results in a slowing down of BOLD fluctuations. For the latter, going from low-ranked to high-ranked regions results in increasingly heavier symmetric tails of the signal amplitude distributions.

b) Perhaps consider including a comparison evolutionary expansion or the participation coefficient (PC2 is likely anti-correlated with these).

We thank the reviewers for their suggestion. We have now added an analysis (Figure 4—figure supplement 2) comparing the two maps with the PCA results. The evolutionary expansion map (Hill et al., 2010; Baum et al., 2020) was obtained through https://github.com/PennLINC/Brain_Organization and parcellated into 400 cortical areas using the Schaefer parcellation (Schaefer et al., 2018). Regional weighted participation coefficients (Brain Connectivity Toolbox; Rubinov and Sporns, 2010) were estimated from functional connectivity graphs with respect to the 7-network partition of intrinsic networks (Yeo et al., 2011; Schaefer et al., 2018). We then compared the maps to the PC1 and PC2 brain score patterns. The evolutionary expansion map was significantly correlated with PC2 topography (rs = 0.52, pspin = 0.0002), but not PC1 topography (rs = 0.01, pspin = 0.95). The regional participation coefficient map was not significantly correlated with either PC1 or PC2 topography (rs = -0.14, pspin = 0.16; and rs = -0.10, pspin = 0.32; respectively).

We have also added the following passage to the revised manuscript (“Results section, “Intrinsic dynamics reflect microscale and macroscale hierarchies” subsection):

“For completeness, we also tested associations with two maps that were previously related to cortical hierarchies: evolutionary expansion (indexing enlargement of cortical areas in the human relative to the macaque) (Hill et al., 2010; Baum et al., 2020) and node-wise functional participation coefficient (indexing the diversity of a node's links) (Bertolero et al., 2017; Baum et al., 2020). PC2 is significantly correlated with evolutionary expansion (rs = 0.52, pspin = 0.0002), but neither component is correlated with participation coefficient (Figure 4—figure supplement 2).”

c) The statistics reporting needs work. Please report all degrees of freedom. Also, please either report the results in the figure that shows the data, or in the text. In many places, the statistics are divided up, and that makes it hard to read. Also, there seem to be many points where a t-test was done, but only a p-value, not a t-value, is reported. Just saying p is close to zero is not informative if I don't know the DoF and the t value. For each finding, report the test type, the test value, the p value (with DoF), and then reference the panel. Or just reference the panel and put those stats in it. Or both! Also, don't put a * if you are not going to explain what it means.

Thank you for pointing this out. We have added the following explanation to the Figure 2 legend to denote what asterisks represent:

“(c, d) Regional time-series features are compared between pairs of cortical areas using their structural and functional connectivity profiles. […] For functional networks, statistical significance of the difference of the mean temporal profile similarity of within and between intrinsic networks is also assessed against a null distribution of differences generated by spatial autocorrelation-preserving label permutation (“spin tests”; Alexander-Bloch et al., 2018) (d, right-most panel).”

We have also included complete information for all t-tests, including degrees of freedom and exact t-values (Results section, “Inter-regional temporal profile similarity reflects network geometry and topology” subsection):

Structurally connected vs. not connected temporal profile similarity (Figure 2C):

t(79,798) = 40.234; p ≅ 0

for distance regressed temporal profile similarity: t(79,798) = 9.916; p ≅ 0

Within vs. between intrinsic network temporal profile similarity (Figure 2D):

t(79,798) = 61.093; p ≅ 0

for distance regressed temporal profile similarity: t(79,798) = 47.112; p ≅ 0

We note that, since inference is performed across a large number of edges, there are correspondingly large degrees of freedom. To ensure that the results are not simply an artefact of sample size, we additionally used nonparametric tests, for which we show null distributions and empirical values in the right-most panels of the figure. Specifically, we used rewired networks for estimating similarity of structurally-connected and -unconnected node pairs, and we used spatial autocorrelation-preserving spin tests to assess the similarity of within- and between-functional network node pairs (both panels shown in Figure 2C, D).

d) Several temporal features do not seem orthogonal (e.g., power spectrum and autocorrelations). I may have missed this information but, were redundant features removed to obtain a balanced set of features describing the resting-state fMRI signal? If similar features were not removed, results might have been biased towards aspects of rest fMRI signals that are over-represented. The Results section suggests that this may be the case.

No features were removed, as the goal of the analysis was to comprehensively characterize the temporal organization of neural activity, without any pre-selection of features. As we outline below, this approach opens fundamentally new ways to analyze neural activity, but potential correlations among features must be explicitly considered and analyzed.

Although some specific hctsa measures attempt to capture similar properties, each one uses a different method and/or parameter settings, yielding 6,441 unique features. At the same time, features could be correlated with each other (e.g. outputs of autocorrelation measures) (Fulcher et al., 2013). This is precisely why we applied PCA to the hctsa output – to identify groups of correlated features. Specifically, we sought to find linear combinations of features that explain maximal variance (principal components) and that span different conceptual types of time-series properties. This approach allows us to discover connections between time-series properties that have different names and come from different literatures, but behave similarly. To clarify this rationale, we have added the following passage to the revised manuscript (Results section, “Two distinct spatial gradients of intrinsic dynamics” subsection):

“The hctsa library generates 6,441 time-series features, with the aim of being comprehensive in coverage across scientific time-series analysis algorithms and, as a result, contains groups of features with correlated outputs (Fulcher et al., 2013). We therefore sought to identify groups of correlated features that explain maximal variance and that span different conceptual types of time-series properties.”

Analyzing all features allowed us to map the temporal organization of neural activity across the brain. However, the downside of using all possible features is that some will be correlated with each other. Although we applied PCA to guard against this, the solutions may be biased towards the more numerous features. To acknowledge this limitation, we have added the following text to the Discussion (Discussion section):

“Fourth, the analysis included all features from hctsa, potentially biasing results towards specific properties of BOLD signals. […] This may obscure the contribution of under-represented feature classes, and should be investigated further in future work.”

e) Patterns of structural connectivity are known to be linked to spontaneous fluctuations in functional connectivity. I suggest extending these results by providing a quantification of the shared/distinct contributions of functional network topology and anatomical factors to temporal profile similarity.

Inter-regional temporal profile similarity could potentially be driven by inter-regional (a) distance, (b) structural connectivity, (c) functional connectivity. To quantify the distinct contributions of each factor, we first constructed a multiple linear regression with distance, structural connectivity and functional connectivity as independent variables and temporal profile similarity as the dependent variable. Observations were edges or node pairs. The analysis was restricted to structurally connected node pairs (k = 4,950 unique edges). We then applied Dominance Analysis, a method for assessing the relative importance of predictors (Budescu, 1993; Azen and Budescu, 2003; https://github.com/dominance-analysis/dominance-analysis).

Briefly, the analysis measures relative importance in a series of pairwise comparisons, where two predictors are compared in the context of all possible models that contain some subset of the other predictors (i.e. for p predictors we have 2p-1 models). The incremental R2 contribution of each predictor to the subset model of all other predictors is then computed. The Total Dominance of each predictor is summarized as the additional contribution of a predictor to all subset models.

The results of the analysis are shown in Supplementary file 1. The total R2 is 0.28 for the complete model, with the following relative importance (last column): distance – 56%, FC – 23.6%, SC – 20.4%. In other words, distance contributes the most to temporal profile similarity, while structural and functional connectivity make distinct but approximately even contributions.

We have also added the following passage (Results section, “Inter-regional temporal profile similarity reflects network geometry and topology” subsection):

“As a final step, we sought to assess the distinct contributions of distance, structural connectivity and functional connectivity to temporal profile similarity. Dominance analysis revealed the relative importance of each predictor (collective R2 = 0.28; distance = 56%, structural connectivity = 20.4%, functional connectivity = 23.6%; Supplementary file 1), suggesting that distance contributes the most to temporal profile similarity, while structural and functional connectivity make distinct but approximately even contributions (Budescu, 1993; Azen and Budescu, 2003; https://github.com/dominance-analysis/dominance-analysis).”

f) Could the authors consider using an estimate of functional connectivity that accounts for shared signal (e.g., partial correlation, multiple regression)? I think this would shed light on the relationships identified in the paper and how much can be related to likely false positives found in correlation FC.

We have repeated all FC-based analyses using partial correlations, as implemented in Nilearn (Abraham et al., 2014). The results are virtually identical to before. Namely, there is a weak relationship between functional connectivity and temporal profile similarity.

We have included the new results as a Figure 2—figure supplement 2 and have added the following sentences to the revised manuscript to refer to this figure (Results section, “Inter-regional temporal profile similarity reflects network geometry and topology” subsection):

“The results are consistent when functional connectivity is estimated using partial correlations (Figure 2—figure supplement 2).”

g) There is little discussion regarding why the variables (e.g., Figure 4) should (or should not) be related to each other (see also reviewer comment a). Given that many of these measures are inter-related, it is not entirely clear to me if they are meaningful. An example of this is the positive/negative relationships observed in Figure 4: Is this bound to happen because of the relationship between PC1 and PC2? Is this an interesting relationship, or just a statistical consequence of PCA?

Why are these relationships significant?

Existing work has highlighted dominant cortical hierarchies but perhaps different types of local computations manifest in different organizational gradients. Our goal was to systematically compare temporal and micro-architectural features. By comprehensively mapping time-series features to other well-studied maps (gene expression, functional gradient, intracortical myelin, cortical thickness, evolutionary expansion), our results can be readily contextualized by other researchers studying brain organization from molecular, cellular or macroscale perspectives. Throughout the analysis, we apply conservative autocorrelation-preserving spin tests to assess associations.

Specifically, we find two dominant but distinct temporal modes that are expressed over the cortex, one indexing signal autocorrelation and the other indexing dynamic range. The first is closely associated with gene expression PC1 (itself closely related to cell type composition, synaptic physiology and cortical cytoarchitecture; Burt et al., 2018), suggesting a molecular and cellular basis for regional differences in temporal autocorrelation. The second is closely associated with the principal functional gradient, as well as with intracortical myelin and cortical thickness, suggesting that the dynamic range of BOLD signals is related to regional variation in macroscale circuit organization. Taken together, we find evidence that molecular and cellular properties (gene expression PC1) relate to regional autocorrelation, while micro-circuit properties (T1w/T2w, cortical thickness) and macroscale network embedding (principal functional gradient) relate to regional dynamic range.

To clarify these relationships, we have included the following passage to the Discussion (Discussion section):

“Applying a data-driven feature extraction method to high-resolution BOLD fMRI, we decompose regional signals into two intrinsic modes, with distinct topographic organization and time-series features. […] How these patterns are related to underlying cell types and subcortical afferent input — in particular, thalamocortical feedback — is an important ongoing question (Abeysuriya et al., 2015; Garrett et al., 2018; Shine et al., 2019; Wang et al., 2019; Muller et al., 2020; Paquola et al., 2020).”

Interestingly, the present results add subtlety to the emerging literature on multimodal cortical gradients. As the reviewers point out, many reports suggest that these gradients are inter-related, with most assuming a unimodal-transmodal hierarchy. In the time-series dynamics literature, this is partly due to the fact that researchers often study specific time-series features of interest. However, the present results show that when multiple time-series features are taken into account, we observe two prominent yet distinct modes of intrinsic dynamics, with different topographies and dynamical signatures. A salient example is the distribution of regional autocorrelation; multiple accounts posit that the unimodal-transmodal hierarchy is expressed as regional differences in autocorrelation or intrinsic timescale. Yet the present results show that the unimodal-transmodal hierarchy is more closely related to regional differences in dynamic range, whereas regional differences in autocorrelation follow a different topographic distribution. By comprehensively studying multiple time-series features and multiple anatomical features, the present report shows that, what is typically assumed to be a single hierarchy, is actually two hierarchies with distinct structural, dynamic and psychological attributes.

Is this a consequence of PCA?

In some instances, PC1 and PC2 scores have correlations with different signs. This is not an artifact of the principal component analysis. Namely, the signs of PC weights are arbitrary; multiplication by -1 (corresponding to axis reflection), yields components with equal effect size (variance accounted for) (McIntosh and Misic, 2013). The only constraint in the analysis is that the eigenvectors (PC weights) are mutually orthogonal. Thus, there is no a priori reason for PC scores to be positively correlated, negatively correlated, or anticorrelated with a third measure of interest, such as gene expression, the principal gradient, etc.

Reference:

McIntosh, A. R., and Mišić, B. (2013). Multivariate statistical analyses for neuroimaging data. Annual review of psychology, 64, 499-525.

https://doi.org/10.7554/eLife.62116.sa2

Article and author information

Author details

  1. Golia Shafiei

    McConnell Brain Imaging Centre, Montréal Neurological Institute, McGill University, Montréal, Canada
    Contribution
    Conceptualization, Resources, Data curation, Formal analysis, Validation, Investigation, Visualization, Methodology, Writing - original draft, Project administration, Writing - review and editing
    For correspondence
    golia.shafiei@mail.mcgill.ca
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-2036-5571
  2. Ross D Markello

    McConnell Brain Imaging Centre, Montréal Neurological Institute, McGill University, Montréal, Canada
    Contribution
    Data curation, Software, Methodology, Writing - review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-1057-1336
  3. Reinder Vos de Wael

    McConnell Brain Imaging Centre, Montréal Neurological Institute, McGill University, Montréal, Canada
    Contribution
    Data curation, Writing - review and editing
    Competing interests
    No competing interests declared
  4. Boris C Bernhardt

    McConnell Brain Imaging Centre, Montréal Neurological Institute, McGill University, Montréal, Canada
    Contribution
    Data curation, Writing - review and editing
    Competing interests
    No competing interests declared
  5. Ben D Fulcher

    School of Physics, The University of Sydney, Sydney, Australia
    Contribution
    Resources, Software, Methodology, Writing - review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-3003-4055
  6. Bratislav Misic

    McConnell Brain Imaging Centre, Montréal Neurological Institute, McGill University, Montréal, Canada
    Contribution
    Conceptualization, Supervision, Funding acquisition, Validation, Investigation, Visualization, Methodology, Writing - original draft, Project administration, Writing - review and editing
    For correspondence
    bratislav.misic@mcgill.ca
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-0307-2862

Funding

Natural Sciences and Engineering Research Council of Canada

  • Golia Shafiei

Natural Sciences and Engineering Research Council of Canada (NSERC Discovery Grant RGPIN #017-04265)

  • Bratislav Misic

Canada First Research Excellence Fund

  • Bratislav Misic

Canada Research Chairs

  • Bratislav Misic

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

We thank Vincent Bazinet, Justine Hansen, Estefany Suarez, Bertha Vazquez-Rodriguez and Zhen-Qi Liu for helpful comments and stimulating discussion. This research was undertaken thanks in part to funding from the Canada First Research Excellence Fund, awarded to McGill University for the Healthy Brains for Healthy Lives initiative. BM acknowledges support from the Natural Sciences and Engineering Research Council of Canada (NSERC Discovery Grant RGPIN #017–04265) and from the Canada Research Chairs Program. GS acknowledges support from the Natural Sciences and Engineering Research Council of Canada (NSERC).

Ethics

Human subjects: Informed consent and consent to publish were obtained during data acquisition process (all data used in this study were obtained from publicly available datasets).

Senior Editor

  1. Chris I Baker, National Institute of Mental Health, National Institutes of Health, United States

Reviewing Editor

  1. Lucina Q Uddin, University of Miami, United States

Reviewers

  1. Lucina Q Uddin, University of Miami, United States
  2. Maxwell Bertolero, University of Pennsylvania, United States

Publication history

  1. Received: August 14, 2020
  2. Accepted: December 16, 2020
  3. Accepted Manuscript published: December 17, 2020 (version 1)
  4. Version of Record published: December 29, 2020 (version 2)

Copyright

© 2020, Shafiei et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

  • 1,950
    Page views
  • 286
    Downloads
  • 7
    Citations

Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Download citations (links to download the citations from this article in formats compatible with various reference manager tools)

Open citations (links to open the citations from this article in various online reference manager services)

Further reading

    1. Neuroscience
    Debora Fusca, Peter Kloppenburg
    Research Article

    Local interneurons (LNs) mediate complex interactions within the antennal lobe, the primary olfactory system of insects, and the functional analog of the vertebrate olfactory bulb. In the cockroach Periplaneta americana, as in other insects, several types of LNs with distinctive physiological and morphological properties can be defined. Here, we combined whole-cell patch-clamp recordings and Ca2+ imaging of individual LNs to analyze the role of spiking and nonspiking LNs in inter- and intraglomerular signaling during olfactory information processing. Spiking GABAergic LNs reacted to odorant stimulation with a uniform rise in [Ca2+]i in the ramifications of all innervated glomeruli. In contrast, in nonspiking LNs, glomerular Ca2+ signals were odorant specific and varied between glomeruli, resulting in distinct, glomerulus-specific tuning curves. The cell type-specific differences in Ca2+ dynamics support the idea that spiking LNs play a primary role in interglomerular signaling, while they assign nonspiking LNs an essential role in intraglomerular signaling.

    1. Neuroscience
    Wanhui Sheng et al.
    Research Article Updated

    Hypothalamic oxytocinergic magnocellular neurons have a fascinating ability to release peptide from both their axon terminals and from their dendrites. Existing data indicates that the relationship between somatic activity and dendritic release is not constant, but the mechanisms through which this relationship can be modulated are not completely understood. Here, we use a combination of electrical and optical recording techniques to quantify activity-induced calcium influx in proximal vs. distal dendrites of oxytocinergic magnocellular neurons located in the paraventricular nucleus of the hypothalamus (OT-MCNs). Results reveal that the dendrites of OT-MCNs are weak conductors of somatic voltage changes; however, activity-induced dendritic calcium influx can be robustly regulated by both osmosensitive and non-osmosensitive ion channels located along the dendritic membrane. Overall, this study reveals that dendritic conductivity is a dynamic and endogenously regulated feature of OT-MCNs that is likely to have substantial functional impact on central oxytocin release.