Complexity of intrinsic brain dynamics

7 The brain is a complex, nonlinear system, exhibiting ever-evolving patterns 8 of activities even without external inputs or tasks. Such intrinsic dynamics play 9 a key role in cognitive functions and psychiatric disorders. A challenge is to link 10 the intrinsic dynamics to the underlying structure, given the nonlinearity. Here 11 we use a biophysically constrained, nonlinear-dynamical model to show how 12 the complexity of intrinsic brain dynamics, manifested as its multistability and 13 temporal diversity, can be sculpted by structural properties across scales. At a 14 local level, multistability and temporal diversity can be induced by sufficient 15 recurrent excitatory connectivity and its heterogeneity. At a global level, such 16 functional complexity can also be created by the synergistic interaction between 17 monostable, locally identical regions. Coordination between model brain 18 regions across attractors in the multistable landscape predicts human functional 19 connectivity. Compared to dynamics near a single attractor, cross-attractor 20 coordination better accounts for functional links uncorrelated with structural 21 connectivity. Energy costs of cross-attractor coordination are modulated by 22 both local and global connectivity, and higher in the Default Mode Network. 23 These findings hint that functional connectivity underscores transitions between 24 alternative patterns of activity in the brain—even more than the patterns 25 themselves. This work provides a systematic framework for characterizing 26 intrinsic brain dynamics as a web of cross-attractor transitions and their energy 27 costs. The framework may be used to predict transitions and energy costs 28 associated with experimental or clinical interventions. 29 Significance Statement 30 The brain is a multifunctional system: different brain regions can coordinate flexibly 31 to perform different tasks. Understanding the regional and global structural con32 straints on brain function is critical to understand cognition. Here, using a unified 33 biophysical network model, we show how structural constraints across scales jointly 34 shape the brain’s intrinsic functional repertoire – the set of all possible patterns that 35 a brain could generate. The modeled functional repertoire is enriched by both the 36 flexibility of individual brain regions and their synergistic interaction in a global 37 network. By carefully examining the modeled repertoire, we found that human 38 resting-state functional connectivity was better predicted by transitions between 39 brain activity patterns rather than any specific pattern per se. 40

both local and global connectivity, and higher in the Default Mode Network. 23 These findings hint that functional connectivity underscores transitions between 24 alternative patterns of activity in the brain-even more than the patterns 25 themselves. This work provides a systematic framework for characterizing 26 intrinsic brain dynamics as a web of cross-attractor transitions and their energy 27 costs. The framework may be used to predict transitions and energy costs 28 associated with experimental or clinical interventions. 29 Significance Statement 30 The brain is a multifunctional system: different brain regions can coordinate flexibly 31 to perform different tasks. Understanding the regional and global structural con-32 straints on brain function is critical to understand cognition. Here, using a unified 33 biophysical network model, we show how structural constraints across scales jointly 34 shape the brain's intrinsic functional repertoire -the set of all possible patterns that 35 a brain could generate. The modeled functional repertoire is enriched by both the 36 flexibility of individual brain regions and their synergistic interaction in a global 37 network. By carefully examining the modeled repertoire, we found that human 38 resting-state functional connectivity was better predicted by transitions between 39 brain activity patterns rather than any specific pattern per se.

