Anatomical and Functional Gradients Shape Dynamic Functional Connectivity in the Human Brain

Large-scale biophysical circuit models can provide mechanistic insights into the fundamental micro-scale and macro-scale properties of brain organization that shape complex patterns of spontaneous brain activity. By allowing local synaptic properties to vary across brain regions, recent large-scale circuit models have demonstrated better fit to empirical observations, such as inter-regional synchrony averaged over several minutes, i.e. static functional connectivity (FC). However, most previous models do not capture how inter-regional synchrony patterns vary over timescales of seconds, i.e., time-varying FC dynamics. Here we developed a spatially-heterogeneous large-scale dynamical circuit model that allowed for variation in local circuit properties across the human cortex. We showed that parameterizing local circuit properties with both anatomical and functional gradients was necessary for generating realistic static and dynamical properties of resting-state fMRI activity. Furthermore, empirical and simulated FC dynamics demonstrated remarkably similar sharp transitions in FC patterns, suggesting the existence of multiple attractors. We found that time-varying regional fMRI amplitude tracked multi-stability in FC dynamics. Causal manipulation of the large-scale circuit model suggested that sensory-motor regions were a driver of FC dynamics. Finally, the spatial distribution of sensory-motor drivers matched the principal gradient of gene expression that encompassed certain interneuron classes, suggesting that heterogeneity in excitation-inhibition balance might shape multi-stability in FC dynamics.

dynamics of each ROI are driven by four components: (1) recurrent (intra-regional) input, (2) 154 inter-regional inputs, (3) external input (potentially from subcortical relays) and (4) neuronal 155 noise. There are "free" parameters associated with each component. First, a larger recurrent 156 connection strength corresponds to stronger recurrent input current. Second, the inter-157 regional inputs depend on the neural activities of other cortical ROIs and the connectional 158 strength between ROIs. The inter-regional connectional strength is parameterized by the SC 159 matrices, scaled by a global scaling constant . Third, is the external input current. Fourth, 160 the neuronal noise is assumed to be Gaussian with standard deviation . 161 In the current study, the recurrent connectional strength , external input current ,    Figure 2B shows a simulated FCD generated by the pMFM using the best model parameters 182 (from the validation set) using SC from the test set. Both empirical and simulated FCD 183 exhibited red off-diagonal blocks representing recurring FC patterns. Across the 10 best 184 candidate sets, KS distance between empirical and simulated FCD was 0.115 ± 0.031 (mean 185 ± std). Correlation between empirical and simulated static FC was 0.66 ± 0.03. As a 186 reference, the correlation between SC and static FC in the test set was 0.28.   Figure 2C shows the simulated FCD using the MFM parameters from our previous 198 study (Wang et al., 2019 Figure 3C).

219
Furthermore, if recurrent connectional strength , external input current , and noise 220 amplitude were allowed to be spatially heterogeneous across brain regions, but not 221 constrained by T1w/T2w or FC gradient (i.e., non-parametric), then simulations could 222 achieve realistic static FC, but not FCD ( Figure S2D). One reason could be the large number 223 of "free" parameters leading to overfitting in the training set.

224
Finally, instead of fitting to both static FC and FCD in the training set, we also tried 225 fitting only to static FC. Not surprisingly, the resulting model yielded unrealistic functional 226 connectivity dynamics ( Figure S3; KS = 0.88 ± 0.004). On the other hand, correlation 227 between static empirical and simulated static FC was 0.74 ± 0.01, which was only slightly 228 better than when optimizing both static FC and FCD ( Figure 2C). This suggests that the goals 229 of generating realistic static FC and FCD were not necessarily contradictory. 230 Overall, these results suggest the importance of parameterizing recurrent connectional 231 strength , external input current , and noise amplitude with spatial gradients that 232 smoothly varied from sensory-motor to association cortex. Furthermore, T1w/T2w and FC 233 gradient are complementary in the sense that combining the two spatial maps led to more   . If our hypothesis were true, we would expect large regional fMRI signal amplitude 301 during the coherent state and small regional fMRI signal amplitude during the incoherent 302 state.  As shown in Figure 5D, the SW-STD was significantly higher during the coherent state than 342 the incoherent state (p = 2.4e-168). Similar results were obtained for the pMFM simulations 343 ( Figure 5E).

