Get 20M+ Full-Text Papers For Less Than $1.50/day. Start a 7-Day Trial for You or Your Team.

Learn More →

Inversion of a large-scale circuit model reveals a cortical hierarchy in the dynamic resting human brain

Inversion of a large-scale circuit model reveals a cortical hierarchy in the dynamic resting... Downloaded from https://www.science.org on August 12, 2025 SCIENCE ADVANCES RESEARCH ARTICLE COGNITIVE NEUROSCIENCE Copyright © 2019 The Authors, some rights reserved; Inversion of a large-scale circuit model reveals a cortical exclusive licensee American Association hierarchy in the dynamic resting human brain for the Advancement of Science. No claim to 1 1 1 1 1 2,3 Peng Wang , Ru Kong , Xiaolu Kong , Raphaël Liégeois , Csaba Orban , Gustavo Deco , original U.S. Government 4 1,5,6,7 Martijn P. van den Heuvel , B.T. Thomas Yeo * Works. Distributed under a Creative We considered a large-scale dynamical circuit model of human cerebral cortex with region-specific microscale Commons Attribution properties. The model was inverted using a stochastic optimization approach, yielding markedly better fit to License 4.0 (CC BY). new, out-of-sample resting functional magnetic resonance imaging (fMRI) data. Without assuming the existence of a hierarchy, the estimated model parameters revealed a large-scale cortical gradient. At one end, sensori- motor regions had strong recurrent connections and excitatory subcortical inputs, consistent with localized processing of external stimuli. At the opposing end, default network regions had weak recurrent connections and excitatory subcortical inputs, consistent with their role in internal thought. Furthermore, recurrent connec- tion strength and subcortical inputs provided complementary information for differentiating the levels of the hierarchy, with only the former showing strong associations with other macroscale and microscale proxies of cortical hierarchies (meta-analysis of cognitive functions, principal resting fMRI gradient, myelin, and laminar- specific neuronal density). Overall, this study provides microscale insights into a macroscale cortical hierarchy in the dynamic resting brain. INTRODUCTION by parameters with physical interpretation (e.g., membrane resting There is converging microscale and macroscale evidence of hierar- potentials), some of which could be empirically measured using chical organization in the primate cerebral cortex, with sensory and cellular neurophysiology or histology. Therefore, large-scale circuit association regions at opposite ends of the hierarchy (1). In terms of models can potentially provide insights into the microscale orga- microscale evidence, cortical regions have been ordered into an an- nization of the dynamic resting brain using only macroscale MRI atomical hierarchy based on the laminar patterns of inter-regional measurements. projections (1, 2). Gene expression patterns and magnetic resonance However, most previous large-scale circuit studies assumed that imaging (MRI) estimates of cortical myelin appeared to vary along local circuit properties are the same across brain regions (13–15). this anatomical hierarchy (3). Macroscale evidence has also arisen Since different brain regions have distinct microscale and macroscale from task-evoked functional MRI (fMRI) and lesion studies (4, 5), properties (16), assuming identical parameters across brain regions as well as from resting-state functional connectivity (RSFC) (6). How- is overly simplistic. Here, we relaxed the influential dynamic mean- ever, there are few studies exploring the links between microscale and field model [MFM; (14)], so that the recurrent connection strength macroscale cortical hierarchies. and excitatory subcortical inputs were free to be different across cor- A powerful approach to bridge microscale and macroscale brain tical regions. A previously published optimization framework based organization is the simulation of large-scale biophysical models of on dynamic causal modeling [DCM; (17, 18)] was then adapted to coupled brain regions (7–9). This approach models the local neural automatically optimize the relaxed MFM (rMFM) parameters. The dynamics in cortical regions by parsimonious but biophysically plau- resulting rMFM yielded a markedly better fit to empirical RSFC from sibleneuralmassorneuralfield models (8, 10). The local models new out-of-sample participants (53% improvement over the original are coupled together by anatomical connections estimated from dif- MFM), suggesting the importance of allowing circuit parameters to fusion MRI or invasive tract tracing. The resulting large-scale circuit differ across brain regions. We then investigated the relationships be- models can be used to simulate complex neural dynamics that are tween the rMFM parameters with other potential proxies of cortical transformed into realistic resting-state fMRI (rs-fMRI) signals via hierarchy using an atlas of resting-state networks (19), a meta-analysis an additional biophysical hemodynamic model (11, 12). In con- of cognitive functions (20), the first principal RSFC gradient (6), T1- trast to statistical algorithms [e.g., independent component analysis weighted/T2-weighted (T1w/T2w) MRI estimate of cortical myelin (ICA)] widely used to interrogate macroscale brain organization, the (21), and histological estimates of cell density/size (22). dynamical properties of these large-scale circuit models are governed The contributions of this study are twofold. First, we automatically inferred the parameters of a biologically plausible MFM with hetero- geneous circuit properties. Although the rMFM made no assumption Department of Electrical and Computer Engineering, ASTAR-NUS Clinical Imaging of the existence of a cortical hierarchy, the resulting model parameters Research Centre, Singapore Institute for Neurotechnology and Memory Networks revealed a large-scale cortical gradient. At one end of the cortical gra- Program, National University of Singapore, Singapore, Singapore. Center for Brain and Cognition, Department of Technology and Information, Universitat Pompeu dient, sensorimotor networks had strong recurrent connections and Fabra, Barcelona, Spain. Institució Catalana de la Recerca i Estudis Avançats, Uni- excitatory subcortical inputs. At the other end of the cortical gradient, versitat Barcelona, Barcelona, Spain. Brain Center Rudolf Magnus, Department of the default network had weak recurrent connections and excitatory Psychiatry, University Medical Center Utrecht, Utrecht, Netherlands. Athinoula A. Martinos Center for Biomedical Imaging, Massachusetts General Hospital, Charles- subcortical inputs. This contrasts with studies that have simply as- town, MA, USA. Centre for Cognitive Neuroscience, Duke-NUS Medical School, sumed local circuit parameters to follow an anatomical hierarchy Singapore, Singapore. NUS Graduate School for Integrative Sciences and Engineer- (23, 24) or relied on statistical approaches not explicitly linked to bio- ing, National University of Singapore, Singapore, Singapore. *Corresponding author. Email: [email protected] physical mechanisms (3, 6). Wang et al., Sci. Adv. 2019; 5 : eaat7854 9 January 2019 1of12 Downloaded from https://www.science.org on August 12, 2025 SCIENCE ADVANCES RESEARCH ARTICLE Second, the two rMFM parameters provided complementary in- atory subcortical input, and (iv) neuronal noise. There are “free” formation for differentiating the levels of the cortical hierarchy. At parameters associated with each component. First, w is the strength one end of the hierarchy, recurrent connection strength differentiated of recurrent connections within a region; a higher w increases the among sensorimotor and attention networks and was strongly asso- amount of recurrent input current relative to other inputs. Second, ciated with other macroscale and microscale proxies of cortical hier- the inter-regional inputs depend on the neural activities of the other archies (cognitive functions, principal RSFC gradient, myelin, and cell cortical ROIs and the connection strength between ROIs. The con- density). At the other end of the cortical hierarchy, subcortical inputs nection strength is determined by the corresponding SC entries differentiated among limbic, control, and default networks and was scaled by a global constant G. Thus, a higher G increases the relative only weakly associated with the principal RSFC gradient. Overall, this strength of inter-regional inputs. Third, I is the excitatory subcortical study provides a computational framework for understanding the mi- input (in nanoampere). Fourth, the neuronal noise is assumed to be croscale foundations of macroscale cortical organization. Gaussian with noise amplitude s. While the original MFM assumed all four parameters to be con- stant across brain regions, here, the model was relaxed to allow the RESULTS recurrent connection w and subcortical input I to be different across Automatic optimization of rMFM parameters significantly ROIs, resulting in 138 parameters. We will refer to this extended improves agreement between simulated and empirical RSFC model as rMFM. Given each ROI’s neural activity, the hemodynamic Four hundred fifty-two participants from the Human Connectome model (11, 12) was used to simulate blood oxygen level–dependent Project (HCP) S500 release (25) were randomly divided into training (BOLD) fMRI at each ROI. The simulated fMRI could then be used (n = 226) and test (n = 226) sets. Following previous work (14), the to compute a 68 × 68 FC matrix. Following previous work (14), agree- Desikan-Killiany cortical parcellation (26) with 68 regions of interest ment between the simulated and empirical FC matrices was defined as (ROIs) was used to compute group-averaged RSFC and structural the Pearson’s correlation between the two matrices (ignoring the di- connectivity (SC) matrices from the training and test sets separately. agonal and nonunique entries) and was used as an index of model fit. Figure 1 (A and B) shows the 68 × 68 group-averaged FC and SC A previous optimization framework for inverting neural mass matrices, respectively, from the test set. models for magnetoencephalography (18) was adapted to automat- Neural dynamics of the 68 ROIs were simulated with the MFM ically estimate the rMFM parameters by maximizing the agreement (14). Here, we highlight the intuitions behind the MFM; details of between simulated and empirical FC (see Materials and Methods for the model are found in Materials and Methods. The MFM assumes details). By applying the optimization algorithm to the training set, that the neural dynamics of each ROI is driven by four components: rMFM parameters (w and I in each ROI, as well as G and s) were (i) recurrent (intraregional) input, (ii) inter-regional inputs, (iii) excit- estimated. The estimated model parameters were used together with the test set SC to generate simulated FC (Fig. 1C). Across 1000 simu- FC (empirical) lations, correlation between simulated and empirical FC in the test SC 0.08 0.4 set (Fig. 1D) was, on average, 0.46, with an SD of 0.016. 0.3 The same procedure was repeated for the MFM (i.e., w and I were 0.06 0.2 constrained to be the same across ROIs). Across 1000 simulations, cor- 0.1 relation between simulated and empirical FC in the test set was, on 0.04 average, 0.30, with an SD of 0.016. This is not better than the baseline −0.1 correlation of 0.30 between SC and empirical FC in the test set. 0.02 −0.2 −0.3 0 Sensorimotor systems exhibit strong recurrent connections 68 cortical regions 68 cortical regions and excitatory subcortical input, while the default network exhibits weak recurrent connections CD FC (relaxed MFM) and excitatory subcortical input 0.4 1 r = 0.46 The above results (Fig. 1) underline the importance of accounting 0.3 0.8 for regional heterogeneity of local circuit properties in biophysical 0.2 0.6 modeling. The question is whether the heterogeneity pattern might 0.1 0.4 be biologically plausible or meaningful. To address this question, 0 0.2 Fig. 2 (A and B) shows the spatial distribution of the estimated re- −0.1 0 current connection strength w and excitatory subcortical input I for −0.2 −0.2 the 68 anatomically defined Desikan-Killiany parcels. The black − 0.3 −0.4 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 lines correspond to the boundaries of seven canonical resting-state 68 cortical regions z-transformed FC networks (19). For reference, Fig. 2C shows the seven resting-state (empirical) networks with the same black boundaries as in Fig. 2 (A and B). Fig. 1. Automatic optimization of rMFM parameters yields stronger agreement While the resting-state network boundaries (black boundaries in between empirical and simulated RSFC. (A) 68 × 68 empirical FC matrix of 68 ROIs Fig. 2) do not exactly align with the anatomically defined parcels [see from HCP test set (n =226). (B) 68 × 68 SC matrix from the HCP test set. (C)Simulated further discussion in (19)], there was a notable correspondence be- 68 × 68 FC matrix using SC matrix from the test set and rMFM parameters estimated tween the resting-state networks and the estimated rMFM parameters. from the HCP training set (n = 226). (D) Correlation between inter-region simulated In particular, the sensorimotor regions appeared to exhibit strong re- FC and inter-region empirical FC (ignoring diagonal elements of the matrices). Cor- current connections (Fig. 2A), while regions of the default network relation between SC and empirical FC in the test set was 0.30. Correlation between simulated and empirical FC was 0.46. appeared to exhibit weak subcortical inputs (Fig. 2B). Wang et al., Sci. Adv. 2019; 5 : eaat7854 9 January 2019 2of12 68 cortical regions 68 cortical regions z-transformed FC (relaxed MFM) 68 cortical regions Downloaded from https://www.science.org on August 12, 2025 SCIENCE ADVANCES RESEARCH ARTICLE A B Recurrent connection w Subcortical input I (nA) 0.44 0.82 0.19 0.38 7-Network parcellation Visual Somatomotor (SomMot) Dorsal attention (DorsAttn) Salience/ventral attention (Sal/VentAttn) Limbic Frontoparietal control Default DE Subcortical input I (nA) Recurrent connection w 0.34 0.7 0.3 0.6 0.26 0.5 0.22 0.18 0.4 Fig. 2. Strength of recurrent connections w and subcortical inputs I in 68 anatomically defined ROIs and their relationships with seven resting-state networks. (A) Strength of recurrent connection w in 68 anatomically defined ROIs. (B) Strength of excitatory subcortical input I in 68 anatomically defined ROIs. Parcels correspond to the 68 Desikan-Killiany ROIs (26). Black boundaries correspond to the boundaries of seven canonical resting-state networks (19). (C) Seven resting-state networks (19). (D) Strength of recurrent connections w in the seven resting-state networks. (E) Strength of subcortical input I in the seven resting-state networks. Regions within sensorimotor systems exhibited strong recurrent connections and excitatory subcortical input, while those within the default network exhibited weak recurrent con- nections and excitatory subcortical input. To quantify this phenomenon, we transferred the estimated rMFM Figure 2D shows the strength of recurrent connection w across the parameters from the 68 Desikan-Killiany parcels to the seven resting- resting-state networks. Visual and somatomotor systems exhibited the state networks. Briefly, we performed the “transfer” as follows. First, strongest recurrent connections, while limbic, control, and default we extracted 51 ROIs from the spatially distributed resting-state net- networks exhibited the weakest recurrent connections. Dorsal atten- works. For example, the default network (red in Fig. 2C) was extracted tion and salience/ventral attention networks exhibited intermediate into 10 ROIs across different lobules. Second, we transferred the rMFM recurrent connections. parameter estimates from the 68 Desikan-Killiany parcels to all vertices Figure 2E shows the strength of the excitatory subcortical input I (from the underlying cortical meshes) comprising each anatomical across the resting-state networks. The default network exhibited parcel. Third, we then averaged the parameter estimates across all the weakest subcortical input, while the visual, somatomotor, dorsal vertices constituting each of the 51 resting-state ROIs. attention, salience/ventral attention, and limbic networks exhibited Wang et al., Sci. Adv. 2019; 5 : eaat7854 9 January 2019 3of12 Visual SomMot DorsAttn Sal/VentAttn Limbic Control Default Visual SomMot DorsAttn Sal/VentAttn Limbic Control Default Downloaded from https://www.science.org on August 12, 2025 SCIENCE ADVANCES RESEARCH ARTICLE the strongest subcortical inputs. The control network exhibited inter- mediate subcortical inputs. Regions with high recurrent connection strength are involved in sensorimotor processing, while those with low recurrent connection strength are involved in higher cognitive functions The correspondences between the estimated rMFM parameters and the resting networks (Fig. 2), as well as the well-known correspon- dence between resting and task networks, suggest that the regional variation in recurrent connection strength and subcortical inputs might reflect cortical processing hierarchies, whereby information might flow from sensory regions to association regions. To make this connection more explicit, we considered 12 cognitive components from a previous large-scale meta-analysis of 10,449 task-evoked functional experiments (20). The cognitive components were asso- ciated with distinct patterns of brain activity and task profiles. On the Low High basis of the top 5 tasks associated with each component [table S1; (20)], these components could be interpreted as being involved in Recurrent connection w various neural processes, ranging from sensorimotor functions (e.g., “hand,”“auditory,” and “visual”) to higher cognitive functions (e.g., “working memory” and “internal mentation”). Normalized activation strength To explore the relationships between the estimated rMFM param- eters and the cognitive components, we grouped the 68 Desikan- Killiany ROIs into 10 zones, ranging from low to high recurrent connection strength w (Fig. 3A). The normalized activation strength associated with each cognitive component was averaged within each of the 10 zones (see details in Materials and Methods). Figure 3B shows the average activity of the 10 zones. The cognitive components were ordered on the basis of the average normalized ac- tivation strength within each of the 10 zones. Zones with high recurrent connection strength exhibited greater brain activity associated with sensory perception and motor actions (e.g., “face,” hand, visual, auditory, etc.). On the other hand, zones with low recurrent connection strength were involved in cognitive functions, including working memory, internal mentation, “reward,” “emotion/limbic,” and “inhibition” (Fig. 3B), which were consistent with known cognitive functions of limbic, control, and default net- works (Fig. 2D). The relationship between subcortical input and cognitive functions (fig. S1) was also consistent with the resting-state network analysis (Fig. 2E). More specifically, zones with high subcortical inputs appeared to exhibit a range of cognitive functions (e.g., emotion/limbic, face, hand, visual, “dorsal attention,” etc.) that were consistent with the fact that sensorimotor, attentional, and limbic networks exhibited uniformly high subcortical inputs. On the other hand, zones with low subcortical inputs were involved in working memory and internal men- tation, which are cognitive functions typically associated with the con- trol and default networks. These results (Figs. 2 and 3) suggest that regional variation in Cognitive components rMFM parameters (especially recurrent connection strength) might Fig. 3. Relationship between recurrent connection strength w and BrainMap be more closely related to traditional processing hierarchies. We will cognitive components. (A) 68 Desikan-Killiany ROIs are grouped into 10 zones now consider two other proxies of cortical processing hierarchies: spanning low to high recurrent connection strength w.(B) Twelve cognitive com- the first principal RSFC gradient (6) and cortical myelin content (3). ponents derived from meta-analysis of 10,449 experiments (20) are ordered on the basis of the average normalized activation strength within each of the 10 zones. Strength of recurrent connections and subcortical inputs is Zones with high recurrent connection strength were involved in sensory percep- associated with the first principal connectivity gradient tion and motor actions (visual, auditory, hand, and face), while those with low re- One proxy of cortical processing hierarchy is the first principal RSFC current connection strength were involved in cognitive functions, such as working gradient, obtained by applying a nonlinear dimensionality reduction memory, internal mentation, and reward. Wang et al., Sci. Adv. 2019; 5 : eaat7854 9 January 2019 4of12 Recurrent connection w Downloaded from https://www.science.org on August 12, 2025 SCIENCE ADVANCES RESEARCH ARTICLE A D First principal gradient Myelin 1.1 1.65 −2 2 B E 2 1.7 r = −0.63 1.5 1.6 −9 P = 7.45 × 10 1 1.5 0.5 1.4 0 1.3 -0.5 1.2 r = 0.58 −1 1.1 −7 P = 1.98 × 10 −1.5 1 0.4 0.5 0.6 0.7 0.8 0.9 0.4 0.5 0.6 0.7 0.8 0.9 Recurrent connection w Recurrent connection w C F 2 1.7 1.6 1.5 1 1.5 0.5 1.4 0 1.3 -0.5 1.2 r = −0.42 r = 0.05 −1 1.1 −4 P = 0.68 P = 3.8 × 10 −1.5 0.15 0.2 0.25 0.3 0.35 0.4 0.15 0.2 0.25 0.3 0.35 0.4 Subcortical input I (nA) Subcortical input I (nA) Fig. 4. Associations of estimated rMFM parameters (strength of recurrent connection w and subcortical input I) with the first principal RSFC gradient and relative myelin content. (A) First principal RSFC gradient obtained by diffusion embedding of the human connectome (6). (B) Association between recurrent con- nection w and first principal gradient. (C) Association between subcortical input I and first principal gradient. (D) T1w/T2w ratio map of estimated myelin content (21). (E) Association between recurrent connection w and myelin. (F) Association between subcortical input I and myelin. algorithm to estimate principal components (“gradients”)that the convention of the original study: low principal gradient values in accounted for the greatest RSFC variance within the cerebral cortex sensorimotor regions and high principal gradient values in the default (6). The first principal gradient was interpreted to represent a pro- network. However, it would be equivalent to reverse the direction of cessing hierarchy anchored by sensorimotor regions at one end and the principal gradient so that sensorimotor regions exhibited high the default network at the other end [Fig. 4A; (6)]. principal gradient values, while the default network exhibited low To investigate the relationship between the first principal RSFC principal gradient values. gradient and the estimated rMFM parameters, the first principal gra- dient values were averaged within each Desikan-Killiany ROI and Recurrent connection strength is positively associated with correlated against the estimated recurrent connection strength w estimates of cortical myelin (Fig. 4B) and subcortical input I (Fig. 4C). The principal gradient The T1w/T2w ratio has been widely used as an estimate of relative was strongly correlated with the recurrent connection strength (r = cortical myelin content [Fig. 4A; (21)]. Burt and colleagues (3) sug- −9 −0.63, P = 7.45 × 10 ) and weakly correlated with the subcortical gest that the T1w/T2w myelin estimate is a good macroscale proxy of −4 input (r = −0.42, P =3.8 ×10 ). microscale cortical processing hierarchy. To investigate the relation- We note that the sign of the correlation is arbitrary because the ship between cortical myelin and the estimated rMFM parameters, we direction of the first principal gradient is arbitrary. Here, we followed averaged the relative myelin content within each Desikan-Killiany Wang et al., Sci. Adv. 2019; 5 : eaat7854 9 January 2019 5of12 First principal gradient First principal gradient Myelin Myelin Downloaded from https://www.science.org on August 12, 2025 SCIENCE ADVANCES RESEARCH ARTICLE ROI and correlated against the estimated recurrent connection strength size (Fig. 5). Table 1 reports all examined correlations. All bolded cor- w (Fig. 4B) and subcortical input I (Fig. 4C). Recurrent connection relations survived a false discovery rate of q < 0.05. The recurrent con- strength w was positively correlated with relative myelin content (r = nection strength w was also significantly correlated with neuronal −7 0.58, P =1.98 × 10 ), while there was no association between sub- density in cortical layers 2, 3, and 6. cortical input I and relative myelin content. Replication with a higher-resolution parcellation and other Strength of recurrent connections is positively associated control analyses with increased neuronal density The above analyses were replicated using a higher-resolution Lausanne2008 Neuronal density (neuron count per cubic millimeter) and neuronal parcellation with 114 cortical ROIs (28). More specifically, across size of the Desikan-Killiany parcels were previously estimated by van 1000 simulations, correlation between simulated and empirical FC den Heuvel and colleagues (22) based on the cytoarchitectonic work in the test set was, on average, 0.46, which was significantly better of von Economo and Koskinas (27). To investigate the relationship than the baseline correlation of 0.34 between SC and empirical FC between neuronal density (or size) and estimated rMFM parameters, in the test set. Similar to before, we found that sensorimotor net- we first note that von Economo and Koskinas did not explicitly dif- works exhibit relatively strong recurrent connections and excitatory ferentiate between left and right hemispheres in their mapping. Vi- subcortical input, while the default network exhibits relatively weak sual inspection of Fig. 2 suggests that estimates of the recurrent recurrent connections and excitatory subcortical input (fig. S2). connection strength w and excitatory subcortical input I appeared Just like in the case of the Desikan-Killiany analyses, Lausanne2008 relatively symmetric across the hemispheres. Therefore, the rMFM parcels with low recurrent connection strength were associated with parameter estimates were averaged between corresponding left and higher cognitive functions, while parcels with high recurrent connection right anatomical parcels and then correlated with estimates of neu- strength were associated with sensory perception and motor actions ronal density and neuronal size. (fig. S3). The relationship between subcortical input and cognitive Figure 5 shows that the recurrent connection strength w was functions was less clear (fig. S4). The recurrent connection strength correlated with neuronal density averaged across all cortical layers was also positively correlated with myelin, while both recurrent connec- (r =0.55, P = 0.00071), but not with neuronal size. The excitatory sub- tion strength and subcortical input were negatively correlated with the cortical input I was not correlated with neuronal density or neuronal first principal gradient (fig. S5). A B 16 1200 r = 0.55 r = 0.19 P = 0.00071 P = 0.28 2 200 0.45 0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85 0.45 0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85 Recurrent connection w Recurrent connection w CD 16 1200 r = 0.05 r = 0.031 14 P = 0.78 P = 0.86 2 200 0.15 0.2 0.25 0.3 0.35 0.4 0.15 0.2 0.25 0.3 0.35 0.4 Subcortical input I (nA) Subcortical input I (nA) Fig. 5. Associations between estimated rMFM parameters (strength of recurrent connection w and subcortical input I) and cytoarchitectonic measures (neuronal density and neuronal size) averaged across all cortical layers. (A) Association between recurrent connection w and neuronal density averaged across all cortical layers. (B) Association between recurrent connection w and neuronal size averaged across all cortical layers. (C) Association between subcortical input I and neuronal density averaged across all cortical layers. (D) Association between subcortical input I and neuronal size averaged across all cortical layers. Wang et al., Sci. Adv. 2019; 5 : eaat7854 9 January 2019 6of12 3 3 Neuronal density per mm Neuronal density per mm Neuronal size Neuronal size Downloaded from https://www.science.org on August 12, 2025 SCIENCE ADVANCES RESEARCH ARTICLE the frontoparietal (control) network shifts its connectivity patterns for adaptive execution of different tasks. Menon and Uddin (29) sug- Table 1. Pearson’s correlation between estimated rMFM parameters gested that the salience network is involved in detecting behaviorally (recurrent connection w and subcortical input I) and cytoarchitectonic data (neuronal cell density and cell size). P values that survived a false relevant stimuli and switching between large-scale networks to facil- discovery rate of q < 0.05 are bolded. itate access to attention and working memory resources and initiate motor systems. Margulies and colleagues (6) have proposed that the wP I P default network is at one end of the cortical processing hierarchy, Layer 1 density −0.13 0.48 0.34 0.053 allowing it to process abstract information not related to direct sen- Layer 2 density 0.50 0.0038 0.23 0.21 sory stimuli. Our results provide evidence for a large-scale hierarchy in the hu- Layer 3 density 0.52 0.0015 0.078 0.66 man cerebral cortex, with the default network and sensorimotor Layer 4 density 0.39 0.036 −0.10 0.60 systems occupying opposing ends. Although the hierarchy is consistent with a recent study (6), our large-scale computational modeling Layer 5 density −0.11 0.54 0.24 0.17 framework enabled further elaborations and insights at the microscale Layer 6 density 0.56 0.00054 0.23 0.20 level. More specifically, regional variation in recurrent connection strength was associated with all hierarchical proxies, including cogni- Cell density averaged across all layers 0.55 0.00071 0.050 0.78 tive functions, first RSFC principal gradient, myelin, and cell density. Layer 1 size −0.29 0.090 0.026 0.88 On the other hand, regional variation in subcortical inputs was only weakly associated with the first principal RSFC gradient. Layer 2 size −0.26 0.15 −0.17 0.37 The recurrent connection strength and subcortical inputs were effec- Layer 3 size 0.20 0.25 −0.042 0.81 tive in differentiating opposite ends of the cortical hierarchy. More spe- cifically, recurrent connection strength strongly differentiated among Layer 4 size 0.31 0.10 −0.052 0.79 visual, somatomotor, dorsal attention, and salience/ventral attention Layer 5 size 0.26 0.14 0.13 0.46 networks (Fig.2D),while limbic,control,and default networkswere Layer 6 size −0.18 0.30 −0.18 0.31 indistinguishable. Previously published hierarchical proxies (e.g., mye- lin and first principal RSFC gradient) exhibited the same issue [see Cell size averaged across all layers 0.19 0.28 0.031 0.86 Figure 3B of Margulies et al.(6)], which is not unexpected, given their strong associations with recurrent connection strength. On the other hand, the subcortical inputs strongly differentiated among limbic, Recurrent connection strength w was correlated with neuronal control, and default networks (Fig. 2E). Thus, the two rMFM para- density averaged across all cortical layers (r =0.49, P = 0.00012), but meters (recurrent connection strength and subcortical inputs) were not with neuronal size averaged across all cortical layers (table S2). complementary in differentiating the levels of the cortical hierarchy. The excitatory subcortical input I was not correlated with neuronal density or neuronal size. The recurrent connection strength w was Large-scale gradients of recurrent connection strength and also significantly correlated with neuronal density in cortical layers subcortical inputs 2, 4, and 6, as well as neuronal size in cortical layer 4 (table S2). At one end of the hierarchy, regions within sensorimotor resting-state Last, we considered an alternative division of the 68 Desikan- networks exhibited strong recurrent connections (Fig. 2). This rela- Killiany ROIs into five cortical types (27). These cortical types were tionship with brain function was made more explicit by our meta- originally defined on the basis of laminar cytoarchitectonic patterns analysis, showing that regions involved in sensorimotor functions ex- and have been functionally interpreted along an action-perception hibited strong recurrent connections (Fig. 3). Furthermore, regions axis, with type 1 representing action and type 5 representing perception. with strong recurrent connection strength also exhibited low RSFC Cortical types 1 to 3 appeared to exhibit low recurrent connection principal gradient values and greater myelin content (Fig. 4). The strength w and subcortical input I (fig. S6). strong positive association between the strength of recurrent connec- tions and neuronal cell density (Fig. 5) suggests an underlying cellular basis for this phenomenon. Previous studies have demonstrated that DISCUSSION brain regions with densely packed neurons are associated with There is a rich literature on hierarchical organization within the pri- specialized local processing (33, 34). Consistent with this interpreta- mate cerebral cortex (1, 29, 30). While traditional accounts of cortical tion, regions within sensorimotor resting-state networks also hierarchies have suggested that sensory streams terminate in pre- appeared to have strong excitatory subcortical inputs (Fig. 2), which frontal areas (5, 31), others have suggested that the apex of processing might correspond to the flow of sensory information from the exter- hierarchy is occupied by parallel, distributed association networks nal environment via subcortical relays. spanning frontal, parietal, temporal, and cingulate cortex (32). Under At the other end of the cortical hierarchy, regions within the this latter account, cortical areas within each association network default network had weak excitatory subcortical inputs (Fig. 2), sug- might occupy the same hierarchical level. gesting the lack of a direct flow of information from the external Recent human brain network literature has largely supported this milieu. This supports the view that the default network is instead network perspective. However, it also proposes that certain associa- more associated with internal processing, as suggested by multiple tion networks might “control” other networks, suggesting (explicitly studies on its involvement in self-generated thought, such as auto- or implicitly) hierarchical differentiation among association net- biographical memory, mind wandering, and thinking about the fu- works (29, 30). For example, Cole and colleagues (30) suggested that ture (35). We also found the default network to have weak recurrent Wang et al., Sci. Adv. 2019; 5 : eaat7854 9 January 2019 7of12 Downloaded from https://www.science.org on August 12, 2025 SCIENCE ADVANCES RESEARCH ARTICLE connections (Fig. 2). Furthermore, regions involved in internal men- very different goal of estimating effective connectivity). This mean- tation exhibited weak recurrent connections (Fig. 3). If strong recur- field reduction is mathematically equivalent to the statistical mechan- rent connections are important for specialized local processing, then ics approach, in which atoms are modeled as populations of atoms, weak recurrent connections might be consistent with the default net- rather than as individual atoms. Similarly, the complex behavior of work’s putative role as a hub of transmodal information integration (6). individual neurons can be mathematically reduced to the behaviors The attentional, limbic, and control networks occupy intermediate of populations of neurons. Obviously, some information has been lost zones in the cortical hierarchy but are further differentiated when both in this reduction, but importantly, a more realistic model of individual recurrent connections and subcortical inputs were considered. For ex- neurons is not better for answering all neuroscientific questions. For ample, the dorsal attention and salience/ventral attention networks example, one could faithfully simulate a single cortical column with exhibited an intermediate recurrent connection strength in contrast thousands of parameters (39), but such a model cannot be easily re- to the limbic and control networks, which exhibited weak recurrent lated to macroscale observations in fMRI (with current approaches). connection strength similar to the default network (Fig. 2). On the In this study, the application of the DCM algorithm to the MFM is other hand, the control network exhibited an intermediate level of an advance along the biological realism spectrum. Nevertheless, it is subcortical inputs, in contrast to the limbic and attentional networks, clear that the specific MFM (14) used in this paper is still a significant which exhibited a strong level of subcortical inputs similar to the early simplification of actual biological processes. For example, the current sensorimotor systems. Thus, we might interpret the two attentional MFM does not account for recurrent excitatory and inhibitory connec- networks as sensorimotor processing systems situated between the tions within the same region (10, 23). As such, the recurrent connec- early sensory (and/or late motor) systems and upstream limbic, control, tion strength in the rMFM can only be thought of as the superposition and default networks. This is consistent with the role of the salience/ of both local excitatory and inhibitory effects. A clear future extension ventral attention network in bottom-up processing of sensory informa- is to consider an rMFM that incorporated both recurrent excitatory tion (36). In the case of the dorsal attention network, the regions com- and inhibitory connections. prising the network include medial temporal complex (MT+), lateral intraparietal (LIP) region, and the frontal eye fields (FEFs). These re- Model fit gions have also been considered as part of a canonical sensorimotor Although the agreement between empirical and simulated FC mark- pathway extending from MT+ to LIP and FEF (1, 19). edly improved from the SF-FC baseline correlation of 0.30 to 0.46, it suggests possible room for improvement from better modeling (pre- Model realism vious section), anatomical connectivity estimates, and/or optimization Statistical approaches (e.g., k-means,ICA,temporalICA,sliding procedures. Given the highly nonlinear nature of the rMFM, we can- window correlations, and hidden Markov models) have been widely not rule out the possibility of other local optima that might yield a used to study brain network organization and dynamics (37). While better fit between empirical and simulated FC. these models have provided significant insights into the human brain, It is difficult to directly compare our results with the quality of fit in they are not meant to “mimic” actual brain mechanisms. On the other the literature on simulating FC from SC. Unfortunately, many papers hand, biophysical models (e.g., spiking neural networks and neural do not report baseline SC-FC correlations, so it is difficult to evaluate mass models) are mechanistically plausible by incorporating model how much “work” the modeling performs beyond the SC-FC correla- parameters linked to microscale characteristics (e.g., membrane rest- tion baseline. For example, Honey and colleagues (13) achieved an ing potentials), some of which could be empirically measured using agreement of 0.70 between empirical and simulated FC. However, cellular neurophysiology or histology (8). Thus, biophysically plausi- they also reported a baseline SC-FC correlation of 0.66. Thus, our ble models can potentially be very useful for linking microscale circuit- modeling results were arguably better, although the correlation be- level mechanisms with macroscale approaches (e.g., fMRI) widely used tween empirical and simulated FC was lower in our study. to study the human brain. Furthermore, most papers reported (in-sample) correlations using To illustrate the difference between statistical and biophysical the same data used to develop their models. Since these papers used models, we note that, although our hierarchy is consistent with the significantly less parameters, overfitting might be less of an issue than first RSFC principal gradient (6), the direction of the principal gradi- in our study. Nevertheless, the true (out-of-sample) correlations be- ent is arbitrary. In Fig. 3A, we followed the convention of the original tween empirical FC and simulated FC are almost certainly going to study: low principal gradient values in sensorimotor regions and high be less than the reported correlations. values in the default network. However, it would be equivalent to have high principal gradient values in sensorimotor regions and low prin- cipal gradient values in the default network. By contrast, we cannot MATERIALS AND METHODS simply invert the recurrent connection strength values and conclude Data that sensorimotor regions have low recurrent connection strength We considered 452 participants from the HCP S500 release (25). All and association regions have high recurrent connection strength. Thus, imaging data were collected on a custom-made Siemens 3T Skyra the rMFM has conceptual advantages over statistical models. scanner using a multiband sequence. Each participant went through It is worthwhile mentioning that biophysical models lie on a two fMRI sessions on two consecutive days. Two rs-fMRI runs were spectrum of model realism. On this spectrum, the MFM is obviously collected in each session. Each fMRI run was acquired in 2-mm iso- less realistic than detailed spiking models but arguably more realistic tropic resolution with a repetition time (TR) of 0.72 s and lasted for than fMRI-DCM (17, 38) in that the MFM is obtained by a mathe- 14 min and 33 s. The structural data consisted of one 0.7-mm iso- matical (mean-field) reduction of a detailed microscopic model of tropic scan for each subject. The diffusion imaging consisted of six neuronal activity, based on integration of higher-order neurons and runs, each lasting approximately 9 min and 50 s. Diffusion weighting realistic synaptic dynamics (keeping in mind that fMRI-DCM has a consisted of three shells of b = 1000, 2000, and 3000 s/mm , with an Wang et al., Sci. Adv. 2019; 5 : eaat7854 9 January 2019 8of12 Downloaded from https://www.science.org on August 12, 2025 SCIENCE ADVANCES RESEARCH ARTICLE approximately equal number of weighting directions on each shell. be a = 270 n/C, b = 108 Hz, and d = 0.154 s. The kinetic parameters Details of the data collection can be found elsewhere (25). The HCP for synaptic activity were set to be r = 0.641 and t =0.1 s. v (t)is s i participants were randomly divided into training (n = 226) and test uncorrelated standard Gaussian noise, and the noise amplitude is (n = 226) sets. controlled by s. The MFM captures the average neural dynamic behaviors of cor- Preprocessing tical regions at the level of neuronal populations using interpretable Details of the HCP anatomical, functional and diffusion preprocessing dynamic variables and physiological parameters, such as population can be found elsewhere (HCP S500 manual). We used the rs-fMRI firing rate and average synaptic gating variable. In comparison to data that the HCP has already projected to the fsLR surface space, more detailed spiking neuronal networks coupled by AMPA, NMDA smoothed by 2 mm, and denoised with ICA-FIX. Following previous (N-methyl-D-aspartate), and GABA (g-aminobutyric acid) synapses, work on large-scale whole-brain computational modeling (7, 13), the the MFM allows a comprehensive study of the relationship between mean cortical signal was further regressed from the fMRI data. For the model parameters and large-scale brain dynamics with a relatively each rs-fMRI run of each participant, the Pearson’s correlation be- low parametric complexity (7, 14). tween the mean fMRI time series of each parcel and every other parcel The simulated neural activities S arefed to theBalloon-Windkessel was computed, resulting in a 68 × 68 FC matrix. The FC matrices were hemodynamic model (17) to simulate the BOLD signals for each ROI. averaged across runs and then averaged across subjects within the Briefly, synaptic activity S in each ROI causes an increase in vasodi- training (or test) set. latory signal z . Inflow f responds in proportion to this signal with i i In the case of the diffusion data, white matter pathways were recon- concomitant changes in blood volume v and deoxyhemoglobin con- structed using generalized Q-sampling imaging, allowing for the re- tent q . The equations relating these biological processes are as follows construction of complex diffusion fiber configurations (i.e., crossing/ kissing fibers) and streamline tractography. The result was a weighted z ¼ S  kz  gð f  1Þ i i i 68 × 68 SC matrix for each subject, where the weight corresponded to the number of tractography streamlines between two regions. Details f ¼ z of the tractography procedure can be found elsewhere (40). To create the group-averaged SC matrix, for each entry of the matrix, if there 1=a tv ¼ f  v i i were less than 50% of the training (or test) subjects with fibers in the hi entry, it was set to zero. Otherwise, the number of streamlines was f 1=f i 1=a1 tq ¼ 1 ð1  rÞ  q v i i averaged across training (or test) subjects with nonzero streamlines. This averaging procedure followed that of previous work (40). Be- cause diffusion tractography was noisy in individual subjects, the where r = 0.34 is the resting oxygen extraction fraction. The kinetic −1 thresholding procedure helped to remove false positives. Further- parameters rate of signal decay k = 0.65 s , rate of elimination g = −1 more, the training (or test) SC matrix was scaled such that its maximal 0.41 s , hemodynamic transit time t = 0.98 s, and Grubb’sexponent a entry was equal to 0.2. = 0.32 followed previous work (14). Given q and v , the BOLD signal is i i given by (11, 12) Dynamic MFM The MFM was derived by mean-field reduction of a detailed spiking BOLD ¼ V k ð1  q Þþ k 1  þ k ð1  v Þ neuronal network model (14) within each brain region to the following i 0 1 2 3 i set of nonlinear stochastic differential equations : S where V = 0.02 is the resting blood volume fraction, and (k , k , k )are 0 1 2 3 ¼ þ rð1  S ÞHðx Þþ sv ðtÞ i i i s a set of parameters dependent on magnetic field strength and a num- ber of acquisition-dependent parameters. The equations of k are as ax  b Hðx Þ¼ follows (12) 1  expðdðax  bÞÞ k ¼ 4:3ϑ rTE 1 0 x ¼ wJS þ GJ∑C S þ I i i ij j k ¼ er rTE 2 0 where x , H(x ), and S denote the total input current, the popula- i i i tion firing rate, and the average synaptic gating variable at the i-th k ¼ 1  e cortical region, respectively. The total input current x is determined by the recurrent connection strength w, the excitatory subcortical where ϑ =28.265B is the frequency offset at the outer surface of 0 0 input I, and inter-region information flow. C corresponds to the magnetized vessels and depends on the main magnetic field strength ij SC between the i-th and j-th cortical regions, thus controlling the B ,which is 3T in the HCP dataset. At a magnetic field strength of 3T, strength of information flow between the two cortical regions. the intravascular relaxation rate is r = 110 Hz, and the ratio between The global constant G scales the strength of information flow from intravascular and extravascular MR signal is e = 0.47 (12). The echo time other cortical regions to the i-th cortical region, relative to the re- TE = 33.1 ms in the HCP dataset. current connection and subcortical inputs. Following previous work The MFM and hemodynamic model were simulated using (14), the value of synaptic coupling J was set to be 0.2609 nA. Euler’s integration with a time step of 10 ms. The starting values Parameter values for the input-output function H(x)wereset to of S in the MFM were randomly initialized. Simulation length i i Wang et al., Sci. Adv. 2019; 5 : eaat7854 9 January 2019 9of12 Downloaded from https://www.science.org on August 12, 2025 SCIENCE ADVANCES RESEARCH ARTICLE for the BOLD signals was 16.4 min. The first 2 min of the BOLD Step 10: Repeat steps 2 to 9 using different random initializations of signals were discarded, and the time series were downsampled to ϕ. Each random initialization required about 36 hours using 30 cen- 0.72 s to have the same temporal resolution as the empirical BOLD tral processing unit (CPU) cores for the Desikan-Killiany parcella- signals in the HCP. tion. Here, we considered 10 random initializations (in addition to the initialization with the prior mean in step 1). The set of model rMFM and automatic estimation of model parameters parameters with the best Pearson’s correlation between the simulated In the MFM, the recurrent connection strength w and excitatory sub- and empirical FC in the training set was selected. cortical input I were assumed to be the same across regions. Here, we relaxed the model, so that each region i had its own recurrent connec- Resting-state networks tion strength w and excitatory subcortical input I . We refer to the Seven resting-state networks from Yeo and colleagues (19)were i i resulting model as rMFM. considered. Each network was spatially distributed, e.g., the red default We optimized 138 rMFM parameters corresponding to w (for network (Fig. 2C) has separate components in the lateral parietal cortex, each ROI), I (for each ROI), global scaling factor G, and noise co- posterior cingulate/precuneus, lateral temporal cortex, medial pre- efficient s by maximizing the similarity between simulated and empir- frontal cortex, and so on. A connected component analysis was used ical FC in the training set. The optimization was based on a previously to extract all the connected “islands” of each network. Small compo- developed algorithm for inverting neural mass models for MEG (18), nents (less than eight vertices) were reassigned to the most correlated which was, in turn, based on the expectation-maximization algorithm networks. Furthermore, vertices whose neighbors belonged to a differ- used in DCM (17). The algorithm is summarized as follows ent component were removed. This resulted in a final set of 51 parcels. Step 1: Initialize model parameters q =[w , I , G, s] = [0.5, 0.3, It is worth noting that the use of the Desikan-Killiany atlas for the 0 i i 1, 0.001]. MFM followed previous applications of the MFM (7, 14). The resting- Step 2: To ensure that the model parameters q remain positive, state ROIs were not directly used for the MFM because some of the ROIs were small and narrow. As such, the tractography estimates the parameters were reparameterized asϕ ¼ ln . Each parameter would not be reliable for these small ROIs. Furthermore, results were ϕ was assumed to be generated from a prior distribution N(0, 0.25), replicated using the Lausanne2008 parcellation, suggesting that the which effectively translated into a prior mean of [0.5, 0.3, 1, 0.001] for results were robust to the exact choice of parcellation. [w , I , G, s] . We note that this prior distribution was sufficiently weak i i to cover the range of parameters (w, I, G, s) used in the literature. BrainMap cognitive components Step 3: Simulate BOLD signals with model parameter q and com- Cognitive components were estimated from a large-scale meta-analysis pute the FC matrix. The upper triangular entries of the 68 × 68 simu- of the BrainMap database (20). Each component was associated with a lated FC matrix were vectorized into a 2278 × 1 vector y = h(ϕ). Let probability of being recruited by a task (of 83 possible task categories). sim y be the corresponding vector from the empirical FC matrix. The top 5 tasks for each component are found in table S1. Each com- emp Step 4: Compute R = y − y ponent was also associated with the probability of activating a voxel. emp sim ∂hðϕÞ Step 5: Numerically approximate Jacobian J ¼ using Newton’s These components were mapped to fsaverage surface space (20). Be- ∂ϕ difference quotient. cause the probabilities must sum to one over all brain voxels, certain Step 6: Initialize error covariance matrix C = exp(l)Q, where Q cognitive components (e.g., reward) that activated both cortical and is the 2278 × 2278 identity matrix and l = −3. subcortical structures would exhibit lower probabilities on the cortical Step 7: Update l in the error covariance matrix C using R and J surface. Therefore, for each component, the activation probabilities computed from previous steps [c.f. M-step in (17)] were normalized to sum to an arbitrary value of 100,000 over the entire cerebral cortex, which we will refer to as normalized activation strength. 1 1 T 1 T 1 To ensure that the numbers of cognitive components and anatomical P ¼ C  C JðJ C JÞ J C e e e e ROIs were comparable, the 68 Desikan-Killiany ROIs were grouped into 10 zones with increasing recurrent connection strength (such that each g ¼0:5tracefPQgþ 0:5R PQPR zone has roughly the same surface area). For each cognitive component, the average normalized activation strength was computed for each of the H ¼0:5tracefPQPQg 10 zones by summing the normalized activation strength within a zone divided by the surface area of the zone. The average normalized activation iþ1 i l ¼ l þ strength is plotted in Fig. 5B. The 10 zones were ordered from low (blue) to high (red) recurrent connection strength, while the 12 components were i+1 The l was used to update the error covariance matrix C = ordered on the basis of the largest normalized activation strength. For ex- i+1 exp(l )Q, and the updated error covariance matrix was again ample, the visual component has the highest normalized activation used in the above equations to update l until convergence. strength in the red anatomical zone (highest recurrent connection Step 8: Update rMFM parameters ϕ [c.f. E-step in (17)] strength) and was therefore assigned to the rightmost column in Fig. 5B. iþ1 i T 1 1 T 1 1 i Myelin and principal gradients ϕ ¼ ϕ þðJ C J þ C Þ ðJ C R þ C ϕ Þ e q e q The relative cortical myelin analyses used the average T1w/T2w ratio where C = 0.25I , where I is a 68 × 68 identity matrix. map released by the HCP (21). The first principal gradient map was q 68×68 68×68 Step 9: Repeat steps 3 to 8 for 512 iterations. The set of model param- provided by D. Margulies (6). In both cases, no further preprocessing eters q = q logϕ resulting in the best Pearson’s correlation between the was performed other than averaging the myelin and principal gradient simulated and empirical FC in the training set was selected. values within each anatomical ROI. Wang et al., Sci. Adv. 2019; 5 : eaat7854 9 January 2019 10 of 12 Downloaded from https://www.science.org on August 12, 2025 SCIENCE ADVANCES RESEARCH ARTICLE 14. G. Deco, A. Ponce-Alvarez, D. Mantini, G. L. Romani, P. Hagmann, M. Corbetta, Resting-state Code and data availability functional connectivity emerges from structurally and dynamically shaped slow linear This study followed the institutional review board guidelines of fluctuations. J. Neurosci. 33,11239–11252 (2013). corresponding institutions. The code used in this paper is publicly 15. A. Zalesky, A. Fornito, L. Cocchi, L. L. Gollo, M. Breakspear, Time-resolved resting-state available at https://github.com/ThomasYeoLab/CBIG/tree/master/ brain networks. Proc. Natl. Acad. Sci. U.S.A. 111, 10341–10346 (2014). 16. K. Zilles, K. Amunts, Centenary of Brodmann’s map—Conception and fate. Nat. Rev. stable_projects/fMRI_dynamics/Wang2018_MFMem. The HCP Neurosci. 11, 139–145 (2010). diffusion MRI, rs-fMRI, and myelin data are publicly available 17. K. J. Friston, L. Harrison, W. Penny, Dynamic causal modelling. Neuroimage 19, 1273–1302 (https://www.humanconnectome.org/). The 51 resting-state ROIs (2003). (derived from the Yeo2011 resting-state networks) can be found 18. P. Wang, T. R. Knösche, A realistic neural mass model of the cortex with laminar-specific here at https://github.com/ThomasYeoLab/CBIG/tree/master/stable_ connections and synaptic plasticity–evaluation with auditory habituation. PLOS ONE 8, e77876 (2013). projects/brain_parcellation/Yeo2011_fcMRI_clustering. The von Econ- 19. B. T. T. Yeo, F. M. Krienen, J. Sepulcre, M. R. Sabuncu, D. Lashkari, M. Hollinshead, omo data can be obtained by directly contacting M. van den Heuvel. J. L. Roffman, J. W. Smoller, L. Zöllei, J. R. Polimeni, B. Fischl, H. Liu, R. L. Buckner, The The BrainMap cognitive components can be downloaded at https:// organization of the human cerebral cortex estimated by intrinsic functional connectivity. surfer.nmr.mgh.harvard.edu/fswiki/BrainmapOntology_Yeo2015. J. Neurophysiol. 106, 1125–1165 (2011). 20. B. T. Yeo, F. M. Krienen, S. B. Eickhoff, S. N. Yaakub, P. T. Fox, R. L. Buckner, C. L. Asplund, M. W. Chee, Functional Specialization and Flexibility in Human Association Cortex. Cereb. Cortex 25,3654–3672 (2015). SUPPLEMENTARY MATERIALS 21. M. F. Glasser, T. S. Coalson, E. C. Robinson, C. D. Hacker, J. Harwell, E. Yacoub, K. Ugurbil, Supplementary material for this article is available at http://advances.sciencemag.org/cgi/ J. Andersson, C. F. Beckmann, M. Jenkinson, S. M. Smith, D. C. Van Essen, A multi-modal content/full/5/1/eaat7854/DC1 parcellation of human cerebral cortex. Nature 536,171–178 (2016). Fig. S1. Relationship between subcortical input I and BrainMap cognitive components. 22. M. P. van den Heuvel, L. H. Scholtens, L. F. Barrett, C. C. Hilgetag, M. A. de Reus, Bridging Fig. S2. Strength of recurrent connections w and subcortical input I in 114 anatomically cytoarchitectonics and connectomics in human cerebral cortex. J. Neurosci. 35, defined ROIs and their relationships with seven resting-state networks. 13943–13948 (2015). Fig. S3. Relationship between recurrent connection strength w and BrainMap cognitive 23. R. Chaudhuri, K. Knoblauch, M.-A. Gariel, H. Kennedy, X.-J. Wang, A large-scale circuit components. mechanism for hierarchical dynamical processing in the primate cortex. Neuron 88, Fig. S4. Relationship between subcortical input I and BrainMap cognitive components. 419–431 (2015). Fig. S5. Associations of estimated rMFM parameters (using the Lausanne2008 parcellation) 24. M. R. Joglekar, J. F. Mejias, G. R. Yang, X.-J. Wang, Inter-areal balanced amplification with relative myelin content and first principal gradient of the human connectome. enhances signal propagation in a large-scale circuit model of the primate cortex. Fig. S6. Relationships between cortical types and estimated rMFM parameters. Neuron 98, 222–234.e8 (2018). Table S1. Top 5 tasks recruiting 12 cognitive components (20). 25. D. C. Van Essen, S. M. Smith, D. M. Barch, T. E. J. Behrens, E. Yacoub, K. Ugurbil; WU-Minn Table S2. Pearson’s correlation between estimated rMFM parameters (recurrent connection w HCP Consortium, The WU-Minn human connectome project: An overview. Neuroimage and subcortical input I) using the Lausanne2008 parcellation and cytoarchitectonic data 80,62–79 (2013). (neuronal cell density and cell size). 26. R. S. Desikan, F. Ségonne, B. Fischl, B. T. Quinn, B. C. Dickerson, D. Blacker, R. L. Buckner, A. M. Dale, R. P. Maguire, B. T. Hyman, M. S. Albert, R. J. Killiany, An automated labeling system for subdividing the human cerebral cortex on MRI scans into gyral based regions of interest. Neuroimage 31, 968–980 (2006). REFERENCES AND NOTES 27. C. F. von Economo, G. N. Koskinas, Die Cytoarchitektonik Der Hirnrinde Des Erwachsenen 1. D. J. Felleman, D. C. Van Essen, Distributed hierarchical processing in the primate cerebral Menschen. (J. Springer, 1925). cortex. Cereb. Cortex 1,1–47 (1991). 28. P. Hagmann, L. Cammoun, X. Gigandet, R. Meuli, C. J. Honey, V. J. Wedeen, O. Sporns, 2. N. T. Markov, J. Vezoli, P. Chameau, A. Falchier, R. Quilodran, C. Huissoud, C. Lamy, Mapping the structural core of human cerebral cortex. PLOS Biol. 6, e159 (2008). P. Misery, P. Giroud, S. Ullman, P. Barone, C. Dehay, K. Knoblauch, H. Kennedy, Anatomy of 29. V. Menon, L. Q. Uddin, Saliency, switching, attention and control: A network model of hierarchy: Feedforward and feedback pathways in macaque visual cortex. J. Comp. insula function. Brain Struct. Funct. 214, 655–667 (2010). Neurol. 522, 225–259 (2014). 30. M. W. Cole, J. R. Reynolds, J. D. Power, G. Repovs, A. Anticevic, T. S. Braver, Multi-task 3. J. B. Burt, M. Demirtas, W. J. Eckner, N. M. Navejar, J. L. Ji, W. J. Martin, A. Bernacchia, connectivity reveals flexible hubs for adaptive task control. Nat. Neurosci. 16, 1348–1355 A. Anticevic, J. D. Murray, Hierarchy of transcriptomic specialization across human cortex (2013). captured by structural neuroimaging topography. Nat. Neurosci. 21, 1251–1259 (2018). 31. E. G. Jones, T. P. Powell, An anatomical study of converging sensory pathways within the 4. M. A. Goodale, A. D. Milner, Separate visual pathways for perception and action. cerebral cortex of the monkey. Brain 93, 793–820 (1970). Trends Neurosci. 15,20–25 (1992). 32. P. S. Goldman-Rakic, Topography of cognition: Parallel distributed networks in primate 5. D. Badre, M. D’esposito, Is the rostro-caudal axis of the frontal lobe hierarchical? Nat. Rev. Neurosci. 10, 659–669 (2009). association cortex. Annu. Rev. Neurosci. 11, 137–156 (1988). 6. D. S. Margulies, S. S. Ghosh, A. Goulas, M. Falkiewicz, J. M. Huntenburg, G. Langs, G. Bezgin, 33. J. H. Kaas, Why is brain size so important: Design problems and solutions as neocortex gets biggeror smaller. Brain Mind 1,7–23 (2000). S. B. Eickhoff, F. X. Castellanos, M. Petrides, E. Jefferies, J. Smallwood, Situating the 34. C. E. Collins, Variability in neuron densities across the cortical sheet in primates. default-mode network along a principal gradient of macroscale cortical organization. Brain Behav. Evol. 78,37–50 (2011). Proc. Natl. Acad. Sci. U.S.A. 113, 12574–12579 (2016). 35. J. R. Andrews-Hanna, J. Smallwood, R. N. Spreng, The default network and self-generated 7. E. C. Hansen, D. Battaglia, A. Spiegler, G. Deco, V. K. Jirsa, Functional connectivity thought: Component processes, dynamic control, and clinical relevance. Ann. N. Y. dynamics: Modeling the switching behavior of the resting state. Neuroimage 105, Acad. Sci. 1316,29–52 (2014). 525–535 (2015). 36. M. Corbetta, G. L. Shulman, Control of goal-directed and stimulus-driven attention in the 8. M. Breakspear, Dynamic models of large-scale brain activity. Nat. Neurosci. 20, 340–352 brain. Nat. Rev. Neurosci. 3, 201–215 (2002). (2017). 37. M. G. Preti, T. A. Bolton, D. Van De Ville, The dynamic functional connectome: State-of- 9. M. Schirner, A. R. McIntosh, V. Jirsa, G. Deco, P. Ritter, Inferring multi-scale neural the-art and perspectives. Neuroimage 160,41–54 (2017). mechanisms with brain network modelling. eLife 7, e28927 (2018). 38. K. E. Stephan, L. Kasper, L. M. Harrison, J. Daunizeau, H. E. M. den Ouden, M. Breakspear, 10. G. Deco, A. Ponce-Alvarez, P. Hagmann, G. L. Romani, D. Mantini, M. Corbetta, How K. J. Friston, Nonlinear dynamic causal models for fMRI. Neuroimage 42, 649–662 (2008). local excitation–inhibition ratio impacts the whole brain dynamics. J. Neurosci. 39. S. Haeusler, W. Maass, A statistical analysis of information-processing properties of 34, 7886–7898 (2014). lamina-specific cortical microcircuit models. Cereb. Cortex 17, 149–162 (2007). 11. K. E. Stephan, N. Weiskopf, P. M. Drysdale, P. A. Robinson, K. J. Friston, Comparing 40. M. P. Van Den Heuvel, O. Sporns, Rich-club organization of the human connectome. hemodynamic models with DCM. Neuroimage 38, 387–401 (2007). J. Neurosci. 31, 15775–15786 (2011). 12. J. Heinzle, P. J. Koopmans, H. E. M. den Ouden, S. Raman, K. E. Stephan, A hemodynamic model for layered BOLD signals. Neuroimage 125, 556–570 (2016). 13. C. Honey, O. Sporns, L. Cammoun, X. Gigandet, J. P. Thiran, R. Meuli, P. Hagmann, Acknowledgments: We thank D. Margulies for sharing with us the principal connectivity Predicting human resting-state functional connectivity from structural connectivity. gradients, and the reviewers for helpful feedback on the paper. Funding: This work Proc. Natl. Acad. Sci. U.S.A. 106, 2035–2040 (2009). was supported by Singapore MOE Tier 2 (MOE2014-T2-2-016), NUS Strategic Research Wang et al., Sci. Adv. 2019; 5 : eaat7854 9 January 2019 11 of 12 Downloaded from https://www.science.org on August 12, 2025 SCIENCE ADVANCES RESEARCH ARTICLE (DPRT/944/09/14), NUS SOM Aspiration Fund (R185000271720), Singapore NMRC authors contributed to drafting and revising the manuscript. Competing interests: The authors (CBRG/0088/2015), NUS YIA, and the Singapore National Research Foundation (NRF) declare that they have no competing interests. Data and materials availability: All data needed Fellowship (Class of 2017). Our research also used resources provided by the Center for to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Functional Neuroimaging Technologies, P41EB015896, and instruments supported by Materials. Additional data related to this paper may be requested from the authors. 1S10RR023401, 1S10RR019307, and 1S10RR023043 from the Athinoula A. Martinos Center for Biomedical Imaging at the Massachusetts General Hospital. Our computational work for Submitted 5 April 2018 this article was partially performed on resources of the National Supercomputing Centre, Accepted 1 December 2018 Singapore (www.nscc.sg). Data were provided by the HCP, WU-Minn Consortium (Principal Published 9 January 2019 Investigators: D. Van Essen and K. Ugurbil; 1U54MH091657) funded by the 16 NIH institutes 10.1126/sciadv.aat7854 and centers that support the NIH Blueprint for Neuroscience Research and by the McDonnell Center for Systems Neuroscience at Washington University. Author contributions: All Citation: P. Wang, R. Kong, X. Kong, R. Liégeois, C. Orban, G. Deco, M. P. van den Heuvel, authors contributed to the conception and design of the work, as well as the analysis and B.T. Thomas Yeo, Inversion of a large-scale circuit model reveals a cortical hierarchy in the interpretation of the data. P.W., R.K., X.K., and R.L. created new software used in this work. All dynamic resting human brain. Sci. Adv. 5, eaat7854 (2019). Wang et al., Sci. Adv. 2019; 5 : eaat7854 9 January 2019 12 of 12 http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Science Advances Unpaywall