41
A fundamental goal of neuroscience is to understand how the structure of the brain 42 constrains its function [1]. The brain's multiscale and nonlinear nature makes it 43 challenging [2][3][4]. Relevant structural factors are scattered across scales, while 44 nonlinearity makes it difficult to predict dynamic properties of the whole from its 45 parts. Nonlinear dynamical models serve as essential tools for bridging structural 46 and functional understanding of the brain [5,6]. When equipped with plausible 47 biophysical and connectivity features of the brain, biophysical network models have 48 successfully provided insights into resting brain dynamics (e.g. [7][8][9]). However, their 49 dynamic capability has not been extensively studied or fully utilized with regard to 50 key nonlinear features such as multistability (c.f. [10]) and rhythmic activities. Here 51 we first provide a systematic analysis of how the dynamic landscape of the model 52 brain, its multistability and temporal diversity, vary with structural properties across 53 scales. Further, we show that functional connectivity patterns as observed in human 54 resting-state fMRI emerge spontaneously from the multistable landscape, reflecting 55 the coordination between brain regions across attractors in the dynamic landscape. 56 Intrinsic brain dynamics have long been observed [11,12], but often treated 57 as a baseline subtracted from task-positive activities. This baseline, however, is 58 more active than meets the eye: it consumes the largest fraction of the brain's 59 energy resources, while task-related consumption adds little [13]. It constrains task 60 performance and related neural activities across multiple time scales [14][15][16], and 61 sustains alteration in neurological and psychiatric disorders [17]. Neuroimaging 62 studies reveal the richness of such intrinsic brain dynamics and its subnetworks 63 [18][19][20][21][22][23][24]. These empirical findings drive the development of large-scale models of 64 resting brain dynamics (see [25] for a review), to reveal its functional purpose [5,26] 65 and as tools for mechanistic diagnostics of psychiatric disorders [27]. 66 Nonlinear dynamical systems are often chosen over linear ones to account for the 67 ubiquitous multistability and coordinated rhythmic activities in biological systems-68 the brain is one of the best examples [28][29][30][31][32]. To say that a system is multistable is 69 to say that multiple stable patterns of activities (attractors) are all achievable by 70 the system. Which pattern is retrieved depends on the external input or intrinsic 71 noise. Multistability of the brain signifies its multi-functionality and its ability 72 to form memory as persistent patterns of activity [2, 28, 33, 34]. The switching 73 dynamics between functional networks in the resting brain is thought to reflect noise-74 driven exploration of the underlying multistable landscape, i.e. the brain's intrinsic 75 functional repertoire [35]. On the other hand, brain dynamics is parsed in time 76 by a hierarchy of diverse rhythmic activities [32,36,37]. The coordination across 77 rhythmic activities at diverse frequencies gives rise to dynamic and flexible integration-78 segregation across multiple scales [32,38]. Disruption of such temporal coordination 79 is associated with neuropsychiatric disorders [39]. Together, multistability and 80 temporal diversity are key features of intrinsic neural dynamics that theorists seek 81 to capture using large-scale nonlinear dynamical models [10,40,41]. 82 In the present work, we use a nonlinear dynamical model that incorporates 83 structural properties of the brain across scales-from neurons, to local populations, 84 and to large-scale networks. The model is a formal unification of the Wilson-Cowan 85 model [42,43] and the reduced Wong-Wang model [7,8,44]. The Wilson-Cowan 86 model [42,43] is a population-level model of brain dynamics, widely used for modeling 87 multistability and rhythmic activity in large-scale brain networks [40,[45][46][47][48][49][50]. It 88 is unconstrained by the biophysical properties of underlying neurons that may be 89 important for predicting the outcome of electrical or pharmacological stimulation in 90 empirical settings (e.g. [51,52]). The reduced Wong-Wang model is a biophysical 91 network model constrained by biologically plausible parameters at the neuronal 92 level [34, 53-57]. Its noise-driven dynamics near certain attractors has been used 93 to capture resting-state functional connectivity (e.g. [8,9]). However, its reduced 94 nonlinearity limits the multistability and makes it less viable for studying oscillations. 95 Our unified model retains both the nonlinearity of the Wilson-Cowan model and 96 the biophysical constraints of the reduced Wong  With this model, we show how the intrinsic dynamic landscape of the brain 98 can be shaped by structural properties across scales, and how multistability of the 99 landscape gives rise to functional connectivity patterns and its energy demands. Our model describes the whole-brain dynamics as the mean-field activity of neuronal 103 populations in each brain region. Each model region contains a pair of excitatory 104 (E) and inhibitory (I) populations, whose activity is described by the local model 105 (Figure 1a Figure 1: Construction of a dynamic mean-field model of the human brain. (a) The model brain (global model) consists of a network of brain regions (local model). The local model (black box) describes the interaction between two local neural populations -one excitatory E and one inhibitory I -via two excitatory connections (red, w EE and w EI ) and two inhibitory connections (blue, w II and w IE ). The excitatory population of each brain region can further receive input (gray arrow, I E ) from other regions, allowed by the large-scale structural connectivity (red dashed connections). (b) Nodes in the global model corresponds to 66 anatomical regions of the human brain, which can be linked together by the human connectome (see text). Regions are index from 1 to 66 (1-33 on the right hemisphere, 34-66 on the left hemisphere is reverse order, following [7]). The local model is described by the equations, (2) The activity of each population has a natural decay time of τ E and τ I respectively. Each population's activity tends to increase with the fraction of closed channels (1 − S p ) and the population firing rate (H p ), scaled by a factor γ p for p ∈ {E, I}. This is described by the second term on the right-hand-side of equation 1-2. H E and H I are transfer functions that map synaptic current input to population firing rate of the excitatory and the inhibitory population respectively (for example, H E shown in Figure 2 as a black curve). In particular, they are sigmoidal functions of the form whose output increases with input monotonically and saturates at r max -the maximal 119 firing rate limited by the absolute refractory period of neurons (around 2 ms in 120 certain cell types [61,62] Figure 2) into a sigmoidal form (black 126 solid line in Figure 2), while retaining the original value of parameters a p , b p , and 127 d p (shown in Table 1). The parameters were chosen to approximate the average 128 response of a population of spiking pyramidal cells (p = E) and interneurons (p = I) 129 respectively, incorporating physiologically plausible parameters [44,57]. Interaction between local populations is modulated by four coupling parameters w pq 0 in equation 1-2, indicating the influence from the local population p to q, where p, q ∈ {E, I} (Figure 1 left box). These coupling parameters reflect the local structural connectivity. The local populations are also capable of responding to external current inputs denoted as I E and I I in equation 1-2, respectively. Importantly, such input can come from other brain regions in a globally connected network ( Figure 1 are the synaptic gating variable of the excitatory and the inhibitory population of the i th brain region respectively, and ξ (i) • is a noise term scaled to an amplitude σ. The state of all excitatory populations is denoted as a vector S E , the i th element of which is S (i) E . The global input to the i th brain region depends on both its connectivity with, and the ongoing state of, other brain regions, where N denotes the total number of brain areas, C ij 0 the long-range structural connectivity from the j th to the i th brain region and G is a global coupling parameter that controls the overall level of interaction across brain regions. Since C ij is only intended to represent long-range connectivity, we let C ij = 0 for any i = j to preclude recurrent connections. For the effects of G and C ij to be independently comparable, here we impose a normalization condition on the matrix norm, Since the global coupling parameter G modulates the level of input to each brain 131 region, one would expect it to have comparable influence on the local dynamics as 132 I E in the local model (equation 1).

133
Next, we discuss its formal connection to two well-studied mean-field models of 134 brain dynamics, namely, the Wilson-Cowan model (Section 2.2) [42,43]  Formally, the above model can be considered a special variant of the Wilson-Cowan model [42,43]. Though the specific interpretation of certain parameters differ, the two models describe similar dynamic mechanisms of population-level interaction. The Wilson-Cowan model, in its initial form [42], concerns the dynamics of a pair of interacting excitatory and inhibitory neuronal populations. The activities of the two populations are denoted as E(t) and I(t)-the proportion of firing excitatory/inhibitory cells averaged over a period of time (the refectory period). The model takes the form τ E and τ I are time constants of the dynamics of the excitatory and inhibitory 138 population respectively. c • 's are the coupling parameters between the two population. 139 Coefficients k • and r • result from a temporal coarse-graining procedure in the 140 initial derivation (see [42] for detail). S • is a sigmoidal transfer function, rising 141 monotonically from 0 to 1 with non-negative input. P and Q are external inputs 142 to their respective populations. If we divide both sides of equation 8-9 by the time 143 constants, we are looking at the same general form as equation 1-2.

144
The main difference is between the respective transfer functions. Wilson and Cowan [42] chose a particular form of S for mathematical analysis: where parameter a determines the maximal slope of the function S and parameter θ 145 the location of the maximal slope. Technically, Wilson and Cowan [42] only requires 146 S to be of a general sigmoidal form. It may reflect the average response of a 147 population of neurons with heterogeneous firing thresholds or heterogeneous afferent 148 connections. The distribution of the said thresholds or connections is reflected in 149 the parameters a and θ.

150
In other words, the choice of the transfer function and the parameters is non-151 specific to a predefined microscopic model. Moreover, Wilson and Cowan [42] took 152 a function-oriented approach to analyzing the model. The key was whether the 153 model was able to produce fundamental behaviors expected from a neural model-154 multistability, hysteresis, and oscillation-for some specific choice of parameters 155 and transfer function. Qualitative conclusions from their analysis depend on the 156 general geometric properties of the transfer function rather than the specific form of 157 equation 10.

158
The transfer function of the present model (equation 3) follows the general 159 geometric properties assumed by Wilson and Cowan [42]. The difference is that the 160 parameters in equation 3 are associated specifically with a microscopic model [57], a 161 network of leaky integrate-and-fire neurons with biologically plausible parameters, 162 as inherited from the reduced Wong-Wang model [8,44]. This choice provides a 163 channel of correspondence between parameters of the models at different scales of 164 description. To expand on this point, we next elaborate on the connection between 165 the present model and the reduced Wong-Wang model. The present model is also a variant of the Wong-Wang model [44] and its highdimensional generalizations, here referred to as the reduced Wong-Wang model [7-9]. In particular, we consider the model of whole-brain dynamics [8,9], following the same notations as in equation 4-6, wherẽ with p ∈ {E, I} denoting the excitatory and the inhibitory population respectively 168 (see Figure 2 dashed line forH E ). The parameters a p , b p and d p were chosen such that 169 H p approximates the average firing rate of an ensemble of leaky integrate-and-fire 170 neurons receiving uncorrelated noisy inputs.

171
More specifically, the sub-threshold dynamics of the membrane potential V (t) of each neuron can be described as where C m is the membrane capacitance, g L the leak conductance, and V L the resting 172 potential of the membrane. The total synaptic input current I syn (t) is a random 173 process with an average µ C and standard deviation σ C . When V (t) reaches a 174 threshold V th , the neuron emits a spike after which the membrane potential returns 175 to a reset voltage V reset and stays there for a duration τ ref , i.e. the refractory 176 period.

177
The average firing rate ν of an ensemble of such neurons can be derived from the Fokker-Planck approximation that describes the evolution of the membrane voltage distribution of an ensemble of neurons (see e.g. [ where τ m = C m /g L is the membrane time constant, σ V = √ τ m σ C /C m the standard deviation of the depolarization, erf(x) the error function and V ss the steady state voltage The transfer function employed by Wong and Wang [44,63], i.e. equation 13 178 with appropriate choice of parameters, is a good approximation of equation 15 when 179 the input level is low.  (Table 2). In other words, it allows one to use empirically 183 measurable quantities at the neuronal level to directly constrain the the transfer 184 function and the entire model. This is a major difference with the Wilson-Cowan 185 model [42,43].  Table 2: Biophysical parameters of a single leaky-integrate-and-fire neuron. If two parameter values are provided in the right column, the first value is for a generic pyramidal cell and the second is for a generic interneuron. Differences between the biophysical parameters of different cell types lead to differences in the transfer functions (equation 13).
According to the first passage equation 15, the firing rate ν is a sigmoidal function of the input, which saturates at r max ≡ 1/τ ref . This is not the case, however, for the transfer functionH (equation 13).To make the transfer function a better approximation of the first passage equation and at the same time retain the mapping between their parameters, we can simply convert the transfer function by substituting the numerator ofH as below Thus, we obtain the transfer function used in the present model (equation 3). As 187 shown in Figure 2, H p matchesH p for low levels of input but flattens out eventually 188 at r max as one would expect from equation 15.

189
In short, the present model is endowed with the geometric properties of the 190 Wilson-Cowan model [42,43] (Figure 1b). The original diffusion imaging (dMRI) data were 199 obtained using a customized Siemens 3T scanner at Washington University in St. 200 Louis, with a standard 32-channel head coil, with TR = 5520 (ms), TE = 89.5 201 (ms), 1.25 (mm) isotropic voxels, b=1000, 2000, 3000 (s/mm 2 ). T1 images were 202 obtained using 3D magnetization-prepared rapid gradient echo sequence (MPRAGE) 203 with TR = 2400 (ms), TE = 2.14 (ms), and 0.7 (mm) isotropic voxels. The HCP 204 minimally processed data were further processed using MRtrix3, including bias-205 field correction, multi-shell multi-tissue constrained spherical deconvolution with 206 a maximum spherical harmonic degree 8. 10 million probabilistic streamlines were 207 generated for each subject using the 2 nd -order Intergration over Fibre Orientation 208 Distributions algorithm (iFOD2) [66] and anatomically-constrained tractography 209 (ACT) [67] (FOD amplitude threshold = 0.06, step size = 0.625 mm). Each streamline 210 was assigned a weight using spherical-deconvolution informed filtering of tractograms 211 (SIFT2) [68]. Connection strengths between ROIs are summed weights of the 212 associated streamlines. Intra-ROI connections are removed. Subjects' connectivity 213 matrices are normalized according to equation 7 before and after averaging. Human functional connectivity used in the present study is estimated using the 216 resting-state fMRI (rfMRI) data of the same 11 unrelated subjects from the Human 217 Connectome Project [59] as the structural connectivity above. rfMRI scans were 218 acquired using EPI sequences with TR = 720 (ms), TE = 33.1 (ms), flip angle = 219 52 • , voxel size = 2.0 (mm, isotropic), multiband factor = 8. Four runs of rfMRI 220 scan were obtained from each subject in 2 separate days (2 runs in each day with 221 opposite phase-encoding direction: RL and LR). Each run last 14 min 33 s (1200 222 TR). . First, a reference volume and its skull-stripped version were 227 generated using a custom methodology of fMRIPrep. A deformation field to cor-228 rect for susceptibility distortions was estimated based on two echo-planar imaging 229 (EPI) references with opposing phase-encoding directions, using 3dQwarp [73] (AFNI 230 20160207). Based on the estimated susceptibility distortion, an unwarped BOLD 231 reference was calculated for a more accurate co-registration with the anatomical 232 reference. The BOLD reference was then co-registered to the T1w reference us-233 ing bbregister (FreeSurfer) which implements boundary-based registration [74]. 234 Co-registration was configured with nine degrees of freedom to account for dis-235 tortions remaining in the BOLD reference. Head-motion parameters with respect 236 to the BOLD reference (transformation matrices, and six corresponding rotation 237 and translation parameters) are estimated before any spatiotemporal filtering using 238 mcflirt [FSL 5.0.9, 75]. The BOLD time-series, were resampled to the fsaverage 239 surface space. Several confounding time-series were calculated including framewise 240 displacement (FD), DVARS and three region-wise global signals. FD and DVARS 241 were calculated for each functional run, both using their implementations in Nipype 242 [following the definitions by 76]. The three global signals were extracted within the 243 CSF, the WM, and the whole-brain masks.  Two types of inter-regional coordination in the model are computed and compared 258 to the human functional connectivity. One estimates how every two model regions 259 move together across different attractors throughout the entire dynamic landscape, 260 namely cross-attractor coordination. The other estimates how two model regions 261 move together during noise-driven exploration near a single attractor, namely within-262 attractor coordination.

263
For cross-attractor coordination, the attractors are first discretized by replacing 264 the state of each region with an integer, determined by a 30-bin histogram of regional 265 states S (i) E (Section S3). This procedure allows us to quantify whether two regions 266 move up and down together across different attractors without considering the 267 distance between attractors and the continuous shift of the position of each attractor. 268 The cross-attractor coordination between region x and y is computed as the Spearman 269 correlation (ρ) of their discrete states across a fixed set of attractors, determined by 270 structural parameters G, w EE , and w EI (see Figure 6 for a conceptual illustration). 271 Note that this procedure only requires the coordinates of each attractor through, 272 e.g. the computation of a bifurcation diagram (Section S1)-simulation of the time 273 series is not required. The distance between attractors is computed separately as 274 the difference between the average state of the excitatory populationS E , which can 275 be interpreted as the energy cost associated with keeping additional x% synaptic 276 channels open. Within-attractor coordination is the Spearman correlation between 277 simulated time series of the excitatory population of region x and y, given an initial 278 condition right at an attractor, a moderate level of noise σ = 0.01 and a duration 279 T = 864 s (14 min 33 s to match the human data).

281
In the following sections, we first examine the dynamic repertoire of isolated brain 282 regions (the local model; Section 3.1) and how nonlinearity in the present model 283 enhances multistability and produces realistic oscillations (Section 3.2). We fur-284 ther show how long-range connectivity between these brain regions interacts with 285 local properties in shaping global multistable landscapes (Section 3.3), human func-286 tional connectivity patterns (Section 3.3), and temporal diversity across regions 287 (Section 3.4). The numerical results are illustrated in the main text while the 288 corresponding analytical supports are provided in the Supplementary Materials. plementary Materials) shows that nonlinearity in the dynamics can essentially be 295 controlled by two local structural properties: the strength of excitatory-to-excitatory 296 connection w EE and the strength of the excitatory-to-inhibitory connection w EI . 297 Geometrically, the two structural properties "twist" the nullclines (dashed lines in 298 Figure 4). Specifically, stronger w EE introduces a deeper twist and fold of the 299 red nullcline (compare Figure 4a and d), whereas stronger w EI introduces a more 300 vertical twist of the blue nullcline (compare Figure 4d and e). These twists are the 301 key sources of dynamic complexity-multistability and oscillation. For example, 302 when w EE is sufficiently large (equation S19), multistability becomes possible: the 303 folded red nullcline allows for multiple intersections with the blue nullcline, and 304 potentially a greater number of attractors (compare Figure 4c to a; see analytical 305 results in Section S7: Multistability). When w EI is sufficiently large (equation S27), 306 oscillatory activity becomes possible (analytical results in Section S7: Oscillation). 307 Moreover, the combination of large w EE and w EI gives rise to sustained oscillation 308 (equation S50). The characteristic frequency of such oscillation further depends on 309 the specific values of w EE and w EI . Note that the general qualitative effects of 310 these two local structural properties are consistent with those of the Wilson-Cowan 311 model, but the specific boundaries at which transitions occur are determined by the 312 biophysical constraints inherited from the reduced Wong-Wang model (see equa-313 tions S27, S50). Analytical results (Section S7) provide detailed quantification of 314 how these boundaries are shifted by different local structural properties.

315
To maintain a sufficient twist in the red nullcline (red dashed line in Figure 4) 316 and associated multistability and oscillation, inhibitory-to-excitatory feedback w IE 317 needs to be proportional to self-excitation w EE (c.f. equation S9). In the present 318 study, we simply let w IE = w EE . This equality is a simpler alternative to the 319 Feedback Inhibition Control adopted in [8] in both numerical and mathematical 320 analyses.
The vector fields (arrows) reflect the underlying dynamics at different points in the state space. Gray trajectories following the vector fields are solutions of the local model (equation 1, 2) given a fixed sets of ten different initial conditions. Nullclines (dashed lines) indicate where the flow of the dynamics is either purely vertical (red) or purely horizontal (blue). The intersections between the nullclines are the fixed points. Different types of fixed points are labeled with different markers (see legend). A fixed point is stable (×) if nearby trajectories converge to it over time, unstable (+) if nearby trajectories diverge from it, or a saddle ( * ) if nearby trajectories approach it in some direction(s) but diverge from it in some other direction(s). A fixed point is said to be a spiral (•) if trajectories near the fixed point rotate either towards the fixed point (damped oscillation) or away from the fixed point (sustained oscillation or limit cycle in the present case). Strong oscillation mainly appears on the ascending branch of the red nullcline. Overall, we see that local connectivity defines the dynamics in each regime essentially by controlling the geometry of the nullclines.  324 44] to more complex scenarios. As expected, the dynamics of the present model 325 match that of the reduced Wong-Wang model for low levels of excitation, i.e. weak 326 local excitatory connectivity (Figure 5a-b). In a regime of stronger local excitatory 327 connectivity, as explored in [9], the two models diverge (Figure 5c-d). In the present 328 model (Figure 5c), all trajectories are well-confined within a physiologically plausible 329 range-state variables S E and S I denote the fraction of open channels, which by 330 definition are between 0 and 1. In contrast, certain trajectories of the reduced 331 Wong-Wang model (Figure 5d) overshoot beyond the physiologically plausible range. 332 The effect of added nonlinearity in the present model manifests through the curvature 333 of the blue nullclines, which confines the flow of oscillatory activities and creates 334 extended multistability (see e.g. Figure 4b). Thus, the present model is more 335 suitable for studying key nonlinear dynamical features in the resting brain.  Table 1). The resulted dynamics are virtually identical. (c) and (d) show a similar comparison between the two models in an oscillatory regime, where the local excitatory connectivity is stronger (w EE = 4, w EI = 1). While the dynamics of the present model (c) is well confined within a realistic range (S E , S I ∈ [0, 1]), it is not the case for the reduced Wong-Wang model (d).
3.3 Multistable landscape of the brain shaped by structural 337 properties across scales 338 In Figure 6, we provide basic intuitions about dynamic landscapes and their repertoire 339 of attractors as a system-level description of intrinsic brain dynamics. An attractor 340 in the global model represents a stable pattern of activation over the whole brain 341 (black boxes in Figure 6). The dynamic landscape of the brain (Figure 6a-c) defines 342 the repertoire of attractors, thereby possible patterns of activity (e.g. a1-a4), paths 343 of transitions (black arrows), and the coordination between brain regions during 344 transitions (2-by-2 matrices). As the landscape itself changes, some attractors may 345 be destroyed (Figure 6, a→b, or a→c) or created (b→a, c→a)-a discrete change of 346 the repertoire called bifurcation. In the present model, local and global structural 347 connectivity control the shape of the landscape, and thereby, the repertoire of 348 attractors, possible transitions, and associated inter-regional coordination.
349 Figure 6: Conceptual illustration of the whole-brain dynamic landscape, bifurcation, and phase transition. A multistable dynamic landscape (a) contains multiple attractors, shown as troughs occupied by purple balls. For the present model, each attractor is interpreted as a distinct pattern of activation over the whole brain (a1-a4). Influenced by external input or intrinsic noise, the model brain may transition from its current state (a1, bright purple ball) to a different one (a2, a3, or a4, dim purple balls), indicated by black arrows. Structural properties of the model brain can alter the shape of the landscape, causing some attractors to appear or disappear through a process mathematically named bifurcation (a→b, a→c, or the reverse). By modifying the repertoire of attractors, bifurcation alters the set of possible transitions and the coordination between regions during transitions: in landscape (a), the left and right hemisphere can be co-activated during a transition (a1→a3), or activated independently through other transitions (a1→a2, or a1→a4); in contrast in landscape (b), the left and right hemisphere can only be co-activated, and in (c), only activated independently. Numerically, a repertoire of attractors can be represented as a matrix, where each row is an attractor and each column is a brain region (e.g. the 4-by-2 matrix below landscape a, for 4 attractors -a1-a4, and 2 regions -left and right hemisphere). Overall cross-attractor coordination between regions can be estimated by rank correlation between the columns of the repertoire matrix. The resulted square, coordination matrix summarizes how brain regions transition together over the entire landscape, serving as a signature of the landscape (e.g. 2-by-2 matrices next to a-c). In more complex landscapes (not shown), there are many more attractors, and they correspond to subtler patterns of activation (d,e; see also Figure 7). The coordination between brain regions during a transition is correspondingly more complex (f=d-e), with some regions co-activated (red) while others co-deactivated (blue). Figure 7 shows such dynamic landscapes and their changes more abstractly as 350 bifurcation diagrams (see Section S1 for computational details). In each bifurcation 351 diagram, the y-coordinate of each colored point indicates the position of an attractor 352 by an order parameter (e.g. average activity of all excitatory populationsS E ), 353 and x-coordinate a control parameter that modulates the shape of the underlying 354 dynamic landscape (e.g. global coupling G in equation 6). Each vertical slice 355 of a bifurcation diagram contains the repertoire of attractors and repellers in a 356 fixed landscape (Figure 7h), corresponding to stable and unstable patterns of brain 357 activity. Horizontal stripes are formed when an attractor changes continuously with 358 the landscape, but destroyed when they merge with a black stripe (unstable patterns) 359 at a bifurcation. Thus, the number of colored stripes reflect the complexity of the 360 landscapes: more stripes, more attractors.