345
Sensory-motor regions drive switching behavior in functional connectivity dynamics 346 In the previous section, we found striking correspondence between the FCD mean correlational spatial map ( Figure 6B).

358
Statistical significance was established using a permutation test (see Methods).

359
Almost all cortical regions were significant after correcting for multiple comparisons (FDR q 360 < 0.05; Figure S5). Across both pMFM simulations and empirically observed data, FCD-STD 361 correlations were the highest in sensory-motor regions and lowest in association cortex.

362
There was strong spatial correspondence between simulated and empirical results (r = 0.87; 363 Figure 6C). We note that the pMFM was optimized to yield realistic FCD with no regard for 364 spatial correspondence, so the high level of spatial correspondence suggests that the pMFM 365 was able to generalize to new unseen properties of FCD.

366
To explore the causal relationship between sensory-motor regions and FCD, we tested  Figure 6B) were then perturbed to increase their amplitude. The perturbation led to 372 the successful transition of the FCD into a more coherent state with higher FCD mean (p = 373 6e-14; Figure 7D). Perturbation of the bottom five FCD-STD regions ( Figure 6B) did not 374 lead to an increase in FCD mean. Figure 7E illustrates example results of the perturbation 375 experiment. Similar results were obtained if we perturbed top 10 and bottom 10 regions.

376
Overall, this suggests that sensory-motor regions were a driver of switching behavior in FCD.   Figure 7D).

423
The recurrent connection strength and noise amplitude were also correlated with 424 the PVALB-SST gene expression map under the spin-test, but not the random-gene-pair tests.

425
This suggests a lack of specificity to PVALB-SST ( Figure 7D). The external input was not 426 correlated with any gene expression pattern. "random gene pair" tested for the specificity of PVALB-SST by randomly sampling pairs of 441 brain-specific genes. P values that survived the false discovery rate (q < 0.05) are bolded.  and default) networks, thus again suggesting potential degeneracy ( Figure S7).

468
In the Schaefer parcellation, time-varying amplitude of sensory-motor time courses 469 tracks switching behavior in time-varying functional connectivity ( Figures S8 and S9).

470
Causal perturbation analysis also confirmed that sensory-motor regions appeared to drive 471 transitions in FCD ( Figure S9). Both simulated and empirical FCD-STD correlation maps 472 were correlated with PVALB-SST gene expression maps (Table S1) (FC) gradients, the resulting mean field model was able to generate dramatically more 504 realistic static FC and FC dynamics than either gradient alone (Figure 3).

505
The optimized mean field model exhibited opposing gradient directions across local 506 circuit parameters (Figure 4). Across all top ten parameter sets, noise amplitude increased from sensory-motor to association cortex, while external input decreased from sensory-motor 508 to association cortex. The higher external input in sensory-motor regions might reflect the  Figure 5A, there are 520 periods of brain activity with strong coherent FC and periods with incoherent FC. We found 521 that the coherent FC state was characterized by larger fMRI signal amplitude across brain 522 regions, while the incoherent FC state was characterized by smaller fMRI signal amplitude 523 ( Figure 5). Intriguingly, transitions in the regional amplitude of sensory-motor regions 524 appeared to track switching behavior in FC dynamics ( Figure 6). Perturbations of the mean 525 field model suggests that this relationship might be causal.