Inversion of a large-scale circuit model reveals a cortical hierarchy in the dynamic resting human brain

Download PDF
12 pages

Loading next page...
 
/lp/unpaywall/inversion-of-a-large-scale-circuit-model-reveals-a-cortical-hierarchy-vcoyvyC4u5

References (42)

Publisher
Unpaywall
ISSN
2375-2548
DOI
10.1126/sciadv.aat7854
Publisher site
See Article on Publisher Site

Abstract

Downloaded from https://www.science.org on August 12, 2025 SCIENCE ADVANCES RESEARCH ARTICLE COGNITIVE NEUROSCIENCE Copyright © 2019 The Authors, some rights reserved; Inversion of a large-scale circuit model reveals a cortical exclusive licensee American Association hierarchy in the dynamic resting human brain for the Advancement of Science. No claim to 1 1 1 1 1 2,3 Peng Wang , Ru Kong , Xiaolu Kong , Raphaël Liégeois , Csaba Orban , Gustavo Deco , original U.S. Government 4 1,5,6,7 Martijn P. van den Heuvel , B.T. Thomas Yeo * Works. Distributed under a Creative We considered a large-scale dynamical circuit model of human cerebral cortex with region-specific microscale Commons Attribution properties. The model was inverted using a stochastic optimization approach, yielding markedly better fit to License 4.0 (CC BY). new, out-of-sample resting functional magnetic resonance imaging (fMRI) data. Without assuming the existence of a hierarchy, the estimated model parameters revealed a large-scale cortical gradient. At one end, sensori- motor regions had strong recurrent connections and excitatory subcortical inputs, consistent with localized processing of external stimuli. At the opposing end, default network regions had weak recurrent connections and excitatory subcortical inputs, consistent with their role in internal thought. Furthermore, recurrent connec- tion strength and subcortical inputs provided complementary information for differentiating the levels of the hierarchy, with only the former showing strong associations with other macroscale and microscale proxies of cortical hierarchies (meta-analysis of cognitive functions, principal resting fMRI gradient, myelin, and laminar- specific neuronal density). Overall, this study provides microscale insights into a macroscale cortical hierarchy in the dynamic resting brain. INTRODUCTION by parameters with physical interpretation (e.g., membrane resting There is converging microscale and macroscale evidence of hierar- potentials), some of which could be empirically measured using chical organization in the primate cerebral cortex, with sensory and cellular neurophysiology or histology. Therefore, large-scale circuit association regions at opposite ends of the hierarchy (1). In terms of models can potentially provide insights into the microscale orga- microscale evidence, cortical regions have been ordered into an an- nization of the dynamic resting brain using only macroscale MRI atomical hierarchy based on the laminar patterns of inter-regional measurements. projections (1, 2). Gene expression patterns and magnetic resonance However, most previous large-scale circuit studies assumed that imaging (MRI) estimates of cortical myelin appeared to vary along local circuit properties are the same across brain regions (13–15). this anatomical hierarchy (3). Macroscale evidence has also arisen Since different brain regions have distinct microscale and macroscale from task-evoked functional MRI (fMRI) and lesion studies (4, 5), properties (16), assuming identical parameters across brain regions as well as from resting-state functional connectivity (RSFC) (6). How- is overly simplistic. Here, we relaxed the influential dynamic mean- ever, there are few studies exploring the links between microscale and field model [MFM; (14)], so that the recurrent connection strength macroscale cortical hierarchies. and excitatory subcortical inputs were free to be different across cor- A powerful approach to bridge microscale and macroscale brain tical regions. A previously published optimization framework based organization is the simulation of large-scale biophysical models of on dynamic causal modeling [DCM; (17, 18)] was then adapted to coupled brain regions (7–9). This approach models the local neural automatically optimize the relaxed MFM (rMFM) parameters. The dynamics in cortical regions by parsimonious but biophysically plau- resulting rMFM yielded a markedly better fit to empirical RSFC from sibleneuralmassorneuralfield models (8, 10). The local models new out-of-sample participants (53% improvement over the original are coupled together by anatomical connections estimated from dif- MFM), suggesting the importance of allowing circuit parameters to fusion MRI or invasive tract tracing. The resulting large-scale circuit differ across brain regions. We then investigated the relationships be- models can be used to simulate complex neural dynamics that are tween the rMFM parameters with other potential proxies of cortical transformed into realistic resting-state fMRI (rs-fMRI) signals via hierarchy using an atlas of resting-state networks (19), a meta-analysis an additional biophysical hemodynamic model (11, 12). In con- of cognitive functions (20), the first principal RSFC gradient (6), T1- trast to statistical algorithms [e.g., independent component analysis weighted/T2-weighted (T1w/T2w) MRI estimate of cortical myelin (ICA)] widely used to interrogate macroscale brain organization, the (21), and histological estimates of cell density/size (22). dynamical properties of these large-scale circuit models are governed The contributions of this study are twofold. First, we automatically inferred the parameters of a biologically plausible MFM with hetero- geneous circuit properties. Although the rMFM made no assumption Department of Electrical and Computer Engineering, ASTAR-NUS Clinical Imaging of the existence of a cortical hierarchy, the resulting model parameters Research Centre, Singapore Institute for Neurotechnology and Memory Networks revealed a large-scale cortical gradient. At one end of the cortical gra- Program, National University of Singapore, Singapore, Singapore. Center for Brain and Cognition, Department of Technology and Information, Universitat Pompeu dient, sensorimotor networks had strong recurrent connections and Fabra, Barcelona, Spain. Institució Catalana de la Recerca i Estudis Avançats, Uni- excitatory subcortical inputs. At the other end of the cortical gradient, versitat Barcelona, Barcelona, Spain. Brain Center Rudolf Magnus, Department of the default network had weak recurrent connections and excitatory Psychiatry, University Medical Center Utrecht, Utrecht, Netherlands. Athinoula A. Martinos Center for Biomedical Imaging, Massachusetts General Hospital, Charles- subcortical inputs. This contrasts with studies that have simply as- town, MA, USA. Centre for Cognitive Neuroscience, Duke-NUS Medical School, sumed local circuit parameters to follow an anatomical hierarchy Singapore, Singapore. NUS Graduate School for Integrative Sciences and Engineer- (23, 24) or relied on statistical approaches not explicitly linked to bio- ing, National University of Singapore, Singapore, Singapore. *Corresponding author. Email: [email protected] physical mechanisms (3, 6). Wang et al., Sci. Adv. 2019; 5 : eaat7854 9 January 2019 1of12 Downloaded from https://www.science.org on August 12, 2025 SCIENCE ADVANCES RESEARCH ARTICLE Second, the two rMFM parameters provided complementary in- atory subcortical input, and (iv) neuronal noise. There are “free” formation for differentiating the levels of the cortical hierarchy. At parameters associated with each component. First, w is the strength one end of the hierarchy, recurrent connection strength differentiated of recurrent connections within a region; a higher w increases the among sensorimotor and attention networks and was strongly asso- amount of recurrent input current relative to other inputs. Second, ciated with other macroscale and microscale proxies of cortical hier- the inter-regional inputs depend on the neural activities of the other archies (cognitive functions, principal RSFC gradient, myelin, and cell cortical ROIs and the connection strength between ROIs. The con- density). At the other end of the cortical hierarchy, subcortical inputs nection strength is determined by the corresponding SC entries differentiated among limbic, control, and default networks and was scaled by a global constant G. Thus, a higher G increases the relative only weakly associated with the principal RSFC gradient. Overall, this strength of inter-regional inputs. Third, I is the excitatory subcortical study provides a computational framework for understanding the mi- input (in nanoampere). Fourth, the neuronal noise is assumed to be croscale foundations of macroscale cortical organization. Gaussian with noise amplitude s. While the original MFM assumed all four parameters to be con- stant across brain regions, here, the model was relaxed to allow the RESULTS recurrent connection w and subcortical input I to be different across Automatic optimization of rMFM parameters significantly ROIs, resulting in 138 parameters. We will refer to this extended improves agreement between simulated and empirical RSFC model as rMFM. Given each ROI’s neural activity, the hemodynamic Four hundred fifty-two participants from the Human Connectome model (11, 12) was used to simulate blood oxygen level–dependent Project (HCP) S500 release (25) were randomly divided into training (BOLD) fMRI at each ROI. The simulated fMRI could then be used (n = 226) and test (n = 226) sets. Following previous work (14), the to compute a 68 × 68 FC matrix. Following previous work (14), agree- Desikan-Killiany cortical parcellation (26) with 68 regions of interest ment between the simulated and empirical FC matrices was defined as (ROIs) was used to compute group-averaged RSFC and structural the Pearson’s correlation between the two matrices (ignoring the di- connectivity (SC) matrices from the training and test sets separately. agonal and nonunique entries) and was used as an index of model fit. Figure 1 (A and B) shows the 68 × 68 group-averaged FC and SC A previous optimization framework for inverting neural mass matrices, respectively, from the test set. models for magnetoencephalography (18) was adapted to automat- Neural dynamics of the 68 ROIs were simulated with the MFM ically estimate the rMFM parameters by maximizing the agreement (14). Here, we highlight the intuitions behind the MFM; details of between simulated and empirical FC (see Materials and Methods for the model are found in Materials and Methods. The MFM assumes details). By applying the optimization algorithm to the training set, that the neural dynamics of each ROI is driven by four components: rMFM parameters (w and I in each ROI, as well as G and s) were (i) recurrent (intraregional) input, (ii) inter-regional inputs, (iii) excit- estimated. The estimated model parameters were used together with the test set SC to generate simulated FC (Fig. 1C). Across 1000 simu- FC (empirical) lations, correlation between simulated and empirical FC in the test SC 0.08 0.4 set (Fig. 1D) was, on average, 0.46, with an SD of 0.016. 0.3 The same procedure was repeated for the MFM (i.e., w and I were 0.06 0.2 constrained to be the same across ROIs). Across 1000 simulations, cor- 0.1 relation between simulated and empirical FC in the test set was, on 0.04 average, 0.30, with an SD of 0.016. This is not better than the baseline −0.1 correlation of 0.30 between SC and empirical FC in the test set. 0.02 −0.2 −0.3 0 Sensorimotor systems exhibit strong recurrent connections 68 cortical regions 68 cortical regions and excitatory subcortical input, while the default network exhibits weak recurrent connections CD FC (relaxed MFM) and excitatory subcortical input 0.4 1 r = 0.46 The above results (Fig. 1) underline the importance of accounting 0.3 0.8 for regional heterogeneity of local circuit properties in biophysical 0.2 0.6 modeling. The question is whether the heterogeneity pattern might 0.1 0.4 be biologically plausible or meaningful. To address this question, 0 0.2 Fig. 2 (A and B) shows the spatial distribution of the estimated re- −0.1 0 current connection strength w and excitatory subcortical input I for −0.2 −0.2 the 68 anatomically defined Desikan-Killiany parcels. The black − 0.3 −0.4 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 lines correspond to the boundaries of seven canonical resting-state 68 cortical regions z-transformed FC networks (19). For reference, Fig. 2C shows the seven resting-state (empirical) networks with the same black boundaries as in Fig. 2 (A and B). Fig. 1. Automatic optimization of rMFM parameters yields stronger agreement While the resting-state network boundaries (black boundaries in between empirical and simulated RSFC. (A) 68 × 68 empirical FC matrix of 68 ROIs Fig. 2) do not exactly align with the anatomically defined parcels [see from HCP test set (n =226). (B) 68 × 68 SC matrix from the HCP test set. (C)Simulated further discussion in (19)], there was a notable correspondence be- 68 × 68 FC matrix using SC matrix from the test set and rMFM parameters estimated tween the resting-state networks and the estimated rMFM parameters. from the HCP training set (n = 226). (D) Correlation between inter-region simulated In particular, the sensorimotor regions appeared to exhibit strong re- FC and inter-region empirical FC (ignoring diagonal elements of the matrices). Cor- current connections (Fig. 2A), while regions of the default network relation between SC and empirical FC in the test set was 0.30. Correlation between simulated and empirical FC was 0.46. appeared to exhibit weak subcortical inputs (Fig. 2B). Wang et al., Sci. Adv. 2019; 5 : eaat7854 9 January 2019 2of12 68 cortical regions 68 cortical regions z-transformed FC (relaxed MFM) 68 cortical regions Downloaded from https://www.science.org on August 12, 2025 SCIENCE ADVANCES RESEARCH ARTICLE A B Recurrent connection w Subcortical input I (nA) 0.44 0.82 0.19 0.38 7-Network parcellation Visual Somatomotor (SomMot) Dorsal attention (DorsAttn) Salience/ventral attention (Sal/VentAttn) Limbic Frontoparietal control Default DE Subcortical input I (nA) Recurrent connection w 0.34 0.7 0.3 0.6 0.26 0.5 0.22 0.18 0.4 Fig. 2. Strength of recurrent connections w and subcortical inputs I in 68 anatomically defined ROIs and their relationships with seven resting-state networks. (A) Strength of recurrent connection w in 68 anatomically defined ROIs. (B) Strength of excitatory subcortical input I in 68 anatomically defined ROIs. Parcels correspond to the 68 Desikan-Killiany ROIs (26). Black boundaries correspond to the boundaries of seven canonical resting-state networks (19). (C) Seven resting-state networks (19). (D) Strength of recurrent connections w in the seven resting-state networks. (E) Strength of subcortical input I in the seven resting-state networks. Regions within sensorimotor systems exhibited strong recurrent connections and excitatory subcortical input, while those within the default network exhibited weak recurrent con- nections and excitatory subcortical input. To quantify this phenomenon, we transferred the estimated rMFM Figure 2D shows the strength of recurrent connection w across the parameters from the 68 Desikan-Killiany parcels to the seven resting- resting-state networks. Visual and somatomotor systems exhibited the state networks. Briefly, we performed the “transfer” as follows. First, strongest recurrent connections, while limbic, control, and default we extracted 51 ROIs from the spatially distributed resting-state net- networks exhibited the weakest recurrent connections. Dorsal atten- works. For example, the default network (red in Fig. 2C) was extracted tion and salience/ventral attention networks exhibited intermediate into 10 ROIs across different lobules. Second, we transferred the rMFM recurrent connections. parameter estimates from the 68 Desikan-Killiany parcels to all vertices Figure 2E shows the strength of the excitatory subcortical input I (from the underlying cortical meshes) comprising each anatomical across the resting-state networks. The default network exhibited parcel. Third, we then averaged the parameter estimates across all the weakest subcortical input, while the visual, somatomotor, dorsal vertices constituting each of the 51 resting-state ROIs. attention, salience/ventral attention, and limbic networks exhibited Wang et al., Sci. Adv. 2019; 5 : eaat7854 9 January 2019 3of12 Visual SomMot DorsAttn Sal/VentAttn Limbic Control Default Visual SomMot DorsAttn Sal/VentAttn Limbic Control Default Downloaded from https://www.science.org on August 12, 2025 SCIENCE ADVANCES RESEARCH ARTICLE the strongest subcortical inputs. The control network exhibited inter- mediate subcortical inputs. Regions with high recurrent connection strength are involved in sensorimotor processing, while those with low recurrent connection strength are involved in higher cognitive functions The correspondences between the estimated rMFM parameters and the resting networks (Fig. 2), as well as the well-known correspon- dence between resting and task networks, suggest that the regional variation in recurrent connection strength and subcortical inputs might reflect cortical processing hierarchies, whereby information might flow from sensory regions to association regions. To make this connection more explicit, we considered 12 cognitive components from a previous large-scale meta-analysis of 10,449 task-evoked functional experiments (20). The cognitive components were asso- ciated with distinct patterns of brain activity and task profiles. On the Low High basis of the top 5 tasks associated with each component [table S1; (20)], these components could be interpreted as being involved in Recurrent connection w various neural processes, ranging from sensorimotor functions (e.g., “hand,”“auditory,” and “visual”) to higher cognitive functions (e.g., “working memory” and “internal mentation”). Normalized activation strength To explore the relationships between the estimated rMFM param- eters and the cognitive components, we grouped the 68 Desikan- Killiany ROIs into 10 zones, ranging from low to high recurrent connection strength w (Fig. 3A). The normalized activation strength associated with each cognitive component was averaged within each of the 10 zones (see details in Materials and Methods). Figure 3B shows the average activity of the 10 zones. The cognitive components were ordered on the basis of the average normalized ac- tivation strength within each of the 10 zones. Zones with high recurrent connection strength exhibited greater brain activity associated with sensory perception and motor actions (e.g., “face,” hand, visual, auditory, etc.). On the other hand, zones with low recurrent connection strength were involved in cognitive functions, including working memory, internal mentation, “reward,” “emotion/limbic,” and “inhibition” (Fig. 3B), which were consistent with known cognitive functions of limbic, control, and default net- works (Fig. 2D). The relationship between subcortical input and cognitive functions (fig. S1) was also consistent with the resting-state network analysis (Fig. 2E). More specifically, zones with high subcortical inputs appeared to exhibit a range of cognitive functions (e.g., emotion/limbic, face, hand, visual, “dorsal attention,” etc.) that were consistent with the fact that sensorimotor, attentional, and limbic networks exhibited uniformly high subcortical inputs. On the other hand, zones with low subcortical inputs were involved in working memory and internal men- tation, which are cognitive functions typically associated with the con- trol and default networks. These results (Figs. 2 and 3) suggest that regional variation in Cognitive components rMFM parameters (especially recurrent connection strength) might Fig. 3. Relationship between recurrent connection strength w and BrainMap be more closely related to traditional processing hierarchies. We will cognitive components. (A) 68 Desikan-Killiany ROIs are grouped into 10 zones now consider two other proxies of cortical processing hierarchies: spanning low to high recurrent connection strength w.(B) Twelve cognitive com- the first principal RSFC gradient (6) and cortical myelin content (3). ponents derived from meta-analysis of 10,449 experiments (20) are ordered on the basis of the average normalized activation strength within each of the 10 zones. Strength of recurrent connections and subcortical inputs is Zones with high recurrent connection strength were involved in sensory percep- associated with the first principal connectivity gradient tion and motor actions (visual, auditory, hand, and face), while those with low re- One proxy of cortical processing hierarchy is the first principal RSFC current connection strength were involved in cognitive functions, such as working gradient, obtained by applying a nonlinear dimensionality reduction memory, internal mentation, and reward. Wang et al., Sci. Adv. 2019; 5 : eaat7854 9 January 2019 4of12 Recurrent connection w Downloaded from https://www.science.org on August 12, 2025 SCIENCE ADVANCES RESEARCH ARTICLE A D First principal gradient Myelin 1.1 1.65 −2 2 B E 2 1.7 r = −0.63 1.5 1.6 −9 P = 7.45 × 10 1 1.5 0.5 1.4 0 1.3 -0.5 1.2 r = 0.58 −1 1.1 −7 P = 1.98 × 10 −1.5 1 0.4 0.5 0.6 0.7 0.8 0.9 0.4 0.5 0.6 0.7 0.8 0.9 Recurrent connection w Recurrent connection w C F 2 1.7 1.6 1.5 1 1.5 0.5 1.4 0 1.3 -0.5 1.2 r = −0.42 r = 0.05 −1 1.1 −4 P = 0.68 P = 3.8 × 10 −1.5 0.15 0.2 0.25 0.3 0.35 0.4 0.15 0.2 0.25 0.3 0.35 0.4 Subcortical input I (nA) Subcortical input I (nA) Fig. 4. Associations of estimated rMFM parameters (strength of recurrent connection w and subcortical input I) with the first principal RSFC gradient and relative myelin content. (A) First principal RSFC gradient obtained by diffusion embedding of the human connectome (6). (B) Association between recurrent con- nection w and first principal gradient. (C) Association between subcortical input I and first principal gradient. (D) T1w/T2w ratio map of estimated myelin content (21). (E) Association between recurrent connection w and myelin. (F) Association between subcortical input I and myelin. algorithm to estimate principal components (“gradients”)that the convention of the original study: low principal gradient values in accounted for the greatest RSFC variance within the cerebral cortex sensorimotor regions and high principal gradient values in the default (6). The first principal gradient was interpreted to represent a pro- network. However, it would be equivalent to reverse the direction of cessing hierarchy anchored by sensorimotor regions at one end and the principal gradient so that sensorimotor regions exhibited high the default network at the other end [Fig. 4A; (6)]. principal gradient values, while the default network exhibited low To investigate the relationship between the first principal RSFC principal gradient values. gradient and the estimated rMFM parameters, the first principal gra- dient values were averaged within each Desikan-Killiany ROI and Recurrent connection strength is positively associated with correlated against the estimated recurrent connection strength w estimates of cortical myelin (Fig. 4B) and subcortical input I (Fig. 4C). The principal gradient The T1w/T2w ratio has been widely used as an estimate of relative was strongly correlated with the recurrent connection strength (r = cortical myelin content [Fig. 4A; (21)]. Burt and colleagues (3) sug- −9 −0.63, P = 7.45 × 10 ) and weakly correlated with the subcortical gest that the T1w/T2w myelin estimate is a good macroscale proxy of −4 input (r = −0.42, P =3.8 ×10 ). microscale cortical processing hierarchy. To investigate the relation- We note that the sign of the correlation is arbitrary because the ship between cortical myelin and the estimated rMFM parameters, we direction of the first principal gradient is arbitrary. Here, we followed averaged the relative myelin content within each Desikan-Killiany Wang et al., Sci. Adv. 2019; 5 : eaat7854 9 January 2019 5of12 First principal gradient First principal gradient Myelin Myelin Downloaded from https://www.science.org on August 12, 2025 SCIENCE ADVANCES RESEARCH ARTICLE ROI and correlated against the estimated recurrent connection strength size (Fig. 5). Table 1 reports all examined correlations. All bolded cor- w (Fig. 4B) and subcortical input I (Fig. 4C). Recurrent connection relations survived a false discovery rate of q < 0.05. The recurrent con- strength w was positively correlated with relative myelin content (r = nection strength w was also significantly correlated with neuronal −7 0.58, P =1.98 × 10 ), while there was no association between sub- density in cortical layers 2, 3, and 6. cortical input I and relative myelin content. Replication with a higher-resolution parcellation and other Strength of recurrent connections is positively associated control analyses with increased neuronal density The above analyses were replicated using a higher-resolution Lausanne2008 Neuronal density (neuron count per cubic millimeter) and neuronal parcellation with 114 cortical ROIs (28). More specifically, across size of the Desikan-Killiany parcels were previously estimated by van 1000 simulations, correlation between simulated and empirical FC den Heuvel and colleagues (22) based on the cytoarchitectonic work in the test set was, on average, 0.46, which was significantly better of von Economo and Koskinas (27). To investigate the relationship than the baseline correlation of 0.34 between SC and empirical FC between neuronal density (or size) and estimated rMFM parameters, in the test set. Similar to before, we found that sensorimotor net- we first note that von Economo and Koskinas did not explicitly dif- works exhibit relatively strong recurrent connections and excitatory ferentiate between left and right hemispheres in their mapping. Vi- subcortical input, while the default network exhibits relatively weak sual inspection of Fig. 2 suggests that estimates of the recurrent recurrent connections and excitatory subcortical input (fig. S2). connection strength w and excitatory subcortical input I appeared Just like in the case of the Desikan-Killiany analyses, Lausanne2008 relatively symmetric across the hemispheres. Therefore, the rMFM parcels with low recurrent connection strength were associated with parameter estimates were averaged between corresponding left and higher cognitive functions, while parcels with high recurrent connection right anatomical parcels and then correlated with estimates of neu- strength were associated with sensory perception and motor actions ronal density and neuronal size. (fig. S3). The relationship between subcortical input and cognitive Figure 5 shows that the recurrent connection strength w was functions was less clear (fig. S4). The recurrent connection strength correlated with neuronal density averaged across all cortical layers was also positively correlated with myelin, while both recurrent connec- (r =0.55, P = 0.00071), but not with neuronal size. The excitatory sub- tion strength and subcortical input were negatively correlated with the cortical input I was not correlated with neuronal density or neuronal first principal gradient (fig. S5). A B 16 1200 r = 0.55 r = 0.19 P = 0.00071 P = 0.28 2 200 0.45 0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85 0.45 0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85 Recurrent connection w Recurrent connection w CD 16 1200 r = 0.05 r = 0.031 14 P = 0.78 P = 0.86 2 200 0.15 0.2 0.25 0.3 0.35 0.4 0.15 0.2 0.25 0.3 0.35 0.4 Subcortical input I (nA) Subcortical input I (nA) Fig. 5. Associations between estimated rMFM parameters (strength of recurrent connection w and subcortical input I) and cytoarchitectonic measures (neuronal density and neuronal size) averaged across all cortical layers. (A) Association between recurrent connection w and neuronal density averaged across all cortical layers. (B) Association between recurrent connection w and neuronal size averaged across all cortical layers. (C) Association between subcortical input I and neuronal density averaged across all cortical layers. (D) Association between subcortical input I and neuronal size averaged across all cortical layers. Wang et al., Sci. Adv. 2019; 5 : eaat7854 9 January 2019 6of12 3 3 Neuronal density per mm Neuronal density per mm Neuronal size Neuronal size Downloaded from https://www.science.org on August 12, 2025 SCIENCE ADVANCES RESEARCH ARTICLE the frontoparietal (control) network shifts its connectivity patterns for adaptive execution of different tasks. Menon and Uddin (29) sug- Table 1. Pearson’s correlation between estimated rMFM parameters gested that the salience network is involved in detecting behaviorally (recurrent connection w and subcortical input I) and cytoarchitectonic data (neuronal cell density and cell size). P values that survived a false relevant stimuli and switching between large-scale networks to facil- discovery rate of q < 0.05 are bolded. itate access to attention and working memory resources and initiate motor systems. Margulies and colleagues (6) have proposed that the wP I P default network is at one end of the cortical processing hierarchy, Layer 1 density −0.13 0.48 0.34 0.053 allowing it to process abstract information not related to direct sen- Layer 2 density 0.50 0.0038 0.23 0.21 sory stimuli. Our results provide evidence for a large-scale hierarchy in the hu- Layer 3 density 0.52 0.0015 0.078 0.66 man cerebral cortex, with the default network and sensorimotor Layer 4 density 0.39 0.036 −0.10 0.60 systems occupying opposing ends. Although the hierarchy is consistent with a recent study (6), our large-scale computational modeling Layer 5 density −0.11 0.54 0.24 0.17 framework enabled further elaborations and insights at the microscale Layer 6 density 0.56 0.00054 0.23 0.20 level. More specifically, regional variation in recurrent connection strength was associated with all hierarchical proxies, including cogni- Cell density averaged across all layers 0.55 0.00071 0.050 0.78 tive functions, first RSFC principal gradient, myelin, and cell density. Layer 1 size −0.29 0.090 0.026 0.88 On the other hand, regional variation in subcortical inputs was only weakly associated with the first principal RSFC gradient. Layer 2 size −0.26 0.15 −0.17 0.37 The recurrent connection strength and subcortical inputs were effec- Layer 3 size 0.20 0.25 −0.042 0.81 tive in differentiating opposite ends of the cortical hierarchy. More spe- cifically, recurrent connection strength strongly differentiated among Layer 4 size 0.31 0.10 −0.052 0.79 visual, somatomotor, dorsal attention, and salience/ventral attention Layer 5 size 0.26 0.14 0.13 0.46 networks (Fig.2D),while limbic,control,and default networkswere Layer 6 size −0.18 0.30 −0.18 0.31 indistinguishable. Previously published hierarchical proxies (e.g., mye- lin and first principal RSFC gradient) exhibited the same issue [see Cell size averaged across all layers 0.19 0.28 0.031 0.86 Figure 3B of Margulies et al.(6)], which is not unexpected, given their strong associations with recurrent connection strength. On the other hand, the subcortical inputs strongly differentiated among limbic, Recurrent connection strength w was correlated with neuronal control, and default networks (Fig. 2E). Thus, the two rMFM para- density averaged across all cortical layers (r =0.49, P = 0.00012), but meters (recurrent connection strength and subcortical inputs) were not with neuronal size averaged across all cortical layers (table S2). complementary in differentiating the levels of the cortical hierarchy. The excitatory subcortical input I was not correlated with neuronal density or neuronal size. The recurrent connection strength w was Large-scale gradients of recurrent connection strength and also significantly correlated with neuronal density in cortical layers subcortical inputs 2, 4, and 6, as well as neuronal size in cortical layer 4 (table S2). At one end of the hierarchy, regions within sensorimotor resting-state Last, we considered an alternative division of the 68 Desikan- networks exhibited strong recurrent connections (Fig. 2). This rela- Killiany ROIs into five cortical types (27). These cortical types were tionship with brain function was made more explicit by our meta- originally defined on the basis of laminar cytoarchitectonic patterns analysis, showing that regions involved in sensorimotor functions ex- and have been functionally interpreted along an action-perception hibited strong recurrent connections (Fig. 3). Furthermore, regions axis, with type 1 representing action and type 5 representing perception. with strong recurrent connection strength also exhibited low RSFC Cortical types 1 to 3 appeared to exhibit low recurrent connection principal gradient values and greater myelin content (Fig. 4). The strength w and subcortical input I (fig. S6). strong positive association between the strength of recurrent connec- tions and neuronal cell density (Fig. 5) suggests an underlying cellular basis for this phenomenon. Previous studies have demonstrated that DISCUSSION brain regions with densely packed neurons are associated with There is a rich literature on hierarchical organization within the pri- specialized local processing (33, 34). Consistent with this interpreta- mate cerebral cortex (1, 29, 30). While traditional accounts of cortical tion, regions within sensorimotor resting-state networks also hierarchies have suggested that sensory streams terminate in pre- appeared to have strong excitatory subcortical inputs (Fig. 2), which frontal areas (5, 31), others have suggested that the apex of processing might correspond to the flow of sensory information from the exter- hierarchy is occupied by parallel, distributed association networks nal environment via subcortical relays. spanning frontal, parietal, temporal, and cingulate cortex (32). Under At the other end of the cortical hierarchy, regions within the this latter account, cortical areas within each association network default network had weak excitatory subcortical inputs (Fig. 2), sug- might occupy the same hierarchical level. gesting the lack of a direct flow of information from the external Recent human brain network literature has largely supported this milieu. This supports the view that the default network is instead network perspective. However, it also proposes that certain associa- more associated with internal processing, as suggested by multiple tion networks might “control” other networks, suggesting (explicitly studies on its involvement in self-generated thought, such as auto- or implicitly) hierarchical differentiation among association net- biographical memory, mind wandering, and thinking about the fu- works (29, 30). For example, Cole and colleagues (30) suggested that ture (35). We also found the default network to have weak recurrent Wang et al., Sci. Adv. 2019; 5 : eaat7854 9 January 2019 7of12 Downloaded from https://www.science.org on August 12, 2025 SCIENCE ADVANCES RESEARCH ARTICLE connections (Fig. 2). Furthermore, regions involved in internal men- very different goal of estimating effective connectivity). This mean- tation exhibited weak recurrent connections (Fig. 3). If strong recur- field reduction is mathematically equivalent to the statistical mechan- rent connections are important for specialized local processing, then ics approach, in which atoms are modeled as populations of atoms, weak recurrent connections might be consistent with the default net- rather than as individual atoms. Similarly, the complex behavior of work’s putative role as a hub of transmodal information integration (6). individual neurons can be mathematically reduced to the behaviors The attentional, limbic, and control networks occupy intermediate of populations of neurons. Obviously, some information has been lost zones in the cortical hierarchy but are further differentiated when both in this reduction, but importantly, a more realistic model of individual recurrent connections and subcortical inputs were considered. For ex- neurons is not better for answering all neuroscientific questions. For ample, the dorsal attention and salience/ventral attention networks example, one could faithfully simulate a single cortical column with exhibited an intermediate recurrent connection strength in contrast thousands of parameters (39), but such a model cannot be easily re- to the limbic and control networks, which exhibited weak recurrent lated to macroscale observations in fMRI (with current approaches). connection strength similar to the default network (Fig. 2). On the In this study, the application of the DCM algorithm to the MFM is other hand, the control network exhibited an intermediate level of an advance along the biological realism spectrum. Nevertheless, it is subcortical inputs, in contrast to the limbic and attentional networks, clear that the specific MFM (14) used in this paper is still a significant which exhibited a strong level of subcortical inputs similar to the early simplification of actual biological processes. For example, the current sensorimotor systems. Thus, we might interpret the two attentional MFM does not account for recurrent excitatory and inhibitory connec- networks as sensorimotor processing systems situated between the tions within the same region (10, 23). As such, the recurrent connec- early sensory (and/or late motor) systems and upstream limbic, control, tion strength in the rMFM can only be thought of as the superposition and default networks. This is consistent with the role of the salience/ of both local excitatory and inhibitory effects. A clear future extension ventral attention network in bottom-up processing of sensory informa- is to consider an rMFM that incorporated both recurrent excitatory tion (36). In the case of the dorsal attention network, the regions com- and inhibitory connections. prising the network include medial temporal complex (MT+), lateral intraparietal (LIP) region, and the frontal eye fields (FEFs). These re- Model fit gions have also been considered as part of a canonical sensorimotor Although the agreement between empirical and simulated FC mark- pathway extending from MT+ to LIP and FEF (1, 19). edly improved from the SF-FC baseline correlation of 0.30 to 0.46, it suggests possible room for improvement from better modeling (pre- Model realism vious section), anatomical connectivity estimates, and/or optimization Statistical approaches (e.g., k-means,ICA,temporalICA,sliding procedures. Given the highly nonlinear nature of the rMFM, we can- window correlations, and hidden Markov models) have been widely not rule out the possibility of other local optima that might yield a used to study brain network organization and dynamics (37). While better fit between empirical and simulated FC. these models have provided significant insights into the human brain, It is difficult to directly compare our results with the quality of fit in they are not meant to “mimic” actual brain mechanisms. On the other the literature on simulating FC from SC. Unfortunately, many papers hand, biophysical models (e.g., spiking neural networks and neural do not report baseline SC-FC correlations, so it is difficult to evaluate mass models) are mechanistically plausible by incorporating model how much “work” the modeling performs beyond the SC-FC correla- parameters linked to microscale characteristics (e.g., membrane rest- tion baseline. For example, Honey and colleagues (13) achieved an ing potentials), some of which could be empirically measured using agreement of 0.70 between empirical and simulated FC. However, cellular neurophysiology or histology (8). Thus, biophysically plausi- they also reported a baseline SC-FC correlation of 0.66. Thus, our ble models can potentially be very useful for linking microscale circuit- modeling results were arguably better, although the correlation be- level mechanisms with macroscale approaches (e.g., fMRI) widely used tween empirical and simulated FC was lower in our study. to study the human brain. Furthermore, most papers reported (in-sample) correlations using To illustrate the difference between statistical and biophysical the same data used to develop their models. Since these papers used models, we note that, although our hierarchy is consistent with the significantly less parameters, overfitting might be less of an issue than first RSFC principal gradient (6), the direction of the principal gradi- in our study. Nevertheless, the true (out-of-sample) correlations be- ent is arbitrary. In Fig. 3A, we followed the convention of the original tween empirical FC and simulated FC are almost certainly going to study: low principal gradient values in sensorimotor regions and high be less than the reported correlations. values in the default network. However, it would be equivalent to have high principal gradient values in sensorimotor regions and low prin- cipal gradient values in the default network. By contrast, we cannot MATERIALS AND METHODS simply invert the recurrent connection strength values and conclude Data that sensorimotor regions have low recurrent connection strength We considered 452 participants from the HCP S500 release (25). All and association regions have high recurrent connection strength. Thus, imaging data were collected on a custom-made Siemens 3T Skyra the rMFM has conceptual advantages over statistical models. scanner using a multiband sequence. Each participant went through It is worthwhile mentioning that biophysical models lie on a two fMRI sessions on two consecutive days. Two rs-fMRI runs were spectrum of model realism. On this spectrum, the MFM is obviously collected in each session. Each fMRI run was acquired in 2-mm iso- less realistic than detailed spiking models but arguably more realistic tropic resolution with a repetition time (TR) of 0.72 s and lasted for than fMRI-DCM (17, 38) in that the MFM is obtained by a mathe- 14 min and 33 s. The structural data consisted of one 0.7-mm iso- matical (mean-field) reduction of a detailed microscopic model of tropic scan for each subject. The diffusion imaging consisted of six neuronal activity, based on integration of higher-order neurons and runs, each lasting approximately 9 min and 50 s. Diffusion weighting realistic synaptic dynamics (keeping in mind that fMRI-DCM has a consisted of three shells of b = 1000, 2000, and 3000 s/mm , with an Wang et al., Sci. Adv. 2019; 5 : eaat7854 9 January 2019 8of12 Downloaded from https://www.science.org on August 12, 2025 SCIENCE ADVANCES RESEARCH ARTICLE approximately equal number of weighting directions on each shell. be a = 270 n/C, b = 108 Hz, and d = 0.154 s. The kinetic parameters Details of the data collection can be found elsewhere (25). The HCP for synaptic activity were set to be r = 0.641 and t =0.1 s. v (t)is s i participants were randomly divided into training (n = 226) and test uncorrelated standard Gaussian noise, and the noise amplitude is (n = 226) sets. controlled by s. The MFM captures the average neural dynamic behaviors of cor- Preprocessing tical regions at the level of neuronal populations using interpretable Details of the HCP anatomical, functional and diffusion preprocessing dynamic variables and physiological parameters, such as population can be found elsewhere (HCP S500 manual). We used the rs-fMRI firing rate and average synaptic gating variable. In comparison to data that the HCP has already projected to the fsLR surface space, more detailed spiking neuronal networks coupled by AMPA, NMDA smoothed by 2 mm, and denoised with ICA-FIX. Following previous (N-methyl-D-aspartate), and GABA (g-aminobutyric acid) synapses, work on large-scale whole-brain computational modeling (7, 13), the the MFM allows a comprehensive study of the relationship between mean cortical signal was further regressed from the fMRI data. For the model parameters and large-scale brain dynamics with a relatively each rs-fMRI run of each participant, the Pearson’s correlation be- low parametric complexity (7, 14). tween the mean fMRI time series of each parcel and every other parcel The simulated neural activities S arefed to theBalloon-Windkessel was computed, resulting in a 68 × 68 FC matrix. The FC matrices were hemodynamic model (17) to simulate the BOLD signals for each ROI. averaged across runs and then averaged across subjects within the Briefly, synaptic activity S in each ROI causes an increase in vasodi- training (or test) set. latory signal z . Inflow f responds in proportion to this signal with i i In the case of the diffusion data, white matter pathways were recon- concomitant changes in blood volume v and deoxyhemoglobin con- structed using generalized Q-sampling imaging, allowing for the re- tent q . The equations relating these biological processes are as follows construction of complex diffusion fiber configurations (i.e., crossing/ kissing fibers) and streamline tractography. The result was a weighted z ¼ S  kz  gð f  1Þ i i i 68 × 68 SC matrix for each subject, where the weight corresponded to the number of tractography streamlines between two regions. Details f ¼ z of the tractography procedure can be found elsewhere (40). To create the group-averaged SC matrix, for each entry of the matrix, if there 1=a tv ¼ f  v i i were less than 50% of the training (or test) subjects with fibers in the hi entry, it was set to zero. Otherwise, the number of streamlines was f 1=f i 1=a1 tq ¼ 1 ð1  rÞ  q v i i averaged across training (or test) subjects with nonzero streamlines. This averaging procedure followed that of previous work (40). Be- cause diffusion tractography was noisy in individual subjects, the where r = 0.34 is the resting oxygen extraction fraction. The kinetic −1 thresholding procedure helped to remove false positives. Further- parameters rate of signal decay k = 0.65 s , rate of elimination g = −1 more, the training (or test) SC matrix was scaled such that its maximal 0.41 s , hemodynamic transit time t = 0.98 s, and Grubb’sexponent a entry was equal to 0.2. = 0.32 followed previous work (14). Given q and v , the BOLD signal is i i given by (11, 12) Dynamic MFM The MFM was derived by mean-field reduction of a detailed spiking BOLD ¼ V k ð1  q Þþ k 1  þ k ð1  v Þ neuronal network model (14) within each brain region to the following i 0 1 2 3 i set of nonlinear stochastic differential equations : S where V = 0.02 is the resting blood volume fraction, and (k , k , k )are 0 1 2 3 ¼ þ rð1  S ÞHðx Þþ sv ðtÞ i i i s a set of parameters dependent on magnetic field strength and a num- ber of acquisition-dependent parameters. The equations of k are as ax  b Hðx Þ¼ follows (12) 1  expðdðax  bÞÞ k ¼ 4:3ϑ rTE 1 0 x ¼ wJS þ GJ∑C S þ I i i ij j k ¼ er rTE 2 0 where x , H(x ), and S denote the total input current, the popula- i i i tion firing rate, and the average synaptic gating variable at the i-th k ¼ 1  e cortical region, respectively. The total input current x is determined by the recurrent connection strength w, the excitatory subcortical where ϑ =28.265B is the frequency offset at the outer surface of 0 0 input I, and inter-region information flow. C corresponds to the magnetized vessels and depends on the main magnetic field strength ij SC between the i-th and j-th cortical regions, thus controlling the B ,which is 3T in the HCP dataset. At a magnetic field strength of 3T, strength of information flow between the two cortical regions. the intravascular relaxation rate is r = 110 Hz, and the ratio between The global constant G scales the strength of information flow from intravascular and extravascular MR signal is e = 0.47 (12). The echo time other cortical regions to the i-th cortical region, relative to the re- TE = 33.1 ms in the HCP dataset. current connection and subcortical inputs. Following previous work The MFM and hemodynamic model were simulated using (14), the value of synaptic coupling J was set to be 0.2609 nA. Euler’s integration with a time step of 10 ms. The starting values Parameter values for the input-output function H(x)wereset to of S in the MFM were randomly initialized. Simulation length i i Wang et al., Sci. Adv. 2019; 5 : eaat7854 9 January 2019 9of12 Downloaded from https://www.science.org on August 12, 2025 SCIENCE ADVANCES RESEARCH ARTICLE for the BOLD signals was 16.4 min. The first 2 min of the BOLD Step 10: Repeat steps 2 to 9 using different random initializations of signals were discarded, and the time series were downsampled to ϕ. Each random initialization required about 36 hours using 30 cen- 0.72 s to have the same temporal resolution as the empirical BOLD tral processing unit (CPU) cores for the Desikan-Killiany parcella- signals in the HCP. tion. Here, we considered 10 random initializations (in addition to the initialization with the prior mean in step 1). The set of model rMFM and automatic estimation of model parameters parameters with the best Pearson’s correlation between the simulated In the MFM, the recurrent connection strength w and excitatory sub- and empirical FC in the training set was selected. cortical input I were assumed to be the same across regions. Here, we relaxed the model, so that each region i had its own recurrent connec- Resting-state networks tion strength w and excitatory subcortical input I . We refer to the Seven resting-state networks from Yeo and colleagues (19)were i i resulting model as rMFM. considered. Each network was spatially distributed, e.g., the red default We optimized 138 rMFM parameters corresponding to w (for network (Fig. 2C) has separate components in the lateral parietal cortex, each ROI), I (for each ROI), global scaling factor G, and noise co- posterior cingulate/precuneus, lateral temporal cortex, medial pre- efficient s by maximizing the similarity between simulated and empir- frontal cortex, and so on. A connected component analysis was used ical FC in the training set. The optimization was based on a previously to extract all the connected “islands” of each network. Small compo- developed algorithm for inverting neural mass models for MEG (18), nents (less than eight vertices) were reassigned to the most correlated which was, in turn, based on the expectation-maximization algorithm networks. Furthermore, vertices whose neighbors belonged to a differ- used in DCM (17). The algorithm is summarized as follows ent component were removed. This resulted in a final set of 51 parcels. Step 1: Initialize model parameters q =[w , I , G, s] = [0.5, 0.3, It is worth noting that the use of the Desikan-Killiany atlas for the 0 i i 1, 0.001]. MFM followed previous applications of the MFM (7, 14). The resting- Step 2: To ensure that the model parameters q remain positive, state ROIs were not directly used for the MFM because some of the ROIs were small and narrow. As such, the tractography estimates the parameters were reparameterized asϕ ¼ ln . Each parameter would not be reliable for these small ROIs. Furthermore, results were ϕ was assumed to be generated from a prior distribution N(0, 0.25), replicated using the Lausanne2008 parcellation, suggesting that the which effectively translated into a prior mean of [0.5, 0.3, 1, 0.001] for results were robust to the exact choice of parcellation. [w , I , G, s] . We note that this prior distribution was sufficiently weak i i to cover the range of parameters (w, I, G, s) used in the literature. BrainMap cognitive components Step 3: Simulate BOLD signals with model parameter q and com- Cognitive components were estimated from a large-scale meta-analysis pute the FC matrix. The upper triangular entries of the 68 × 68 simu- of the BrainMap database (20). Each component was associated with a lated FC matrix were vectorized into a 2278 × 1 vector y = h(ϕ). Let probability of being recruited by a task (of 83 possible task categories). sim y be the corresponding vector from the empirical FC matrix. The top 5 tasks for each component are found in table S1. Each com- emp Step 4: Compute R = y − y ponent was also associated with the probability of activating a voxel. emp sim ∂hðϕÞ Step 5: Numerically approximate Jacobian J ¼ using Newton’s These components were mapped to fsaverage surface space (20). Be- ∂ϕ difference quotient. cause the probabilities must sum to one over all brain voxels, certain Step 6: Initialize error covariance matrix C = exp(l)Q, where Q cognitive components (e.g., reward) that activated both cortical and is the 2278 × 2278 identity matrix and l = −3. subcortical structures would exhibit lower probabilities on the cortical Step 7: Update l in the error covariance matrix C using R and J surface. Therefore, for each component, the activation probabilities computed from previous steps [c.f. M-step in (17)] were normalized to sum to an arbitrary value of 100,000 over the entire cerebral cortex, which we will refer to as normalized activation strength. 1 1 T 1 T 1 To ensure that the numbers of cognitive components and anatomical P ¼ C  C JðJ C JÞ J C e e e e ROIs were comparable, the 68 Desikan-Killiany ROIs were grouped into 10 zones with increasing recurrent connection strength (such that each g ¼0:5tracefPQgþ 0:5R PQPR zone has roughly the same surface area). For each cognitive component, the average normalized activation strength was computed for each of the H ¼0:5tracefPQPQg 10 zones by summing the normalized activation strength within a zone divided by the surface area of the zone. The average normalized activation iþ1 i l ¼ l þ strength is plotted in Fig. 5B. The 10 zones were ordered from low (blue) to high (red) recurrent connection strength, while the 12 components were i+1 The l was used to update the error covariance matrix C = ordered on the basis of the largest normalized activation strength. For ex- i+1 exp(l )Q, and the updated error covariance matrix was again ample, the visual component has the highest normalized activation used in the above equations to update l until convergence. strength in the red anatomical zone (highest recurrent connection Step 8: Update rMFM parameters ϕ [c.f. E-step in (17)] strength) and was therefore assigned to the rightmost column in Fig. 5B. iþ1 i T 1 1 T 1 1 i Myelin and principal gradients ϕ ¼ ϕ þðJ C J þ C Þ ðJ C R þ C ϕ Þ e q e q The relative cortical myelin analyses used the average T1w/T2w ratio where C = 0.25I , where I is a 68 × 68 identity matrix. map released by the HCP (21). The first principal gradient map was q 68×68 68×68 Step 9: Repeat steps 3 to 8 for 512 iterations. The set of model param- provided by D. Margulies (6). In both cases, no further preprocessing eters q = q logϕ resulting in the best Pearson’s correlation between the was performed other than averaging the myelin and principal gradient simulated and empirical FC in the training set was selected. values within each anatomical ROI. Wang et al., Sci. Adv. 2019; 5 : eaat7854 9 January 2019 10 of 12 Downloaded from https://www.science.org on August 12, 2025 SCIENCE ADVANCES RESEARCH ARTICLE 14. G. Deco, A. Ponce-Alvarez, D. Mantini, G. L. Romani, P. Hagmann, M. Corbetta, Resting-state Code and data availability functional connectivity emerges from structurally and dynamically shaped slow linear This study followed the institutional review board guidelines of fluctuations. J. Neurosci. 33,11239–11252 (2013). corresponding institutions. The code used in this paper is publicly 15. A. Zalesky, A. Fornito, L. Cocchi, L. L. Gollo, M. Breakspear, Time-resolved resting-state available at https://github.com/ThomasYeoLab/CBIG/tree/master/ brain networks. Proc. Natl. Acad. Sci. U.S.A. 111, 10341–10346 (2014). 16. K. Zilles, K. Amunts, Centenary of Brodmann’s map—Conception and fate. Nat. Rev. stable_projects/fMRI_dynamics/Wang2018_MFMem. The HCP Neurosci. 11, 139–145 (2010). diffusion MRI, rs-fMRI, and myelin data are publicly available 17. K. J. Friston, L. Harrison, W. Penny, Dynamic causal modelling. Neuroimage 19, 1273–1302 (https://www.humanconnectome.org/). The 51 resting-state ROIs (2003). (derived from the Yeo2011 resting-state networks) can be found 18. P. Wang, T. R. Knösche, A realistic neural mass model of the cortex with laminar-specific here at https://github.com/ThomasYeoLab/CBIG/tree/master/stable_ connections and synaptic plasticity–evaluation with auditory habituation. PLOS ONE 8, e77876 (2013). projects/brain_parcellation/Yeo2011_fcMRI_clustering. The von Econ- 19. B. T. T. Yeo, F. M. Krienen, J. Sepulcre, M. R. Sabuncu, D. Lashkari, M. Hollinshead, omo data can be obtained by directly contacting M. van den Heuvel. J. L. Roffman, J. W. Smoller, L. Zöllei, J. R. Polimeni, B. Fischl, H. Liu, R. L. Buckner, The The BrainMap cognitive components can be downloaded at https:// organization of the human cerebral cortex estimated by intrinsic functional connectivity. surfer.nmr.mgh.harvard.edu/fswiki/BrainmapOntology_Yeo2015. J. Neurophysiol. 106, 1125–1165 (2011). 20. B. T. Yeo, F. M. Krienen, S. B. Eickhoff, S. N. Yaakub, P. T. Fox, R. L. Buckner, C. L. Asplund, M. W. Chee, Functional Specialization and Flexibility in Human Association Cortex. Cereb. Cortex 25,3654–3672 (2015). SUPPLEMENTARY MATERIALS 21. M. F. Glasser, T. S. Coalson, E. C. Robinson, C. D. Hacker, J. Harwell, E. Yacoub, K. Ugurbil, Supplementary material for this article is available at http://advances.sciencemag.org/cgi/ J. Andersson, C. F. Beckmann, M. Jenkinson, S. M. Smith, D. C. Van Essen, A multi-modal content/full/5/1/eaat7854/DC1 parcellation of human cerebral cortex. Nature 536,171–178 (2016). Fig. S1. Relationship between subcortical input I and BrainMap cognitive components. 22. M. P. van den Heuvel, L. H. Scholtens, L. F. Barrett, C. C. Hilgetag, M. A. de Reus, Bridging Fig. S2. Strength of recurrent connections w and subcortical input I in 114 anatomically cytoarchitectonics and connectomics in human cerebral cortex. J. Neurosci. 35, defined ROIs and their relationships with seven resting-state networks. 13943–13948 (2015). Fig. S3. Relationship between recurrent connection strength w and BrainMap cognitive 23. R. Chaudhuri, K. Knoblauch, M.-A. Gariel, H. Kennedy, X.-J. Wang, A large-scale circuit components. mechanism for hierarchical dynamical processing in the primate cortex. Neuron 88, Fig. S4. Relationship between subcortical input I and BrainMap cognitive components. 419–431 (2015). Fig. S5. Associations of estimated rMFM parameters (using the Lausanne2008 parcellation) 24. M. R. Joglekar, J. F. Mejias, G. R. Yang, X.-J. Wang, Inter-areal balanced amplification with relative myelin content and first principal gradient of the human connectome. enhances signal propagation in a large-scale circuit model of the primate cortex. Fig. S6. Relationships between cortical types and estimated rMFM parameters. Neuron 98, 222–234.e8 (2018). Table S1. Top 5 tasks recruiting 12 cognitive components (20). 25. D. C. Van Essen, S. M. Smith, D. M. Barch, T. E. J. Behrens, E. Yacoub, K. Ugurbil; WU-Minn Table S2. Pearson’s correlation between estimated rMFM parameters (recurrent connection w HCP Consortium, The WU-Minn human connectome project: An overview. Neuroimage and subcortical input I) using the Lausanne2008 parcellation and cytoarchitectonic data 80,62–79 (2013). (neuronal cell density and cell size). 26. R. S. Desikan, F. Ségonne, B. Fischl, B. T. Quinn, B. C. Dickerson, D. Blacker, R. L. Buckner, A. M. Dale, R. P. Maguire, B. T. Hyman, M. S. Albert, R. J. Killiany, An automated labeling system for subdividing the human cerebral cortex on MRI scans into gyral based regions of interest. Neuroimage 31, 968–980 (2006). REFERENCES AND NOTES 27. C. F. von Economo, G. N. Koskinas, Die Cytoarchitektonik Der Hirnrinde Des Erwachsenen 1. D. J. Felleman, D. C. Van Essen, Distributed hierarchical processing in the primate cerebral Menschen. (J. Springer, 1925). cortex. Cereb. Cortex 1,1–47 (1991). 28. P. Hagmann, L. Cammoun, X. Gigandet, R. Meuli, C. J. Honey, V. J. Wedeen, O. Sporns, 2. N. T. Markov, J. Vezoli, P. Chameau, A. Falchier, R. Quilodran, C. Huissoud, C. Lamy, Mapping the structural core of human cerebral cortex. PLOS Biol. 6, e159 (2008). P. Misery, P. Giroud, S. Ullman, P. Barone, C. Dehay, K. Knoblauch, H. Kennedy, Anatomy of 29. V. Menon, L. Q. Uddin, Saliency, switching, attention and control: A network model of hierarchy: Feedforward and feedback pathways in macaque visual cortex. J. Comp. insula function. Brain Struct. Funct. 214, 655–667 (2010). Neurol. 522, 225–259 (2014). 30. M. W. Cole, J. R. Reynolds, J. D. Power, G. Repovs, A. Anticevic, T. S. Braver, Multi-task 3. J. B. Burt, M. Demirtas, W. J. Eckner, N. M. Navejar, J. L. Ji, W. J. Martin, A. Bernacchia, connectivity reveals flexible hubs for adaptive task control. Nat. Neurosci. 16, 1348–1355 A. Anticevic, J. D. Murray, Hierarchy of transcriptomic specialization across human cortex (2013). captured by structural neuroimaging topography. Nat. Neurosci. 21, 1251–1259 (2018). 31. E. G. Jones, T. P. Powell, An anatomical study of converging sensory pathways within the 4. M. A. Goodale, A. D. Milner, Separate visual pathways for perception and action. cerebral cortex of the monkey. Brain 93, 793–820 (1970). Trends Neurosci. 15,20–25 (1992). 32. P. S. Goldman-Rakic, Topography of cognition: Parallel distributed networks in primate 5. D. Badre, M. D’esposito, Is the rostro-caudal axis of the frontal lobe hierarchical? Nat. Rev. Neurosci. 10, 659–669 (2009). association cortex. Annu. Rev. Neurosci. 11, 137–156 (1988). 6. D. S. Margulies, S. S. Ghosh, A. Goulas, M. Falkiewicz, J. M. Huntenburg, G. Langs, G. Bezgin, 33. J. H. Kaas, Why is brain size so important: Design problems and solutions as neocortex gets biggeror smaller. Brain Mind 1,7–23 (2000). S. B. Eickhoff, F. X. Castellanos, M. Petrides, E. Jefferies, J. Smallwood, Situating the 34. C. E. Collins, Variability in neuron densities across the cortical sheet in primates. default-mode network along a principal gradient of macroscale cortical organization. Brain Behav. Evol. 78,37–50 (2011). Proc. Natl. Acad. Sci. U.S.A. 113, 12574–12579 (2016). 35. J. R. Andrews-Hanna, J. Smallwood, R. N. Spreng, The default network and self-generated 7. E. C. Hansen, D. Battaglia, A. Spiegler, G. Deco, V. K. Jirsa, Functional connectivity thought: Component processes, dynamic control, and clinical relevance. Ann. N. Y. dynamics: Modeling the switching behavior of the resting state. Neuroimage 105, Acad. Sci. 1316,29–52 (2014). 525–535 (2015). 36. M. Corbetta, G. L. Shulman, Control of goal-directed and stimulus-driven attention in the 8. M. Breakspear, Dynamic models of large-scale brain activity. Nat. Neurosci. 20, 340–352 brain. Nat. Rev. Neurosci. 3, 201–215 (2002). (2017). 37. M. G. Preti, T. A. Bolton, D. Van De Ville, The dynamic functional connectome: State-of- 9. M. Schirner, A. R. McIntosh, V. Jirsa, G. Deco, P. Ritter, Inferring multi-scale neural the-art and perspectives. Neuroimage 160,41–54 (2017). mechanisms with brain network modelling. eLife 7, e28927 (2018). 38. K. E. Stephan, L. Kasper, L. M. Harrison, J. Daunizeau, H. E. M. den Ouden, M. Breakspear, 10. G. Deco, A. Ponce-Alvarez, P. Hagmann, G. L. Romani, D. Mantini, M. Corbetta, How K. J. Friston, Nonlinear dynamic causal models for fMRI. Neuroimage 42, 649–662 (2008). local excitation–inhibition ratio impacts the whole brain dynamics. J. Neurosci. 39. S. Haeusler, W. Maass, A statistical analysis of information-processing properties of 34, 7886–7898 (2014). lamina-specific cortical microcircuit models. Cereb. Cortex 17, 149–162 (2007). 11. K. E. Stephan, N. Weiskopf, P. M. Drysdale, P. A. Robinson, K. J. Friston, Comparing 40. M. P. Van Den Heuvel, O. Sporns, Rich-club organization of the human connectome. hemodynamic models with DCM. Neuroimage 38, 387–401 (2007). J. Neurosci. 31, 15775–15786 (2011). 12. J. Heinzle, P. J. Koopmans, H. E. M. den Ouden, S. Raman, K. E. Stephan, A hemodynamic model for layered BOLD signals. Neuroimage 125, 556–570 (2016). 13. C. Honey, O. Sporns, L. Cammoun, X. Gigandet, J. P. Thiran, R. Meuli, P. Hagmann, Acknowledgments: We thank D. Margulies for sharing with us the principal connectivity Predicting human resting-state functional connectivity from structural connectivity. gradients, and the reviewers for helpful feedback on the paper. Funding: This work Proc. Natl. Acad. Sci. U.S.A. 106, 2035–2040 (2009). was supported by Singapore MOE Tier 2 (MOE2014-T2-2-016), NUS Strategic Research Wang et al., Sci. Adv. 2019; 5 : eaat7854 9 January 2019 11 of 12 Downloaded from https://www.science.org on August 12, 2025 SCIENCE ADVANCES RESEARCH ARTICLE (DPRT/944/09/14), NUS SOM Aspiration Fund (R185000271720), Singapore NMRC authors contributed to drafting and revising the manuscript. Competing interests: The authors (CBRG/0088/2015), NUS YIA, and the Singapore National Research Foundation (NRF) declare that they have no competing interests. Data and materials availability: All data needed Fellowship (Class of 2017). Our research also used resources provided by the Center for to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Functional Neuroimaging Technologies, P41EB015896, and instruments supported by Materials. Additional data related to this paper may be requested from the authors. 1S10RR023401, 1S10RR019307, and 1S10RR023043 from the Athinoula A. Martinos Center for Biomedical Imaging at the Massachusetts General Hospital. Our computational work for Submitted 5 April 2018 this article was partially performed on resources of the National Supercomputing Centre, Accepted 1 December 2018 Singapore (www.nscc.sg). Data were provided by the HCP, WU-Minn Consortium (Principal Published 9 January 2019 Investigators: D. Van Essen and K. Ugurbil; 1U54MH091657) funded by the 16 NIH institutes 10.1126/sciadv.aat7854 and centers that support the NIH Blueprint for Neuroscience Research and by the McDonnell Center for Systems Neuroscience at Washington University. Author contributions: All Citation: P. Wang, R. Kong, X. Kong, R. Liégeois, C. Orban, G. Deco, M. P. van den Heuvel, authors contributed to the conception and design of the work, as well as the analysis and B.T. Thomas Yeo, Inversion of a large-scale circuit model reveals a cortical hierarchy in the interpretation of the data. P.W., R.K., X.K., and R.L. created new software used in this work. All dynamic resting human brain. Sci. Adv. 5, eaat7854 (2019). Wang et al., Sci. Adv. 2019; 5 : eaat7854 9 January 2019 12 of 12

Journal

Science Advances – Unpaywall

Published: Jan 4, 2019

There are no references for this article.