361
In the simple case where all brain regions are connected to each other at the same 362 strength (uniform connectivity; Figure 7d-f), stronger local excitatory connection 363 (e,f) produces a more complex landscape (3 attractor stripes) than a weak one (d; 2 364 attractor stripes). These bifurcation diagrams are very similar to those of a single 365 brain region (Figure 7a-c), in terms of the number of attractors and the presence 366 oscillation. In fact, the whole brain ( Figure 7e) moves up and down together between 367 discrete state of activation very much like a single region (Figure 7b). One remarkable 368 difference between the uniformly connected global models and the corresponding 369 local models is that, with the absence of persistent input, the global model can retain 370 memories of prior input while the local model cannot. That is, when input I E = 0, 371 the local model is monostable (bottom red stripes in Figure 7a-c), i.e. returning 372 to the same state regardless of prior input. The global model (equation 4-6) by 373 definition does not receive external input; yet the model is multistable for a sufficient 374 amount of global coupling, e.g. G > 1 (Figure 7d-f). We further substantiate and 375 generalize this result analytically and numerically (Section S8 Multistability) to 376 the case where each isolated brain region is monostable for any I E . These findings 377 suggest that the coupling between brain regions can synergistically create a functional 378 repertoire or memory capacity that isolated brain regions do not possess.