526
Regional fMRI amplitude has been previously linked with the differential expression  Here, we found that PVALB-SST gene expression map correlates with the spatial distribution 533 of sensory-motor drivers whose time-varying amplitude tracks functional connectivity 534 dynamics (Figure 7). 535 However, we note that this association cannot be solely attributed to PVALB-SST  Details of the HCP preprocessing can be found in the HCP S1200 manual. We utilized 561 rs-fMRI data, which had already been projected to fsLR surface space, denoised with ICA-  where , ( ) and denote the average synaptic gating variable, population firing rate 593 and total input current of the -th cortical ROI. The total input current is the superposition 594 of three inputs. The first input, the intra-regional input, is controlled by the recurrent 595 connection strength . The second input, the inter-regional input, is controlled by the SC 596 matrix ( is the SC between regions and ), as well as a global scaling factor . The 597 third input is the external input current , which might include inputs from subcortical relays.  In this study, recurrent connection strength , external input current and noise 621 amplitude were allowed to vary across brain regions, while was kept as a constant.  To minimize the cost function in the training set, we seek to compute an "average" 651 FCD matrix. We note that FCD matrices could not be directly averaged across rs-fMRI runs 652 and participants because there was no temporal correspondence across runs during the 653 resting-state. Because the goal here was to compute the KS distance, we simply averaged the 654 pdfs from the FCD matrices all the runs of all participants within the training set, which we 655 referred to as average FCD pdf. When evaluating KS distance in the validation and test sets, 656 average FCD pdfs were also computed using the same approach.

657
Optimization procedure 659 To optimize the cost function, we considered three algorithms: covariance matrix was iterated 500 times, generating 500 candidate parameter sets. This procedure was repeated 666 10 times, yielding 5000 candidate parameter sets. For each algorithm, the 5000 candidate 667 parameter sets were evaluated in the validation set to obtain top 10 candidate parameter sets.

668
Across the three algorithms, CMA-ES performed the best in the validation set ( Figure S1), so 669 this study focused on CMA-ES.

670
The top 10 candidate parameter sets from CMA-ES were then applied to the HCP test were selected, yielding 300 time segments. Low FCD mean was defined as being less than 693 0.6.

694
Perturbation was applied to the neural signals (synaptic gating variable ) of the top 695 5 regions whose SW-STD correlated with FCD ( Figure 6B). We note that during the value. This was repeated 1000 times yielding a complete null distribution.

720
To test the specificity of PVALB-SST, we performed a random-gene-pair tests. A 721 random pair of genes was selected from the 2413 brain-specific genes (Burt et al., 2018).

722
Gene expression difference between the random gene pairs was computed and correlated with 723 the STD-FCD correlation maps generating a null correlation value. This was repeated 10,000 724 times yielding a complete null distribution.

726
Code and data availability 727 This study followed the institutional review board guidelines of corresponding institutions.

728
The HCP diffusion MRI, rs-fMRI and T1w/T2w data are publicly available 729 (https://www.humanconnectome.org/). The code used in this paper is publicly available at  Figure S2. Comparison between the original pMFM (main text) and (A) constraining recurrent connection strength to be constant across ROIs, (B) constraining external input to be constant across ROIs, (C) constraining noise amplitude to be the same across ROIs, and (D) allowing local circuit parameters to vary independent (i.e., not parameterized by anatomical and/or functional gradients). Across all panels, agreement between simulated and empirical static FC was measured using Pearson's correlation, while disagreement between simulated and empirical FCD was measured using KS distance. Across all analyses, top ten model parameter sets were selected from the validation set and applied to the test set. The error bars correspond to standard error across the 10 parameter sets. Across all four panels, the original pMFM yielded the best results.  Figure 2B). In terms of KS distance, there is a large improvement when optimizing both static FC and FCD (KS = 0.12 versus 0.88). However, when optimizing only static FC, the resulting simulated static FC was only slightly better than the original pMFM (r = 0.74 versus 0.66). This suggests that the goals of generating realistic static FC and FCD were not necessarily contradictory. We note that across all analyses, top ten model parameter sets were selected from the validation set and applied to the test set. The error bars correspond to standard error across the 10 parameter sets.   Table S1. Table of correlations between FCD-STD correlational spatial maps and two gene expression maps: PVALB-SST and first principal component of gene expression (Burt et al., 2018;Anderson et al., 2020b). P values that survived the false discovery rate (q < 0.05) are bolded. Standard deviations reported in the table were obtained by bootstrapping.