379
Given a realistic global structural connectivity (human connectome; Figure 7g-i), 380 the complexity of the whole-brain dynamic landscape is dramatically increased: 171 381 attractor stripes in (g), 610 in (h), and 682 in (i) (per single-linkage clustering). 382 Correspondingly, the patterns of activation (Figure 7f) are also more complex, with 383 greater differentiation between regions; the coordination between brain regions across 384 attractors is consequently more flexible and subtle (Figure 6f). The heterogeneous 385 nature of the human connectome breaks the large-scale spatial symmetry of the 386 model brain, creating more functional differentiation between brain regions and a 387 greater functional complexity for the whole brain. In short, the complexity of the 388 global dynamical landscape is a joint product of strong local excitatory connection 389 and complex topology of the large-scale network. See Section S8 for additional 390 analytical supports. (c) w EE = 2.8 and w EI = 1. Overall, local connectivity increases from (a) to (c). The activity of the excitatory population S E is used as an order parameter, indicating the location of each attractor. The external input I E is used as a control parameter. Each point in the diagram indicates the location of a particular fixed point. The color denotes the type of each fixed point: non-black points represent attractors, black points unstable fixed points that are not associated with a limit cycle. Horizontal stripes indicate that the attractors are changing continuously with the control parameter for a certain range. All (a)-(c) have a upper stripe and a lower stripe. (b)-(c) have an additional stripe in the middle, where the brain region oscillates. Insets of (b) and (c) show the oscillation frequency of the brain region as a function of the input current. Each stripe corresponds to a discrete level of activation for a single brain region (circled brains in b; color indicates discrete S E levels, shown in circled legend). (d)-(f) show the corresponding bifurcation diagrams for three uniform global networks, i.e. the large-scale structural connectivity C ij 's are identical between any two brain regions i = j (equation 6). The average activity of all excitatory populations (S E ) is used as an order parameter and the global coupling G (equation 6) as an control parameter. Each attractor stripe corresponds to a pattern of activation over the whole brain (circled brains in (e) show S (i) E 's on the left hemisphere). Similarly, (g)-(i) show the corresponding bifurcation diagrams for three realistic global networks, i.e. C ij 's reflect the human connectome (see text for details). Here each vertical slice (gray line in h) contains the attractor repertoire of a fixed dynamic landscape shaped by the human connectome. Each attractor repertoire is associated with a matrix describing the coordination between brain regions across attractors (e.g. Figure 8b). See Figure 6 for a cartoon illustration of an attractor repertoire and its associated cross-attractor coordination matrix.
What do these complex landscapes say about the dynamics in the human brain? 392 Below we show that large-scale functional connectivity patterns in human fMRI data 393 (Figure 8a) can be aptly explained by how brain regions transition together across 394 attractors-cross-attractor coordination (Figure 8b), better than how they coordinate 395 within any single basin of attraction-within-attractor coordination (Figure 8c).

396
Human functional connectivity (Figure 8a) is calculated from the resting fMRI 397 data from the Human Connectome Project [59], averaged across the same subjects 398 whose structural connectome is incorporated in the model (Figure 8e; see Section 2.4 399 and Section S4 for more details). The human functional connectivity exhibits large-400 scale symmetry across two hemispheres. That is, functional connectivity within the 401 right hemisphere is similar to the connectivity within the left hemisphere (compare 402 upper-left and lower-right block in Figure 8a), and similar to the connectivity across 403 the left-right hemispheres (compare upper-left and lower-left block in Figure 8a). 404 This dominant feature of human resting brain dynamics is well preserved in the 405 model by the cross-attractor coordination between model regions (Figure 8b; see 406 Section 2.4.3 for detail, Figure 6a-c for intuition regarding cross-attractor coordina-407 tion). Quantitatively, this is also reflected as the consistent model fit for both intra-408 and inter-hemispheric connectivity (dashed lines in Figure 8d; black dashed line: 409 whole-brain fit, ρ = 0.591, p < 0.001, R 2 = 0.349; red dashed line: intra-hemispheric 410 fit, ρ = 0.582, p < 0.001, R 2 = 0.338; blue dashed line: inter-hemispheric fit, 411 ρ = 0.604, p < 0.001, R 2 = 0.364), and across model parameters (Figure 9a-c).

412
The above finding paints a rather different picture than previous theories in 413 which the resting brain is thought to explore the dynamic landscape near a single 414 attractor, delicately poised near a bifurcation (e.g. [7,9]). For comparison, we 415 simulated functional connectivity within each of the individual attractors involved 416 in cross-attractor coordination. The best-fit within-attractor coordination matrix 417 is shown in Figure 8c (ρ = 0.463, p < 0.001, R 2 = 0.190), clearly missing the 418 symmetry between intra-and inter-hemispheric connectivity observed in humans 419 (Figure 8a) and cross-attractor coordination of the model (Figure 8b). Moreover, 420 within-attractor coordination does not capture inter-hemispheric connectivity in 421 humans as well as intra-hemispheric connectivity (solid lines in Figure 8d), and overall 422 fits worse than cross-attractor coordination (dashed lines in Figure 8d; for comparison, 423 best-fit intra-hemispheric within-attractor coordination gives ρ = 0.534, p < 0.001, 424 R 2 = 0.285). The distinction between within-and cross-attractor coordination is 425 more drastic when the contribution of structural connectivity is controlled using 426 partial correlation (Figure 8f; contrast best-fit within-attractor intra-hemispheric 427 coordination: ρ = 0.218, p < 0.001, ∆R 2 = 0.0269, F (1, 1053) = 52.3, v.s. cross-428 attractor intra-hemispheric coordination: ρ = 0.417, p < 0.001, ∆R 2 = 0.0991, 429 F (1, 1053) = 221). These findings suggest that cross-attractor coordination better 430 captures features of human resting brain dynamics that are unexplained, linearly, 431 by the structure. 432 Figure 8: Human resting functional connectivity better captured by cross-attractor, rather than withinattractor, coordination between model brain regions. Human functional connectivity matrix (a) is calculated using the resting fMRI data from the Human Connectome Project [59], averaged over the same subjects whose average structural connectome defines the long-range structural connectivity in the model (e). Regions (columns and rows) are ordered symmetrically for the left and right hemisphere (see Figure 1b) to reveal the large-scale symmetry of resting brain dynamics. White lines delineate the matrix (a) into four blocks, describing the functional connectivity within the right hemisphere (upper left block), within the left hemisphere (lower right), and between two hemispheres (lower left/upper right). Functional connectivity within the hemispheres are similar to each other (symmetric along anti-diagonal), and also similar to inter-hemispheric connectivity (symmetric along white lines). This large-scale functional symmetry is well captured by inter-regional coordination in the model brain (w EE = 2, w EI = 1, G = 2.22) across attractors (b; 97 attractors shown in Figure 7h for G = 2.22). Such symmetry is not captured by the coordination within any of the said attractors (c; best fit within-attractor coordination matrix). In (d), solid line shows the distribution of correlation coefficients between human functional connectivity (a) and within-attractor coordination matrices (not shown except c). All of which are lower than the correlation (dashed lines in d) between human functional connectivity (a) and the cross-attractor coordination matrix (b). (f) shows a similar comparison as (d) using partial correlation controlling the contribution of structural connectivity (e). (*** p<0.001, Bonferroni corrected) Figure 9: Model-human similarity consistent for within-and cross-hemispheric coordination, affected by energy constraints. (a-c) Cross-attractor coordination in the model is highly similar to human functional connectivity across different local excitatory connectivity: (a) w EE = 0.7 and w EI = 0.35, (b) w EE = 2 and w EI = 1, (c) w EE = 2.8 and w EI = 1, as in Figure 7. Red solid lines indicate the model-human similarity (Spearman correlation) for intra-hemispheric connectivity, blue ones for inter-hemispheric connectivity, black ones for whole-brain connectivity. Dashed lines and shaded areas indicate the means and 95% confidence intervals of the null distributions, constructed using 1000 random permutations of ROI labels (red, blue, black distributions are largely overlapping not visually distinguishable). Black solid lines in (d-f) are reproduced from (a-c) indicating the model-human similarity over the whole brain. Dashed lines in (d-f) indicate the model-human similarity when the model is not allowed to cross the largest energy gap between attractors. The shaded area (∆ρ) indicates the loss of similarity between the model and the humans.
As alluded to, cross-attractor coordination in the model captures human func-433 tional connectivity for a wide range of structural configurations (Figure 9a-c; see also 434 Figure S3). Then, what roles do local (w EE , w EI ) and large-scale structural proper-435 ties (e.g. global coupling G) play in producing human-like functional connectivity? 436 The answer lies in the energy gap between different attractors. In the computation 437 of cross-attractor coordination matrices (e.g. Figure 8b, Figure 6a-c), we considered 438 only whether two brain regions move up and down together across the dynamic 439 landscape, not how difficult the movements are. In fact, the average of the pattern 440 change (∆S E , e.g. Figure 6f) between attractors (e.g. Figure 6, d→e) reflects an 441 energy gap -the energy needed to keep x% additional synaptic channels open. 442 The similarity to human functional connectivity drops greatly if the model brain is 443 not allowed to traverse the largest energy gap in the dynamic landscape (dashed 444 lines in Figure 9d-f). This energy constraint has a greater impact on model-human 445 similarity when the local structural connectivity is weaker (area of the shaded region 446 grows from Figure 9d to f; bars in Figure 10) and the global structural connectivity 447 is stronger (height of shaded regions grows with G in Figure 9d-f). The loss of 448 similarity grows with the maximal gap size, especially for a gap size greater than 449 0.2 ( Figure 10). Thus, local and global structural connectivity both influence the 450 energy costs associated with cross-attractor coordination, and thereby how likely 451 it is for human-like functional connectivity patterns to emerge under energy con-452 straints. Overall, stronger local connection reduces the energy gaps (Figure 11a-c) 453 and stronger global connection (G) increases the energy gaps (Figure 11a)-local 454 and global structural connectivity pull the energy cost in different directions. 455 Figure 10: Loss of model-human similarity (∆ρ/ρ all ) due to energy constraints depends on the maximal energy gap and local excitatory connectivity. Each point in the scatter plot represents the loss of model-human similarity due to the inability to cross the maximal energy gap given a specific combination of global coupling G, and local connectivity w EE and w EI . Overall, the loss increases with the gap size, which in turn depends on G (Figure 11a). The average loss (bars) decreases with increasing local excitatory connectivity (w EE , w EI ). (*** p<0.001 with Tukey HSD; error bars are standard errors throughout the text) Figure 11: Local (w EE , w EI ) and global (G) structural connectivity jointly shape the energy cost of cross-attractor coordination. (a) Overall, the maximal (solid lines) and average energy gaps (dashed lines) increase with global coupling G, though there is a transient decrease when maximal energy gap is less than 0.2. Both types of gaps decrease with increasing local connectivity w EE , w EI (b,c). (*** p<0.001 with Tukey HSD) Is the energy demand for cross-attractor coordination uniform over the whole 456 brain? Here we compare the energy demand within the Default Model Network 457 (DMN)-a signature of the resting brain-and other brain regions (Figure 12). 458 Cross-attractor coordination is more energy demanding within DMN than the rest of 459 the brain, measured by the ratio between the maximal energy gap within DMN v.s. 460 others (gap ratio>1, above dashed line in Figure 12a; or equivalently log gap ratio>0 461 in Figure 12b). DMN dominates the energy consumption primarily when the overall 462 energy consumption of the whole brain is low (max energy gap≤ 0.2, Figure 12b), 463 and when local excitatory connections are stronger (left bars in Figure 12b). (a) DMN generally contains a greater maximal energy gap than other networks (gap ratio>1), especially when the maximal energy gap of the whole brain (x-axis) is small (e.g. ≤ 0.2). (b) On average, the log gap ratio between DMN and other networks is significantly greater than zero (i.e. gap ratio>1) when the whole-brain maximal energy gap ≤ 0.2 (left three bars), but very close to zero when the whole-brain maximal energy gap > 0.2 (right three bars). The dominance of DMN's energy demand increases with local excitatory connectivity (left three bars: blue > green > red). (* p<0.05, *** p<0.001 with Tukey HSD)

Local and global causes of temporal diversity 465
Now we turn to the structural constraints on temporal diversity. In particular, we 466 show how spectral properties of the simulated neural activities and correspond-467 ing hemodynamic responses are affected by both the diversity of local structural 468 properties and the structure of the large-scale connectome.

469
Given a uniform global network, temporal diversity across the whole brain can be 470 induced by the diversity of local excitatory-to-excitatory connection (w EE ), as shown 471 in Figure 13a. Brain regions with relatively weak w EE (blue) have low characteristic 472 frequencies around 10 Hz (alpha range), while brain regions with strong w EE (red) 473 have higher characteristic frequencies around 30 Hz (beta/gamma range). In other 474 words, the characteristic frequency of the oscillation increases monotonically with 475 w EE (see also Figure S4a). This is expected from the behavior of isolated brain 476 regions (Figure 3d). In addition to the expected diversity, signs of coordination 477 between regions can be seen as the wide-spread alpha peaks (Figure 13a). In 478 contrast, regions with a higher characteristic frequency (beta/gamma range) are 479 not as influential to other regions. That is, low-frequency oscillations, rather than 480 high-frequency ones, are responsible for global coordination.

481
The above observations concern high-frequency dynamics typically measured 482 using, e.g. electroencephalography (EEG) and magnetoencephalography (MEG). For 483 low-frequency dynamics typical for functional magnetic resonance imaging (fMRI), 484 we examine the low-frequency content (0.01-0.1 Hz) of the normalized power spectra 485 of BOLD activities, derived from the same simulated neural dynamics (see Section S2 486 in Supplementary Materials for details). The result is shown in Figure 13b: there is 487 no significant dependency of low-frequency power on w EE (Spearman correlation 488 ρ = −0.029, p = 0.81). In short, we find differential effects of local structural 489 diversity on neural dynamics at the time scales typical for different neural imaging 490 modalities.  E for i = 1, · · · , N . The spectrum for each brain region is color coded by the rank of w EE -blue to red indicate the smallest to the largest w EE . The peak frequency of these spectra clearly increases with w EE . (b) shows the rank of the low-frequency power of the corresponding BOLD signal, integrated over the frequency range [0.01, 0.1] Hz (see Section S2 for details), which depends little on the rank of w EE . (c) and (d) show results of similar analyses but for the second simulated trial, where the individual brain regions are identical (w (i) EE = 2 for all i) but the global structural connectivity is realistic, i.e. C ij here reflects the human connectome [59, 60] with the global coupling G = 2.5. Both low-frequency (d) and high-frequency (c) activities are highly affected by the degree of the brain region in the global network (rank color-coded).
On the other hand, temporal differentiation does not mandate the brain regions 492 themselves to be structurally different. As shown in Figure 13c-d, locally identical 493 brain regions can behave very differently due to the topology of the large-scale 494 network (human connectome as in Section 3.3). The influence of large-scale structural 495 connectivity on temporal diversity is manifested in both the high-frequency neural 496 dynamics (Figure 13c; Figure S5a) and the low-frequency power of the BOLD signals 497 (Figure 13d; Figure S5b). Specifically, the low-frequency power is inversely related 498 to the degree of each node (brain region) in the large-scale network (Spearman 499 correlation ρ = −0.584, p < 10 −6 ).

500
In Section S6, we demonstrate that the above effects are robust over 200 simulated 501 trials of the same parameter settings. Overall, both local (Figure 13a,b; Figure S4) 502 and large-scale structural connectivity (Figure 13c,d; Figure S5) contribute to the 503 diversification of local dynamics. The contribution of local structural differences is 504 stonger in a higher-frequency range (Figure 13d; Figure S4a), while the contribution 505 of global structural connectivity is stronger in a very-low frequency range (Figure 13d; 506 Figure S5b). Modeling real neural dynamics requires considering both sides of the 507 spectrum. The present work shows systematically how multistability and temporal diversity 510 of the brain can be shaped by structural constraints across scales using a biophysi-511 cally constrained nonlinear dynamical model. We show that human-like functional 512 connectivity emerges spontaneously from the multistable landscape as brain regions 513 transition across attractors in a coordinated manner. The work suggests a transi-514 tion-centered view of human functional connectivity, over a attractor -centered one. 515 The theoretical and empirical implications are discussed below.

516
The rich dynamics of a single isolated brain region can be effectively controlled 517 by two key local structural properties: local excitatory-to-excitatory connectivity 518 (self-excitation) and local excitatory-to-inhibitory connectivity. In the real brain, 519 local excitatory-to-excitatory connections are particularly abundant [78], and in our 520 model, they contribute indispensably to multistability (Section S7). Multistability 521 is a key source of biological complexity from molecular to social levels  At the large-scale network level, multistability can be created or amplified by the 531 synergistic interaction between mono-or multi-stable brain regions (Section S7 Mul-532 tistability). Different large-scale network structures have dramatically different 533 capability at amplifying local complexity: a realistic global network (Figure 7 g-i) is 534 much more powerful than a uniform one (Figure 7 d-f). The human connectome 535 breaks the spatial symmetry of the global model, whereas symmetry breaking is 536 often a key to complex dynamics [2, 38, 85-88]. The human connectome is also 537 endowed with specific features such as modularity, small-worldness, and multiscale 538 characteristics [1,[89][90][91]. A systematic study of how these features alter the geome-539 try of the global dynamic landscape is worthy of further theoretical investigation 540 (see Section S8).

541
Within the multistable landscape sculpted by the human connectome, coordina-542 tion between model brain regions across different attractors gives rise to human-like 543 functional connectivity (Figure 8b). Importantly, such cross-attractor coordination 544 better captures human functional connectivity than within-attractor coordination-545 synchronization between brain regions near a single attractor (Figure 8c,d). This 546 raises the possibility that functional connectivity patterns reflect transitions between 547 stable brain states (Figure 6f) more than the brain states themselves (Figure 6d,e). 548 This transition-centered view offers alternative perspectives on several theoretical 549 and empirical issues. First, with a within-attractor approach [7], similar models 550 exhibited much lower inter-hemispheric connectivity than that of humans. It was 551 attributed to an underestimation of structural connections across hemispheres [7, 552 9]. However, strong functional connectivity in humans is known to exist between 553 regions that are not directly connected [92]. Thus, strong inter-hemispheric func-554 tional connectivity may reflect the nonlinearity in the dynamical system despite 555 the weak structural connectivity. Cross-attractor approach takes into account such 556 nonlinear effect. In contrast, within-attractor coordination can be approximated 557 by a linear dynamical system, which closely depends on the structural connectivity 558 [7, 8] (c.f. Figure 8d v.s. f). Second, cross-attractor coordination in the present 559 study is measured over the entire dynamic landscape, which itself is time invariant. 560 Empirically observed stability and convergence of human functional connectivity [93, 561 94] may reflect this invariance of the underlying landscape. The static landscape also 562 support a dynamic view on functional connectivity [95], since at any given time the 563 possible transitions depend on the current attractor (stable state). However, because 564 transitions within two different sets of attractors may be similar, similar dynamic 565 functional connectivity patterns, under a transition-centered view, may come from 566 distinct places in the landscape. Thus, clustering of dynamic functional connectivity 567 patterns may be difficult to interpret in dynamical system terms. Finally, although 568 the pattern of cross-attractor coordination can remain similar as the landscape 569 changes (Figure 9a-c), it is not the case for its energy cost ( Figure 11). This means 570 that the potential for exhibiting normal functional connectivity patterns may always 571 be there, but different structural constraints impose different energy costs. In our 572 model, the Default Mode Network dominates energy consumption over other parts of 573 the brain when the energy gaps are overall small. This resonates with the empirical 574 finding that aerobic glycolysis, accommodating small and rapid energy demands, is 575 significantly elevated in the Default Mode Network [96]. 576 The temporal diversity of the model brain is also affected by both local and 577 global structural constraints. In the local model, oscillatory activity requires a suffi-578 ciently strong excitatory-to-inhibitory connection. The oscillation may be damped 579 or sustained at various characteristic frequencies, contingent upon the strength of 580 excitatory-to-excitatory connection (see Figure 3a-d and Section S7 Oscillation). 581 The importance of inhibitory neurons and their interaction with pyramidal cells 582 for generating rhythmic activity has been well demonstrated in both theoretical 583 and empirical studies [97][98][99][100][101]. For multiple oscillatory processes to form complex 584 spatiotemporal patterns, it often requires the coexistence of diverse time scales [31, 585 38, 102]. In the present model, temporal differentiation can be caused by local 586 structural differences, i.e. the strength of local excitatory-to-excitatory connection 587 (Figure 13ab). It has been shown that incorporating such local structural diversity 588 in the reduced Wong-Wang model better describes real neural dynamics [9], demon-589 strating its empirical relevance. On the other hand, temporal differentiation can 590 also be induced solely by the structure of the global network-the whole defining 591 the parts (Figure 13cd). The diversity of node degree is a key contributor to the 592 spectral diversity in the low-frequency range (Figure 13d), which has been observed 593 empirically (e. g. [103]). It resonates with earlier findings that slow dynamics are 594 more reflective of the large-scale network structure (e. g. [104]). These multiscale 595 structural sources of temporal diversity may influence each other through their joint-596 action on brain synchronization and activity-based plasticity. Further theoretical 597 investigation of such cross-scale interaction may shed light on how structural and 598 dynamical properties stabilize each other across scales during brain development 599 (see [105]).

600
In conclusion, complex dynamic features such as multistability and temporal 601 diverse are both supported by local structural connectivity and further synergized 602 by the structure of the large-scale network. Cross-attractor coordination between 603 model brain regions in the multistable landscape provides a stronger explanation for 604 human functional connectivity than within-attractor or near-equilibrium coordination. 605 Energy costs of such coordination is further shaped by structural constraints across 606 scales and prominent in the Default Mode Network. These findings suggest that 607 empirically observed resting functional connectivity reflects transitions more than 608 stable brain states. Lowering the energy costs of such transitions may be a means to 609 restore normal functional connectivity. Constrained by the biophysics, the modeling 610 approach may be used to predict stimulation-induced transitions in experimental 611 and clinical settings. The computation of bifurcation diagrams (Figure 7) was carried out in MATLAB, 876 utilizing the build-in function fsolve. Given a proper initial guess, fsolve provides 877 the coordinates of a nearby fixed point and the Jacobian matrix at the fixed point. 878 The spectrum {λ k } 2N k=1 of the Jacobian matrix is used to classify the fixed point, 879 where N is the number of brain regions in the model. The fixed point is a stable 880 equilibrium if λ k is real and negative for all k. The fixed point is associated with 881 damped oscillation if Re λ k < 0 for all k and Im λ k = 0 for some k. The fixed 882 point is associated with a limit cycle if Re λ k > 0 and Im λ k = 0 for some k with 883 the additional criteria that after a small perturbation from the fixed point, the 884 time-average of the solution remains close to the fixed point. All other types of fixed 885 points are classified as unstable. For damped oscillation and limit cycles in the local 886 model, the frequency of the oscillation (Figure 3) is defined as | Im λ k |/(2π).

887
For the local model, a 2D dynamical system, the complete characterization of 888 all fixed points is relatively easy by searching exhaustively through a grid of initial 889 guesses (as for Figure 7a-c). This approach becomes unfeasible when it comes to 890 the global model due to the high dimensionality. Thus, for the global model, we 891 implemented a recursive search: for each value of G, (1) find zeros of equation 4-6 892 (main text) given a set of initial guesses that includes, if any, the zeros for G − δG 893 (δG = 0.01 for the present study); (2) sort the list of zeros obtained from (1) by 894 the average of S (i) E 's; (3) use the middle points between consecutive zeros in the 895 sorted list as initial guesses; (4) continue to use middle points between past initial 896 guesses as new initial guesses recursively until at least one new zero is found or the 897 recursion has reached a certain depth; (5) append the new zero(s) to the list of zeros 898 and repeat (2)-(5) until the number of identified zeros exceeds a certain value. In 899 the present study, we limit the maximal depth in (4) to 8 and the maximal number 900 of zeros in (5) to 200.

902
In the present study, we are interested in not only the high-frequency activity 903 measurable by, for example, EEG recordings but also low-frequency fluctuations 904 that are often a subject of investigation in fMRI studies. Therefore, we simulated 905 the BOLD activities induced by the underlying neural dynamics and examine their 906 low-frequency properties.

907
BOLD (Blood-oxygen-level-dependent) activities are computed using the Balloon-Windkessel model [1][2][3][4]. The hemodynamic response of the i th brain area takes the formṡ (S1) where the interpretation and value of the parameters are given in Table S1. The 908 initial condition is where ω N is the Nyquist frequency (approximately 0.7 Hz for the chosen sampling interval). The low-frequency power is defined as S3 Discretization of regional states 913 Although the dynamic landscape of the global model can be quite complex (Figure 7g-914 i), each region in the globally connected network still falls into discrete states 915 ( Figure S1) very much like in the local model (disconnected stripes in Figure 7a-c)-916 it is the combination of regional states that produces a great variety of attractors at 917 the global level. Discretized regional states (number on black disks in Figure S1) 918 thus give rise to discretized attractors in the global model. Rank correlation 919 (Spearman) between regional states across these discretized attractors are used to 920 quantify cross-attractor coordination in the model brain ( Figure 8b in the main text). 921 Using discretized attractors, we quantify how much two regions move up and down 922 together across attractors without considering the distance between the attractors. 923 The distance between attractors are considered separately as the energy cost that 924 constrains such transitions ( Figure 11 in the main text). The average functional connectivity is highly similar between Day 1 (Figure 8a in 927 the main text) and Day 2 ( Figure S2a). Linear regression analysis indicates that 928 the Day-1 matrix is highly predictive of the Day-2 matrix (red x in Figure S2b; 929 β 1 = 0.98, t(2143)=203, p<0.001, R 2 = 0.951; only elements below the diagonal are 930 compared due to the symmetry of the matrix). Moreover, the functional connectivity 931 matrix ( Figure 8a in the main text) obtained using Spearman correlation is highly 932 predictive of the corresponding Pearson correlation coefficients ( Figure S2c; β 1 = 1.03, 933 t(2143)=897, p<0.001, R 2 = 0.997), which itself is highly consistent across two days 934 (blue o in S2b; β 1 = 0.97, t(2143)=211, p<0.001, R 2 = 0.954). Figure S2: Average human functional connectivity highly reliable across two days and the types of correlation analysis. (a) shows the functional connectivity averaged across the same subjects as in Figure 8a (main text) but using the resting fMRI data (two runs) from Day 2. The matrices from the two days are highly correlated (b; red x). (c) shows the functional connectivity estimated using Pearson correlation averaged over all runs in Day 1, which is also consistent across two days (b; blue o), and is not markedly different from Figure 8a (main text). See text for statistical information. In Figure 13 of the main text, we illustrate with two simulated trials how high-939 frequency and low-frequency dynamics depend on local excitatory-to-excitatory 940 connectivity w EE and the topology of the global network. To show that these effects 941 are not incidental, we simulated 200 trials for each of the conditions: (1) the global 942 network is uniform but local connectivity w EE is diverse (as in Figure 13a,b), and 943

935
(2) local connectivity w EE is identical but the global network follows the human 944 connectome (as in Figure 13c,d). We characterize the high-frequency content of a 945 spectrum as its peak frequency, i.e. the frequency at which the spectral power is 946 the highest (e.g. peaks in Figure 13a,c); the low-frequency content as the integral 947 of the power between 0.01 and 0.1 Hz (Section S2). The dependency of these 948 features on local (w EE ) and global structural properties (node degree) is quantified 949 using Spearman correlation. The distributions of the correlation coefficients (ρ) and 950 corresponding p-values are shown in Figure S4 and Figure S5 for condition 1 and 2 951 respectively. Figure S4 shows that local structural connectivity w EE strongly affects 952 the peak frequency of the brain region (a,c) but not so much the low-frequency power 953 (b,d). The stronger the local connectivity, the higher the peak frequency. Figure S5 954 shows that the node degree of the global network has a strong and negative effect 955 on the low-frequency power, and a weak and positive effect on the peak frequency. 956 Figure 13 illustrates such dependencies using typical trials (median correlation 957 coefficients) from the distributions ( Figure S4-S5).
958 Figure S4: Dependency of peak frequency and low-frequency power on local excitatory-to-excitatory connectivity. 200 trials are simulated following the same parameter setting as Figure 13a,b, where the global network is uniform but the local connectivity w EE 's spread between 1 to 2 for different brain regions. The noise terms in equation 4-5 make these trials different realizations of the same noisy process. The peak frequency of the spectra, e.g. from 10 to 30 Hz in Figure 13a, strongly depends on local connectivity w EE (a: ρ's all close to 1; c: p-values all less than 0.05). In contrast, low-frequency power does not significantly depend on w EE (b: ρ's distribute around zero; d: p-values spread between 0 and 1). Figure 13b shows this lack of dependency in an example trial that corresponds to the median of the distribution (b). Figure S5: Dependency of peak frequency and low-frequency power on node degree in the global network. 200 trials are simulated following the same parameter setting as Figure 13c,d, where the local connectivity w EE 's are the same across brain regions but the global-network reflects the human connectome (see main text). The peak frequency of the spectra, e.g. between 0 and 30 Hz in Figure 13c, moderately increases with the node degree of each region (a: positive ρ's around 0.32; c: p-values all less than 0.05). Low-frequency power decreases more significantly with node degree (b: ρ's distribute around -0.6; d: p-values all less than 0.05). Figure 13d illustrates this dependency with an example trial that corresponds to the median of the distribution (b).
nullclines are confined within the interval (S11) Within this interval S I = f (S E ) (equation S9; red nullcline), overall, goes down 960 from +∞ to −∞, while S E = g(S I ) (equation S10; blue nullcline) goes up from −∞ 961 to +∞. This results from the dominant effect of H −1 p for a very large or very small 962 input.

963
In between these extremes, the effect of the linear term is more pronounced. This is especially the case for S I = f (S E ) (red nullcline): the linear term monotonically increases with S E , counteracting the descending trend of −H −1 E . Given a sufficiently strong excitatory-to-excitatory connection w EE (self-excitation), the linear term "twists" the nullcline counterclockwise, creating an ascending branch in the middle. If we balance the level of self-excitation with inhibitory-feedback-let w EE = w IEequation S9 becomes In this simplified case, increasing self-excitation w EE reduces the influence of H −1 E 964 such that the slope of middle branch approaches 1.

965
For S E = g(S I ) (equation S10), the linear term and the H −1 I term increase 966 together, so that S E = g(S I ) (blue nullcline) is always monotonically increasing. 967 Given a fixed w II , S E = g(S I ) increases with S I at an overall slower rate for larger 968 w EI , or more conveniently seen as S I = g −1 (S E ) increasing faster with S E for larger 969 w EI . Intuitively, increasing w EI twists S E = g(S I ) counterclockwise, seen as the 970 middle segment of the blue nullcline becoming more vertical. 971 We have discussed above how local connectivity w EE and w EI influence the gross 972 geometry of the nullclines-twisting the middle segment of the curve counterclockwise. 973 But how are these geometric changes going to affect the dynamics? We show below 974 that they critically control the multistability and oscillation in the local model.

984
Since g −1 (x) is always monotonically increasing and the existence of multistability 985 requires the existence of multiple intersections between g −1 (x) and f (x), a monoton-986 ically decreasing f (x) implies that the system cannot be multistable. In other words, 987 if the system is multistable, then f (x) cannot be monotonically decreasing.

988
This result highlights the importance of self-excitation w EE in equation S9 989 and equation S12-multistability can only occur when w EE is sufficiently large. 990 Correspondingly in the numerical result (Figure 3), the region of multistability 991 appears only for larger w EE 's.

992
Note that the above argument is not restricted to the present model, but applicable to models that share the geometry form of the Wilson-Cowan model in general. Nevertheless, one would hope to know how large a w EE is large enough for multistability to be possible, and this depends on the specific formulation of the transfer function (equation 3) and the underlying assumptions about neuronal level properties (equation 15). Ideally, to know the minimal w EE , one need to find the minimal slope of H −1 E (u(S E )) with respective to S E , where u(S E ) := S E / (τ E γ E (1 − S E )). The exact solution is, however, rather perplexing to calculate. Here we provide a rough, but simple, estimation instead. The slope of interest is (S14) Instead of finding the minimum of equation S14, we aim to find a representative 993 point S * E such that equation S14 is relatively small.

994
One option is to use the minimum of the numerator. The minimum of the numerator . By design, H E reaches its maximal slope a E at the inflection pointx, where H E (v) = r max /2. That is, we need But note here that, in the case where r max is a large number, the representative point 995 S * E is very close to one, which further results in a small denominator in equation S14 996 and a large slope for H −1 E . Thus, the inflection point of H E (v) is not a very good 997 choice.

998
To avoid the small denominator problem for equation S14, we need to choose a S * E as small as possible while H E (v) remains close to the line a E v − b E . For this purpose, we take v * to be the intersection between the line a E v − b E and the horizontal axis, (approximate values can be obtained from the Taylor expansion ofĤ E near v * ). Given equation S15, we need Now for the nullcline S I = f (S E ) to have a positive slope at S * E , one simply needs 999 w EE > h E . (S19) Here h E is approximately 0.2 based on the present parameter choices, inherited from 1000 Wong and Wang's initial derivation [7]. This result is confirmed numerically by the 1001 bifurcation diagrams (Figure 7a-c vs. Figure S6a) of the local model-multistability 1002 exists for some level of input I E when w EE > 0.2.

1003
Oscillation. Now we look for the conditions for oscillation to emerge. Here we 1004 are mainly concerned with the oscillation occurring on the ascending segment of 1005 S I = f (S E ) (red nullcline). Following a similar argument as Wilson and Cowan 1006 [8], one notice that for the flow around a fixed point-an intersection between the 1007 nullclines-to have consistent rotation, the nullcline g −1 (S E ) (blue) must have a 1008 greater slope than f (S E ) (red nullcline). Qualitatively, one would expect oscillation 1009 to be induced by increasing w EI , which twists g(S I ) (blue nullcline) counterclockwise. 1010 This expectation is confirmed by the numerical results in Figure 3a-d: oscillation 1011 emerges for sufficiently large w EI for fixed points on the ascending branch of 1012 S I = f (S E ) (Figure 4a-d). 1013 Quantitatively, we consider the derivative of the two nullclines at a respective representative point. First, we extend the results in equation S17-S18 to the second nullcline S E = g(S I ) (blue): For parameters used in the present study, h I ≈ 0.4. We have the slope of the two nullclines at their respective representative points, and we need With balanced inhibitory feedback w IE = w EE , we have For very large w EE , one simply need 1014 w EI > w II + h I .
Given the present parameter choices, we need w EI > 0.45 to induce oscillation for 1015 some level of input I E and I I . This is in line with the numerical results in Figure 3. 1016 For h E > 0, as assumed here, lowering w EE also lowers the threshold for oscillation. 1017 Linear stability analysis. In addition to the presence of oscillation, one would also want to know if such oscillation is sustainable or damped. Here we extend the above analysis by linearizing the system near a specific fixed point. A fixed point is where the two nullclines (equations S9-S10) intersect. Conveniently, we let them intersect at their respective representative points (S * E , f (S * E )) and (g(S * I ), S * I ) (see equation S17 and equation S20), The two equations can be satisfied by the appropriate choice of I E and I I . The fixed point of our choice (S * E , S * I ) inherits a couple of properties from the above analysis, which we shall soon see. First, we define For simplicity, let parameters for p ∈ {E, I}. Note that by definition, both α p and β p are positive. Given 1018 parameters used in the present study, we have α E ≈ 14, α I ≈ 111, β E ≈ 71, and 1019 β I ≈ 276. 1020 We write the Jacobian matrix as The eigenvalues of the Jacobian are Assuming that the system is already oscillatory near the fixed point, i.e. tr 2 J < 4 det J, to have sustained oscillation (limit cycle), we need Given the parameters used in the present study, the emergence of limit cycles requires 1021 w EE > 2. Correspondingly in the numerical results (Figure 3), equation S50 provides 1022 an estimate of the lower bound of the Hopf bifurcation (gray dashed line). Note 1023 that stronger inhibitory-to-inhibitory connection w II increases the minimal w EE 1024 required to induce sustained oscillation. Overall, these analyses show that sustained 1025 oscillation requires both strong self-excitation and a sufficiently active inhibitory 1026 population.

1028
In summary, we have shown analytically how structural connectivity w EE and 1029 w EI critically shape the dynamics-in this very low-dimensional parameter space, 1030 the system can easily switch between qualitatively different behavior. In particu-1031 lar, excitatory-to-inhibitory connectivity w EI controls the emergence of oscillation; 1032 excitatory-to-excitatory connectivity w EE controls both the emergence of multistabil-1033 ity and sustained oscillation. The qualitative description of the system only depends 1034 on the gross geometric form of the Wilson-Cowan model, but the exact bound-1035 aries between regimes depend on the specific transfer function and the associated 1036 biophysical constraints.   Multistability. Following a similar argument as for the local model, we first show 1046 that, without global interaction (i.e. G = 0), the system cannot be multistable, if 1047 w EE is sufficiently small such that f (i) ( S E ) monotonically decreases with S (i) E for 1048 all i. As shown above, the monotonicity condition implies that each local node by 1049 itself is not multistable.

1050
Proof. Assume there are at least two distinct fixed points of the system: S * and S * * , 1051 where S = ( S E , S I ) and S p = (S (1) p , · · · S (i) p , · · · S (N ) Figure S6: Synergistic memory between monostable nodes. Three bifurcation diagrams are shown for local parameters w EE = 0.1 and w EI = 0.35. They correspond to Figure 7a, d, g but with a lower w EE such that each local node by itself is monostable for any level of input (a). While each local node is completely monostable (no memory capacity), once there is sufficient global coupling G between them, the whole brain acquires memory capacity (b, c) that cannot be attributed to the parts alone-synergistic memory. Nevertheless, the size of the global memory capacity is still fundamentally constrained by the complexity of the local node (42 attractor branches in (c), very small compare to Figure 7g, h, i). See text for further discussion.
What we have not addressed in the above analyses is to what extent the global 1076 system is multistable-what is the number of stable states, or the size of the memory 1077 capacity-and what are the contributions from local self-excitation and global 1078 network connectivity. An analytical approach to this problem is difficult; thus, it 1079 is mainly addressed numerically (c.f. Figure S6 and Figure 7). Nevertheless, we 1080 provide an intuitive argument below as to how local and global connectivity affects 1081 the relevant geometrical properties of the dynamical system.

1082
Local origin of geometrical complexity. At an intuitive level, the number of 1083 intersections between these hypersurfaces (equation S54-S55) is likely to increase with 1084 the number of folds of each surface. In the present case, the folding of hypersurfaces 1085 entails the temporary reversal of the sign of its partial derivative along a certain 1086 direction. Observe equation S54 and see that global coupling cannot create any 1087 folding of the surfaces. Thus, the geometrical complexity of the nullclines purely 1088 depends on the local properties of each node, in particular, the folding effect of 1089 self-excitation w way dependent on the structure connectivity C ij . This may remove or introduce new 1093 intersections between the surfaces without changing the geometrical complexity of 1094 these surfaces. Thus, global coupling allows system-level multistability to be created 1095 synergistically, given appropriate structural connectivity C ij .

1097
In summary, local and global coupling produce different geometrical effects on 1098 the system and jointly affect the number of possible stable states.