Cross-attractor repertoire provides new perspective on structure-function relationship in the brain

The brain is a complex system exhibiting ever-evolving activity patterns without any external inputs or tasks. Such intrinsic dynamics (or lack thereof) are thought to play crucial roles in typical as well as atypical cognitive functioning. Linking the ever-changing intrinsic dynamics to the rather static anatomy is a challenging endeavor. Dynamical systems models are important tools for understanding how structure and function are linked in the brain. Here, we provide a novel modeling framework to examine such structure-function relations. Our deterministic approach complements previous modeling frameworks, which typically focus on noise-driven (or stochastic) dynamics near a single attractor. We examine the overall organizations of and coordination between all putative attractors. Using our approach, we first provide evidence that examining cross-attractor coordination between brain regions could better predict human functional connectivity than examining noise-driven near-attractor dynamics. Further, we observed that structural connections across scales modulate the energy costs of such cross-attractor coordination. Overall, our work provides a systematic framework for characterizing intrinsic brain dynamics as a web of cross-attractor transitions and associated energy costs. The framework may be used to predict transitions and energy costs associated with experimental or clinical interventions.

functional connectivity than examining noise-driven near-attractor dynamics. 23 Further, we observed that structural connections across scales modulate the 24 energy costs of such cross-attractor coordination. Overall, our work provides a 25 systematic framework for characterizing intrinsic brain dynamics as a web of 26 cross-attractor transitions and associated energy costs. The framework may 27 be used to predict transitions and energy costs associated with experimental 28 or clinical interventions. 29 1 Introduction 30 A fundamental goal of neuroscience is to understand how the structure of the brain 31 constrains its function [1]. The advent of neuroimaging techniques has enabled 32 detailed, quantitative examination of the structure-function relation, often by com-33 paring the structural and functional connectivity between brain regions [2, 3]. The 34 resting-state functional connectivity is of particular interest, given its relevance to a 35 variety of cognitive functions, neurological diseases, and psychiatric disorders [3,4]. 36 Though structural connectivity and functional connectivity are linearly correlated, 37 the former does not entirely predict the latter-strong functional coupling exists 38 between regions with only weak or indirect structural connections [2,5]. Nonlinear 39 dynamical models have been used to provide a mechanistic understanding of the 40 structure-function relation [6-8] and provided many insights (e.g. [9][10][11][12][13][14][15][16]). While 41 successful, previous dynamical system approaches open focus on dynamics near a 42 single stable state, i.e., an attractor. However, biological systems such as the brain 43 are often multistable [17][18][19], i.e., multiple attractors can coexist in the brain's 44 dynamical landscape. Such multistability begs the question of whether examining 45 the overall layout of brain's attractor states could better inform or complement what 46 we know about the structure-function relation than explorations around a single 47 state. 48 Intrinsic brain dynamics have long been observed [20, 21] but often treated 49 as a baseline subtracted from task-positive activities. This baseline, however, is 50 more active than meets the eye: it consumes the largest fraction of the brain's 51 energy resources, while task-related consumption adds little [22]. It constrains 52 task performance and related neural activities across multiple time scales [23-53 25], and sustains alteration in neurological and psychiatric disorders [26][27][28][29][30]. In 54 contrast to the restless dynamics is the (relatively) static structure-the anatomical 55 connections between brain regions, which can be estimated non-invasively using 56 large-scale tractography on diffusion-weighted images [31]. How can one compare 57 the ever-changing with the unchanging? From a statistical perspective, one may 58 compute the time-averaged features of the dynamics, such as the correlation between 59 signals generated by two brain regions across time-a common measure of functional 60 connectivity. Such functional connectivity patterns can be directly compared to 61 structural ones through linear correlation [2]. From a dynamics perspective [6,8,32], 62 the strength of anatomical connections can be incorporated as constant parameters 63 in a system of differential equations, i.e. a dynamical system. The dynamical system, 64 in turn, describes how the state of a model brain, endowed with realistic anatomy, 65 would evolve over time. The time series generated by the model brain and the 66 derived functional connectivity patterns can then be fitted to that of the real brain. 67 Thus, a dynamical system naturally bridges between the unchanging structure and 68 the ever-changing dynamics. 69 The brain's nonlinearity gives rise to multistability and sustained oscillation [17-70 19, 33-36]. In a multistable landscape, multiple attractors coexist. That is, a brain 71 may persistent in one of many qualitatively distinct patterns of activity depending 72 on the environment and intrinsic factors. One popular modeling approach is to 73 simulate the noise-driven dynamics near a chosen attractor, such as the low activity 74 ground state, and compare it to the human resting brain dynamics (see [37] for a 75 summary of different approaches). Such noise-driven exploration of a single attractor 76 has been shown to exhibit key features of human resting brain dynamics, especially 77 near criticality (e.g. [10,12,16]). On the other hand, noise-driven exploration 78 beyond a single attractor-across multiple attractors or "ghost" attractors-has 79 been shown to capture non-stationary resting brain dynamics and the switching 80 between different dynamic functional connectivity patterns [11,14,15,38]. The best 81 fit to empirical data is often found near the onset of multistability [14,38]. These 82 observations suggest that examining the layout of the attractor repertoire over the 83 entire multistable landscape could be crucial for understanding the organization of 84 resting brain dynamics (c.f. [14]). 85 Complementing existing single-attractor approaches, the present work focuses on 86 the deterministic features of the multistable landscape and examines their empirical 87 relevance. Specifically, we systematically study the organization of the attractor 88 repertoire as a window into the overall shape of the dynamic landscape. To do so, 89 we use a biophysical network model that formally combines the reduced Wong-Wang 90 model [12,13,39] and the Wilson-Cowan model [40,41]. We first should that 91 the model exhibits extensive multistability, i.e. a large repertoire of attractors 92 to serve as landmarks of the landscape. Further, the model also allows us to 93 examine, computationally and analytically, how structural features across scales 94 shape this repertoire. It is important to note that here we used a broader definition 95 structural features and not only include large-scale structural connectivity between 96 brain regions, but also local recurrent connectivity within regions, and biophysical 97 constraints at the cellular level. Using this modeling framework and a small dataset 98 from the Human Connectome Project (HCP; n=100), we provide evidence with 99 regards to how the cross-attractor relations in the repertoire could better capture 100 key features of human resting functional connectivity and how such features are 101 shaped by structural features across scales. Finally, we provide a novel framework to 102 analyze the energy constraints for such cross-attractor coordination across different 103 local and global structures. Whole-brain dynamics are modeled as the mean-field activity of neuronal populations 107 in each brain region. We use an adapted version of the Wong-Wang model [13, 108 39] with a sigmoidal transfer function (equation S11). The adaption improves 109 the biological plausibility and multistability upon the original model. Here, we 110 briefly introduce the model; an extensive analysis of the numeric and mathematical 111 properties of the model is provided in the Supplementary Materials (Section S4-S6 112 for numeric results, Section S13-S14 for analytical results). Each model region 113 contains a pair of excitatory (E) and inhibitory (I) populations, whose activity is 114 described by the local model ( Figure 1a Figure 1: 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). The two populations are coupled 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 via long-range structural connections (red dashed curves). (b) Nodes in the global model correspond to 66 anatomical regions of the human brain, which can be linked together by the human connectome (see text). Regions are indexed from 1 to 66 (1-33 on the right hemisphere, 34-66 on the left hemisphere in reverse order, following [12]). Specific region names are listed in Table S1.
In the present work, the local and global models are used in two ways: (1) to 130 compute the repertoire of attractors using zero-finding algorithms and (2) to be 131 numerically integrated to generate simulated brain dynamics. The former is used to 132 characterize the overall organization of the model dynamic landscape. The latter 133 is used to characterized local explorations of the dynamic landscape. Both aspects 134 are compared to the human data to demonstrate the empirically relevant features. 135 Below, we first illustrate the concept of a dynamic landscape and cross-attractor 136 coordination using a toy example, which is followed by more realistic models to 137 examine how structural properties across scales affect the dynamic landscape. 138 2.2 Multistable landscape of the brain shaped by structural 139 properties across scales 140 Here, using a toy example (Figure 2), we first intuitively illustrate the idea of a 141 dynamic landscape and how it may serve as a system-level description of intrinsic 142 brain dynamics. An attractor in the toy model represents a stable pattern of 143 activation over the whole brain, shown as boxed brains in Figure 2. The dynamic 144 landscape determines a repertoire of attractors with possible paths of transitions 145 between them. Such a landscape is akin to the topographical map of a land, with 146 hills and valleys, where the attractors could be presented as valleys while the paths 147 between the valleys can be thought of as transitions. Figure 2a represents such a toy 148 landscape with four attractors (i-iv), each with a different whole-brain activation 149 map. The landscape does not have to be static -as the landscape changes, some 150 attractors (valleys) could be destroyed or created, causing a discrete change of the 151 attractor repertoire, a.k.a. a bifurcation. Further, changes in the attractor repertoire 152 can alter how the brain regions coordinate with each other. Two example bifurcations 153 are shown in Figure 2b-c. In Figure 2b, only two out of four attractors are left, 154 such that the brain can now only transition between attractor (i) and (iii), thereby 155 leading both hemispheres to be in sync (on or off together). This coordination of 156 brain regions (or hemispheres) can be captured by estimating the cross-attractor 157 coordination matrix (shown in Figure 2). Similarly, in Figure 2c, three out of four 158 attractors are left after bifurcation, leading to more complex coordination between 159 brain regions (or hemispheres). Figure 2d-f presents a slightly more complex example, 160 where each brain region can now take three activation values (instead of just on or 161 off), resulting in more complex spatial activation patterns across the whole brain as 162 well as complex coordination between brain regions. See Section 4.3 for more details 163 about how the coordination matrix is estimated.
164 Figure 2: 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. Each attractor corresponds to a distinct pattern of activation over the whole brain (i-iv). Influenced by external input or intrinsic noise, the model brain may transition from its current state (attractor i, bright purple ball) to a different one (ii, iii, or iv, dim purple balls), indicated by black arrows. Structural features of a model brain can alter the shape of the landscape, causing some attractors to appear or disappear through a process mathematically named a 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. For example, in landscape (a), the left and right hemisphere can be co-activated during a transition (i→iii), or activated separately through other transitions (i→ii, or i→iv). In contrast, in landscape (b), the left and right hemisphere can only be co-activated, and in (c), only activated separately. Numerically, a repertoire of attractors can be represented as a matrix, where each row represents an attractor and each column represents a brain region (repertoire matrix a, b, c, with entries shown as blue/red color blocks). The overall inter-regional coordination across attractors can be estimated by the 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 (coordination matrix a, b, c, shown to the right of each repertoire). In more complex landscapes (not shown), there are many more attractors, and they correspond to subtler patterns of activation (d,e; see also Figure 3). The coordination between brain regions during a transition is correspondingly more complex (f=e-d), with some regions co-activated (red) while others co-deactivated (blue).
Next, we present a set of more realistic examples to show how local as well as 165 global structural connectivity of the brain can shape the dynamic landscape, its 166 associated repertoire of attractors and their transitions. Here, we depict dynamic 167 landscapes and their changes as bifurcation diagrams (Figure 3; see Section 4.2 for 168 computational details). Figure 3 shows nine different bifurcation diagrams, across 169 its three rows and columns. The rows correspond to bifurcation diagrams from: 170 (first row: a-c) a single brain region (local model); (second row: d-f) the entire brain 171 with uniform connectivity across all brain regions; and (third row: g-i) the entire 172 brain with realistic connectivity across all brain regions. Thus, the second and third 173 rows of Figure 3 aims to depict the effect of changes in global structural connectivity 174 on the dynamical landscape. The columns, on the other hand, in Figure 3 aims to 175 depict the effect of changes in local connectivity, i.e., level of excitation within the 176 individual brain regions, on the dynamic landscape.

177
In each bifurcation diagram, the y-coordinate of each colored point indicates 178 the position of an attractor: here we use the average activity of all excitatory 179 populationsS E (Figure 3d-i; black points are repellers). The x-coordinate indicates 180 the value of a control parameter, which modulates the shape of the underlying 181 dynamic landscape: here we use the overall strength of long-range connections-the 182 global coupling G (equation 6). Further, each vertical slice of a bifurcation diagram 183 contains the repertoire of attractors and repellers in a fixed landscape (an example 184 slice is shown in Figure 3h), corresponding to stable and unstable patterns of brain 185 activity respectively. Lastly, an attractor traces out a horizontal "stripe" as it changes 186 continuously with the landscape (examples shown in Figure 3d-i). A colored stripe 187 merges with a black stripe at a bifurcation, where an attractor is annihilated by a 188 repeller. The number of colored stripes indicates the complexity of the landscapes, 189 i.e., more stripes indicate more attractors.

190
Locally wihtin each model brain region, the dynamics are controlled by local 191 structural connectivity (w's in Figure 1a; see Section S4 for detailed numeric results 192 and Section S13 for analytical results). In particular, a single model region can 193 switch between a rich set of dynamic regimes by varying the excitatory-to-excitatory 194 connections (w EE ) and the excitatory-to-inhibitory connections (w EI ; Figure S1). 195 Figure 3a-c show the bifurcation diagrams for three local connectivity settings from 196 three distinct dynamic regimes (regime e, d, and a respectively in Figure S1). With 197 an overall increase of local excitatory connectivity from (a) to (c), a single region 198 becomes more complex-more attractors and stronger oscillatory activities.

199
To understand the effects of global connectivity, we first examined the brain 200 dynamical landscape where all regions are uniformly connected with each other 201 (Figure 3d-f). Here, stronger local excitatory connections (e,f) produce a more 202 complex landscape (3 attractor stripes) than weak ones (d; 2 attractor stripes; see 203 Figure S10 for even weaker local connections). These bifurcation diagrams are very 204 similar to those of a single brain region (Figure 3a-c), in terms of the number of 205 attractors and the presence of oscillation. In fact, the whole brain (Figure 3e) moves 206 up and down together between discrete states of activation, very much like a single 207 region ( Figure 3b). Note that, for the global model (equation 4-6) to be multistable, 208 a minimal amount of global coupling is required, i.e., G > 1 (Figure 3d-f). If the 209 brain regions act independently (G = 0), both the individual regions (Figure 3a-c at 210 I E = 0) and the whole brain (Figure 3d-f at G = 0) are monostable-there is only 211 one stable pattern of activity, where the gating variables are all close to zero. This 212 result indicates that a functionally complex brain can emerge out of the synergistic 213 interaction between simple regions. Additional analytical and numerical results are 214 provided in Section S14 (Multistability) for further validation and generalization. 215 Next, we show that in addition to the global coupling G, the details of inter-216 regional connections matter too (C ij in equation 6). Given a realistic global structural 217 connectivity (human connectome; Figure 3g-i), the complexity of the whole-brain 218 dynamic landscape increases dramatically: 171 attractor stripes in (g), 610 in 219 (h), and 682 in (i) (per single-linkage clustering). Correspondingly, the patterns 220 of activation (Figure 3h) are more complex, with greater differentiation between 221 regions; the coordination between brain regions across attractors is consequently 222 more flexible and subtle (as depicted in Figure 2f). The heterogeneous nature 223 of the human connectome breaks the large-scale spatial symmetry of the model 224 brain, creating more functional differentiation between brain regions and greater 225 functional complexity for the whole brain. In short, the complexity of the global 226 dynamical landscape is a joint product of strong local excitatory connection and 227 complex topology of the large-scale network. See Section S14 for additional analytical 228 supports.

229
Figure 3: Local and global structural properties jointly determine the complexity of whole-brain dynamics. (a-c) show the bifurcation diagrams of the local model for three different types of 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 (representatives of distinct dynamic regimes of the local model, Figure S1). 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 an 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 (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 a 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 4b). See Figure 2 for a cartoon illustration of attractor repertoires and the associated cross-attractor coordination matrices.
2.3 Cross-attractor coordination reveals large-scale symme-230 try of human brain functional connectivity 231 In this section, using data from the Human Connectome Project (HCP) [43] we 232 present both qualitative and quantitative results to show how cross-attractor coordi-233 nation could better relate with key features of human resting functional connectivity 234 than noise-driven within-attractor coordination. For qualitative results, we used 235 averaged structural and functional connectivity across subjects from a smaller HCP 236 cohort (n=11 individuals), whereas for the quantitative results, we used individual 237 structural and functional connectivity estimates across all 100 unrelated individuals 238 from the HCP data.

239
Typically, human functional connectivity (FC) is calculated from the fMRI data 240 by estimating co-fluctuations across brain regions. The estimated resting state FC 241 matrix usually reflects large-scale symmetry across the two hemispheres, such that 242 brain regions across the two hemispheres highly co-fluctuate (i.e., high antidiagonal 243 values in Figure 4a). For qualitative results, we used averaged structural and 244 functional connectivity matrices across a small subset of the HCP n=11 unrelated 245 individuals (see Section 4.5 and Section S9 for more details). Using averaged 246 structural connectivity matrix, a dynamical landscape (Figure 3h) was generated 247 and cross-attractor coordination was estimated for a chosen G = 2.2. At the 248 selected G, 97 attractors were found. Mathematically, each attractor is represented 249 by a row vector denoting the activity level of all brain regions (equation 8-9). 250 Thus, for a chosen G, using the (#attractors × #regions) matrix we estimated 251 the cross-attractor coordination between regions (equation 10), such that regions 252 that co-fluctuate across attractors tend to show high cross-attractor coordination 253 (or similarity). See Section 4.2 for mathematical details and Figure 2a-c for an 254 intuition regarding the estimation of cross-attractor coordination matrix. As shown 255 in Figure 4b, the dominant feature of human resting state FC, i.e., large-scale 256 symmetry across brain regions, is well preserved in the attractor cross-coordination 257 matrix. Such inter-hemispheric symmetry is not seen in stochastic within-attractor 258 coordination ( Figure 4d). Moreover, the pattern of within-attractor coordination is 259 a closer reflection of the human structural connectivity (Figure 4c) than the human 260 functional connectivity (Figure 4a). 261 Figure 4: Cross-attractor coordination captures large-scale symmetry of human functional connectivity and its nonlinear dependency on structure better than within-attractor coordination. An example of human functional connectivity matrix (a) is calculated using the resting fMRI data from the Human Connectome Project [43], averaged over 11 unrelated subjects. The average structural connectivity of the same subjects is shown in (c). Regions (columns and rows) are ordered symmetrically for the left and right hemispheres (see Figure 1b) to reveal the large-scale symmetry of resting brain dynamics (ordering follows [12]; see the complete list of region names in Table S1). White dashed 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 patterns within the hemispheres are similar to each other and similar to inter-hemispheric connectivity patterns. The symmetry between intra-and inter-hemispheric connectivity is well captured by inter-regional coordination in the model brain across attractors (b) (w EE = 2, w EI = 1, G = 2.22 for local maximum model-human similarity; c.f. Figure 3h gray slice). Such symmetry is not captured by the coordination within any of the said attractors (d; best fit within-attractor coordination matrix). When examined across individuals (using n=100 unrelated HCP participants), a significantly better model-human similarity was obtained when crossattractor coordination was used instead of within-attractor coordination (e); the difference between cross-attractor coordination and within-attractor coordination with respect to model-human similarity is even greater when partial correlation is used to control for the contribution of structural connectivity (f). Comparing intra-hemisphere and inter-hemisphere functional connectivity separately (g), we found that cross-attractor coordination captures human intra-and inter-hemisphere functional connectivity equally well (red bars), while within-attractor coordination is better at capturing intra-hemisphere than inter-hemisphere functional connectivity (blue bars). The distributions of correlation coefficients were obtained through a model fitting procedure with the same local parameters (w EE = 2, w EI = 1) while allowing the G parameter to vary from 1.7 to 3.0 by steps of 0.1.
Next, to examine the model fit on an individual basis, we ran a quantitative anal-262 ysis using the n=100 unrelated-individuals cohort of the HCP data [43]. Our results 263 confirm that cross-attractor coordination can better predict human functional connec-264 tivity than within-attractor coordination. The maximum model-human correlation 265 for cross-attractor coordination has an average Spearman's Rho of 0.50 (0.09), while 266 the correlation for within-attractor coordination has an average of 0.34 (0.06) (Fig-267 ure 4d). Based on within-subject paired t-test, cross-attractor coordination provided 268 significantly better fit of FC (t(99) = 20.2, p < 10 −36 ; Figure 4e). To understand how 269 the two types of models relate to the underlying structural connectivity, we calculate 270 the partial correlation between model coordination matrices and human functional 271 connectivity, controlling for the linear contribution of structural connectivity. With 272 the contribution of structural connectivity controlled, the model-human correlation 273 for cross-attractor coordination has an average Spearman's Rho of 0.39 (0.10), which 274 is significantly greater than that of the within-attractor coordination 0.17 (0.07) 275 (Figure 4f; t(99) = 21.9, p < 10 −39 ). Furthermore, the difference between the partial 276 correlation coefficients for the cross-attractor coordination and the within-attractor 277 coordination (0.22 ± 0.10) is significantly greater than that of the regular correlation 278 coefficients (0.16 ± 0.08) with t(99) = 9.5 and p < 10 −14 (Figure 4f vs. Figure 4e), 279 which suggests that within-attractor coordination is more linearly dependent on 280 the structural connectivity. Finally, we show that the model-human similarity for 281 inter-hemisphere coordination is significantly lower than that of intra-hemisphere co-282 ordination for within-attractor coordination but not for cross-attractor coordination 283 (Figure 4g). In other words, unlike within-attractor coordination, cross-attractor 284 coordination captures human intra-and inter-hemisphere functional connectivity 285 equally well.  In the above section, we computed the cross-attractor coordination matrices (equa-297 tion 10 in Section 4.3), which only concerns whether two brain regions move up and 298 down together across the dynamic landscape, but not how difficult or metabolically 299 expensive such movements are. In this section, we examine the "energy gaps" between 300 the attractors, and how they are shaped by local and structural properties of the 301 model. Figure 5a gives a conceptual illustration of the relation between attractors 302 and the energy gaps between them. Each attractor is associated with an average 303 level of activity or energy (S E ; equation 11) given a fixed parameter G. There is an 304 energy gap between each pair of adjacent attractors (equation 12; see a full technical 305 description in Section 4.3). In the subject-average model (n=11, HCP; Figure Table S5-S6 for corresponding confidence intervals from the permutation tests. Error bars are standard errors.) Next, we examine the effect of energy constraints on model fit. To incorporate 312 the effect of energy constraints (see Section 4.3 for detailed methods), we split an 313 ordered attractor repertoire into sub-repertoires between which the energy gap is 314 considered too high. Thus, we obtain a sub-repertoire above the maximum energy 315 gap, equation 13, and a sub-repertoire below the maximum energy gap, equation 14. 316 Cross-attractor coordination matrices computed within the sub-repertoires can 317 be considered as energy-constrained coordination patterns between brain regions. 318 Figure 6a-c shows that such energy-constrained cross-attractor coordination (dashed 319 lines) is more sensitive to different structural features in its ability to capture 320 human functional connectivity. The energy constraint inflicts a greater loss of 321 model-human similarity when the local structural connectivity is weak (area of the 322 shaded region shrinks from Figure 6a to c, and bars in d decreases with strong 323 local connectivity from left to right) and the global structural connectivity is strong 324 (height of shaded regions grows with G in Figure 6a-c). The loss of similarity grows 325 with the maximum gap size (Figure 6d; ρ = 0.96, p<10 −100 for w EE = 0.7 and 326 w EI = 0.35; ρ = 0.92, p<10 −100 for w EE = 2 and w EI = 1; ρ = 0.85, p<10 −100 327 for w EE = 2.8 and w EI = 1). In short, local and global structural connectivity 328 both influence the energy costs associated with cross-attractor coordination, and 329 thereby the model-human similarity under energy constraints. When cross-attractor 330 coordination matrices were fitted to n=100 unrelated HCP subjects (Figure 4e), 331 the optimal G is 2.5 ± 0.28 (Figure 6e). The corresponding maximum energy gaps 332 average to 0.12 ± 0.08 (Figure 6f), and the corresponding mean energy gaps average 333 to 0.009 ± 0.015 ( Figure 6f). Thus, the optimal cross-attractor models in the present 334 study are located in the regime least affected by energy constraints (Figure 6d, 335 G<0.2). 336 Figure 6: The loss of model-human similarity due to energy constraints depends on local and global structural features. (a-c) The model-human similarity for cross-attractor coordination (no energy constraint, black solid lines) is stable with respect to varying global coupling G and local excitatory connectivity w EE and w EI (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). Dashed lines indicate the model-human similarity when the model is energy-constrained, i.e., it does not cross the maximum energy gap between attractors. This "energy-constrained" similarity is computed by splitting the attractor repertoire into two subrepertoires, one above the maximum energy gap and one below it (c.f. Figure 5a). The cross-attractor coordination within each subrepertoire is compared to the human functional connectivity through Spearman's correlation. The dashed lines indicate the greater correlation coefficient (ρ) among the two subrepertoires. The shaded area (∆ρ) indicates the loss of model-human similarity due to the energy constraint. (d) Each point in the scatter plot represents the percent loss of model-human similarity for not crossing the maximum energy gap, given a specific combination of global coupling G and local connectivity w EE and w EI . Overall, the loss increases with the maximum gap size, which in turn depends on G (Figure 5b). The average loss (bars in d) decreases with increasing local excitatory connectivity (w EE , w EI ). When the cross-attractor model was fitted to individual subjects (using n=100 unrelated HCP participants; c.f. Figure 4e), the optimal global connectivity G (e) and the corresponding maximum energy gap (f) and mean energy gap (g) are low, where the cross-attractor coordination is least affected by energy constraints. (*** p<0.001 with Bonferroni correction; see Table S4 for corresponding confidence intervals from the permutation tests. Error bars are standard errors.)

337
The present work examines how the brain's multistable dynamic landscape can 338 be shaped by structural features across scales and what features of the landscape 339 are relevant to empirical observations. Complementing the previous stochastic 340 noise-driven exploration approach, the present work focuses on the deterministic 341 features of the multistable landscape and examines their empirical relevance. We 342 demonstrate that large-scale symmetries of human functional connectivity patterns 343 and their nonlinear dependency on the structure could be better explained by the 344 relation between attractors in the landscape than the property of any individual 345 attractor. Thus, the present work offers a novel cross-attractor perspective on resting 346 brain dynamics, equipped with a computational framework to produce empirical 347 relevant summaries of the attractor repertoire in full as well as in parts.

348
The functional complexity of the model brain is controlled by both local and 349 global structural connectivity. At the level of a single isolated brain region, the 350 dynamic repertoire can be effectively controlled by two key local structural properties: 351 local excitatory-to-excitatory connectivity (self-excitation) and local excitatory-to-352 inhibitory connectivity. In the real brain, local excitatory-to-excitatory connections 353 are particularly abundant [45], and in the model brain, they contribute indispensably 354 to multistability (Section S13). Multistability is a key source of biological complexity 355 from molecular to social levels [18,19], often tied to self-excitation or positive feedback 356 In the present work, this critical value is a constant (equation S29), which depends 365 on cellular-level properties such as the membrane time constant and the gain of the 366 input-output response. Thus, manipulating such microscopic properties can induce, 367 or remove, multistability from a single brain region.

368
At the large-scale network level, multistability can be created or amplified by 369 the synergistic interaction between mono-or multi-stable brain regions. Different 370 large-scale network structures have dramatically different capabilities at amplifying 371 local complexity: a realistic global network is much more powerful than a uniform 372 one. The human connectome breaks the spatial symmetry of the global model, 373 whereas symmetry breaking is often a key to complex dynamics [33, 54-58]. On the 374 other hand, the human connectome is endowed with more specific features such as 375 modularity, small-worldness, and multiscale characteristics [1,[59][60][61]. A systematic 376 study of how these features alter the geometry of the global dynamic landscape is 377 worthy of further theoretical investigation (see Section S14).

378
Within the multistable landscape sculpted by the human connectome, inter-379 regional coordination across attractors exhibits key features of human functional 380 connectivity patterns. Such cross-attractor coordination better captures human 381 functional connectivity than within-attractor coordination-synchronization between 382 brain regions within the same basin of attraction. This finding raises the possibility 383 that functional connectivity patterns reflect transitions between stable brain states 384 more than the brain states themselves. A transition-based, or cross-attractor, view 385 on functional connectivity has several theoretical and empirical implications. 386 First, it provides an explanation for the large-scale symmetry of human FC, i.e., 387 the similarity between intra-and inter-hemispheric connectivity patterns. It has 388 been noted that within-attractor dynamics of similar models lack such symmetry, 389 exhibiting weak inter-hemispheric FC [12,16]. The weak inter-hemispheric FC has 390 been attributed to an underestimation of structural connections across hemispheres 391 using diffusion-weighted imaging. This explanation is reasonable given that within-392 attractor dynamics can be approximated by a linear dynamical system, which 393 closely depends on the structural connectivity [12,13]. Nevertheless, an alternative 394 explanation could be that human FC implicates far-from-equilibrium dynamics where 395 the nonlinearity cannot be ignored [15]. It is a signature of nonlinear systems that a 396 small input does not necessarily produce a small effect. Indeed, strong functional 397 connectivity in humans is known to exist between regions that are not directly 398 connected [5]. As we have shown, cross-attractor coordination takes into account 399 such nonlinear effects. Mathematically, our cross-attractor approach amounts to 400 studying the relation between the zeros of a nonlinear function-a problem that 401 does not admit a linear approximation. The symmetry between intra-and inter-402 hemispheric connectivity is likely to reflect a symmetry of the set of all zeros, i.e. the 403 (roughly) invariance of the zero set under the exchange of variable indices between 404 the homologous regions of the left and right hemispheres. The invariance of the zero 405 set is, in turn, consequent to the invariance of the differential equations under such 406 a left-right reflection. A full mathematical treatment of the problem is beyond the 407 scope of the present study. Nevertheless, it invites theoretical investigations of the 408 symmetry groups of nonlinear neural dynamical models (c.f. [55]).

409
The second implication is that the cross-attractor view is compatible with 410 treating functional connectivity as both stationary and dynamical. Cross-attractor 411 coordination, when measured over the entire dynamic landscape, is itself time-412 invariant. Empirically observed stability and convergence of human functional 413 connectivity [62, 63] may reflect this invariance of the underlying landscape. The 414 static landscape can also support dynamic functional connectivity (dFC) [15,64]. At 415 any given time, the possible transitions depend on the attractor currently dwelled 416 upon. Thus, in a short time window, cross-attractor coordination is confined to 417 a subset of attractors. In this perspective, dFC reflects the transitions between 418 attractors in a subset of the repertoire. As a consequence, a state-based and a 419 dFC-based representation of neural dynamics may diverge-subsets of attractors 420 that are close in the state space may have distinct patterns of transitions, and 421 subsets of attractors that are far apart in the state space may have similar patterns 422 of transitions. In other words, precaution may be used when treating dFC patterns 423 as brain states.

424
Finally, the cross-attractor view attaches the concept of energy costs to functional 425 connectivity patterns. Mathematically, cross-attractor coordination matrices and 426 energy gaps are different depictions of inter-attractor relations. Cross-attractor 427 coordination matrices describe in which direction the attractors fall in line with each 428 other; the energy gaps depict the spacing between attractors in a predefined direction 429 (such a direction can represent the whole brain or a specific subnetwork). Although 430 the pattern of cross-attractor coordination can remain similar as the landscape 431 changes, it is not the case for its energy cost. In other words, the potential for 432 exhibiting normal functional connectivity patterns may always be there, but the 433 energy costs modulate the difficulty for such potential to be realized.

434
In summary, the present work examines intrinsic brain dynamics in terms of an 435 underlying landscape and the repertoire of stable activity patterns it affords. Model-436 based analyses reveal that empirically observed functional connectivity patterns may 437 reflect transitions between activity patterns more than the patterns themselves. The 438 work outlines a modeling framework that emphasizes the relation between stable 439 activity patterns. It is thus suitable for examining systemic changes in the brain 440 that result in interrelated improvement or impairment in multiple cognitive and 441 affective functions, such as in development and psychiatric disorders.
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. In particular, they are sigmoidal functions of the form whose output increases with input monotonically and saturates at r max -the maximum firing rate limited by the absolute refractory period of neurons (around 2 ms in certain cell types [65, 66]). The specific shape of each transfer function is determined by three additional parameters a p , b p and d p (a p and b p determine the location and slope of the near-linear segment in the middle; d p determines the smoothness of the corners bordering the said near-linear segment). This transfer function is converted from Wong and Wang's original formulation [39, 67] (a soft rectifier function, equation S6) into a sigmoidal form, while retaining the original value of parameters a p , b p , and d p (shown in Table 1). The parameters were chosen to approximate the average response of a population of spiking pyramidal cells (p = E) and interneurons (p = I) respectively, incorporating physiologically plausible parameters [39,68]. 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 ( 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 445 region, one would expect it to have comparable influence on the local dynamics as 446 I E in the local model (equation 1).
inhibitory-to-inhibitory coupling 0.05 (nA) I E external input to excitatory population ∼ (nA) I I external input to inhibitory population 0.1 (nA) G global coupling ∼ (nA) C ij structural connectivity between brain regions ∼ σ noise amplitude ∼  of the Jacobian matrix is used to classify the fixed points and identify which ones 454 are attractors. The fixed point is a stable equilibrium if λ k is real and negative 455 for all k. The fixed point is associated with damped oscillation if Re λ k < 0 for 456 all k and Im λ k = 0 for some k. The fixed point is associated with a limit cycle if 457 Re λ k > 0 and Im λ k = 0 for some k with the additional criteria that after a small 458 perturbation from the fixed point, the time-average of the solution remains close to 459 the fixed point. The above three types of fixed points-a stable equilibrium, a stable 460 spiral (damped oscillation), a fixed point associated with a limit cycle (sustained 461 oscillation)-represent attractors in the present study. All other types of fixed points 462 are classified as unstable. For damped oscillation and limit cycles in the local model, 463 the frequency of the oscillation ( Figure S1) is defined as | Im λ k |/(2π).

464
For the local model, a 2D dynamical system, the complete characterization of 465 all fixed points is relatively easy by searching exhaustively through a grid of initial 466 guesses (as for Figure 3a-c). This approach becomes unfeasible when it comes to the 467 global model due to the high dimensionality. Thus, for the global model (Figure 3d-i), 468 we implemented a recursive search: for each value of G, (1) find zeros of equation 4-6 469 using fsolve given a set of initial guesses that includes, if any, the zeros for G − δG 470 (δG = 0.01 for the present study) in addition to a fixed set of grid points; (2) sort 471 the list of zeros obtained from (1) by the average of S (i) E 's denoted asS E ; (3) use 472 the middle points between consecutive zeros in the sorted list as initial guesses; (4) 473 continue to use middle points between past initial guesses as new initial guesses 474 recursively until at least one new zero is found or the recursion has reached a certain 475 depth; (5) append the new zero(s) to the list of zeros and repeat (2)-(5) until the 476 number of identified zeros exceeds a certain value. In the present study, we limit the 477 maximum depth in (4) to 8 and the maximum number of zeros in (5) to 200. The 478 set of zeros so obtained are the fixed points of the dynamical system. Each fixed 479 point is further classified using the respective Jacobian matrix as described above to 480 identify the attractors-a subset of the fixed points forming the attractor repertoire. 481 For each set of structural parameters (G, C, w EE , w EI ), we represent the 482 attractor repertoire as a M-by-N matrix, where M is the number of attractors, and N is the number of brain regions.

484
As parameters vary (e.g., G in Figure 3d-i), the attractors form discrete connected 485 components (e.g., stripes in Figure 3) in the product of the state space and the 486 parameter space. Attractors within the same connected component can be considered 487 qualitatively equivalent, as they morph into each other under continuous parameter 488 change. Understanding the relation between these connected components are critical 489 to characterizing phase transitions and bifurcations. It is thus meaningful to define 490 a discretized version of the attractor repertoire, E,i is positive integer representing a discrete level of activation which S  Given a discretized attractor repertoireÂ (equation 9), the cross-attractor coordi-501 nation matrix is an N-by-N matrix, whereÂ ·,j denotes the j th column ofÂ and ρ(x, y) the Spearman's correlation 503 between variables x and y (see Figure 2 for toy examples, Figure 4b for a more 504 elaborate example). Spearman's correlation is chosen to reflect the ordinal nature of 505 the variableŜ gives the level of cross-attractor coordination 506 between model brain region i and j. The use of the discretized attractor repertoire 507 A ensures that coordination matrix P is invariant within the same dynamic regime 508 and only changes during a bifurcation. Thus, matrix P connects the change of 509 brain coordination patterns to dynamical systems concepts such as bifurcation-a 510 qualitative change in the dynamic landscape of the model brain.

511
The coordination matrix P by itself does not explicitly concern how difficult or 512 energy consuming these cross-attractor movements are. To incorporate energetic 513 properties, we equip each attractor repertoire with a sequence of energy gaps. We 514 first order the rows of the attractor repertoire matrix A so that the row averages 515 descend with the row index. The row averages of the ordered repertoire matrix 516 provide a sequence of energy levels, whereS E,i = N −1 N j=1 A ij ,S E,i >S E,i+1 for any i < M , and M is the number of 518 attractors. The corresponding energy gaps are where ∆e i =S E,i−1 −S E,i is the energy gap between the (i − 1) th and i th attractor 520 in the repertoire. Physically, each energy gap ∆e i can be interpreted as the energy 521 cost associated with keeping additional x% synaptic channels open. The sequence of 522 energy gaps can be used to partition the attractor repertoire into submatrices. For 523 example, if ∆e i is the maximum energy gap, one can split A (and its discretized 524 versionÂ) into a repertoire above the energy gap and a repertoire below the energy gap where A i,· is the i th row of the ordered repertoire matrix A. Each of the sub-527 repertoires A + and A − is associated with its own cross-attractor coordination 528 matrix, say, P + and P − (equation 10). P + and P − can be considered as the 529 "energy-constrained" coordination patterns where the model brain is not allowed to 530 cross the energy gap ∆e i . The same line of analysis is applicable to a sub-network 531 of the model brain (e.g., the default mode network) by constructing a reduced 532 repertoire matrix that only contains a subset of the columns in A. These selected 533 columns map to the brain regions within the sub-network. This series of analysis 534 applied to attractor repertoire matrices and sub-matrices provides a systematic and 535 simple way to characterize brain dynamic landscape and inter-regional coordination.   The human structural connectome used in the present study is from the S1200 552 Release from the Human Connectome Project (HCP) [43]. The average connectome 553 of 11 unrelated subjects were used for the first qualitative analysis (it has been 554 shown that averaging over 5 subjects is sufficient [12,15] (Figure 1b). The original diffusion imaging (dMRI) 559 data were obtained using a customized Siemens 3T scanner at Washington University 560 in St. Louis, with a standard 32-channel head coil, with TR = 5520 (ms), TE = 561 89.5 (ms), 1.25 (mm) isotropic voxels, b=1000, 2000, 3000 (s/mm 2 ). T1 images were 562 obtained using 3D magnetization-prepared rapid gradient echo sequence (MPRAGE) 563 with TR = 2400 (ms), TE = 2.14 (ms), and 0.7 (mm) isotropic voxels. The HCP 564 minimally processed data were further processed using MRtrix3, including bias-565 field correction, multi-shell multi-tissue constrained spherical deconvolution with 566 a maximum spherical harmonic degree 8. 10 million probabilistic streamlines were 567 generated for each subject using the 2 nd -order Intergration over Fibre Orientation 568 Distributions algorithm (iFOD2) [70] and anatomically-constrained tractography 569 (ACT) [71] (FOD amplitude threshold = 0.06, step size = 0.625 mm). Each streamline 570 was assigned a weight using spherical-deconvolution informed filtering of tractograms 571 (SIFT2) [72]. Connection strengths between ROIs are summed weights of the 572 associated streamlines. Intra-ROI connections are removed. Subjects' connectivity 573 matrices are normalized according to equation 7 before and after averaging. Human functional connectivity used in the present study is estimated using the 576 resting-state fMRI (rsfMRI) data of the same subjects from the Human Connectome 577 Project [43] as aforementioned ones from the structural connectivity. rsfMRI scans 578 were acquired using EPI sequences with TR = 720 (ms), TE = 33.1 (ms), flip angle 579 = 52 • , voxel size = 2.0 (mm, isotropic), multiband factor = 8. Four runs of rsfMRI 580 scan were obtained from each subject in 2 separate days (2 runs in each day with 581 opposite phase-encoding direction: RL and LR). Each run last 14 min 33 s (1200 582 TR). . First, a reference volume and its skull-stripped 587 version were generated using a custom methodology of fMRIPrep. A deformation 588 field to correct for susceptibility distortions was estimated based on two echo-589 planar imaging (EPI) references with opposing phase-encoding directions, using 590 3dQwarp [77] (AFNI 20160207). Based on the estimated susceptibility distortion, 591 an unwarped BOLD reference was calculated for a more accurate co-registration 592 with the anatomical reference. The BOLD reference was then co-registered to the 593 T1w reference using bbregister (FreeSurfer) which implements boundary-based 594 registration [78]. Co-registration was configured with nine degrees of freedom to 595 account for distortions remaining in the BOLD reference. Head-motion parameters 596 with respect to the BOLD reference (transformation matrices, and six corresponding 597 rotation and translation parameters) are estimated before any spatiotemporal filtering 598 using mcflirt [FSL 5.0.9, 79]. The BOLD time-series were resampled to the fsaverage 599 surface space. Several confounding time-series were calculated including framewise 600 displacement (FD), DVARS and three region-wise global signals. FD and DVARS 601 were calculated for each functional run, both using their implementations in Nipype 602 [following the definitions by 80]. The three global signals were extracted within the 603 CSF, the WM, and the whole-brain masks.  For the second analysis, an alternative preprocessing pipeline that is standard for 617 individual rsfMRI preprocessing was applied. Specifically, minimally preprocessed 618 rsfMRI volumetric data from the HCP database was downloaded and further denoised 619 through application of a published spatial independent component rejection procedure 620 [81-83]. The resulting rsfMRI timeseries was then converted to ROI space by 621 averaging voxels part of the Desikan-Killiany atlas in volume space [84]. Relative 622 to the list in Table S1, there were six missing ROIs, which are the left and right 623 temporal pole, the left and right frontal pole, and the left and right bank of the 624 superior temporal sulcus. Similar to the case with the 11 subjects, there were 4 625 runs per subject (2 runs with opposite encoding for 2 separate days). To calculate 626 the functional connectivity, ROI time series were z-scored and then correlated with 627 Spearman correlation between each pair of ROIs. The resulting 4 connectivity 628 matrices were averaged to obtain one final connectivity matrix per subject. To find the optimal human-model fitting on an individual basis, simulations were 631 run with one representative set of local parameters (W EE = 2, W EI = 1), while the 632 global coupling parameter was varied from 1.7 to 3.0 with step sizes of 0.1. This 633 set of parameter configurations were chosen with considerations for the run time 634 while capturing the individual variability for the optimal value of G. Additional 635 details of the simulation, including computation of the within and cross attractor 636 coordination, and energy gaps between attractors, were identical to the analysis of 637 the 11 subject averaged data.

638
Spearman's correlation was used to quantify the model-human correlation. The 639 optimal G parameter for each individual was determined based on the maximum 640 correlation for the cross-attractor coordination. Within-attractor coordination 641 matrices were determined for individual attractors at that optimal G and used to find 642 the highest model-human correlation for within-attractor coordination. To compare 643 the model fit between within and cross-attractor coordination on an individual level, 644 a paired t-test was applied between the model-human correlation for the two across 645 the 100 individuals.        Table S1: Names and indices of brain regions. As in [1], the brain is parcellated into 66 regions as shown in Figure 1b. Here is a list of the specific region names corresponding to each region index. Region 1-33 are on the right hemisphere, and region 66-34 are the homologous regions on the the left hemisphere.

S2 Relation to the Wilson-Cowan model
τ E and τ I are time constants of the dynamics of the excitatory and inhibitory 861 population respectively. c • 's are the coupling parameters between the two population. 862 Coefficients k • and r • result from a temporal coarse-graining procedure in the initial 863 derivation (see [2] for detail). S • is a sigmoidal transfer function, rising monotonically 864 from 0 to 1 with non-negative input. P and Q are external inputs to their respective 865 populations. If we divide both sides of equation S1-S2 by the time constants, we are 866 looking at the same general form as equation 1-2.

867
The main difference is between the respective transfer functions. Wilson and Cowan [2] chose a particular form of S for mathematical analysis:

II S
(i) following the same notations as in equation 4-6, wherẽ with p ∈ {E, I} denoting the excitatory and the inhibitory population respectively. 891 The parameters a p , b p and d p were chosen such thatH p approximates the average 892 firing rate of an ensemble of leaky integrate-and-fire neurons receiving uncorrelated 893 noisy inputs.

894
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 895 potential of the membrane. The total synaptic input current I syn (t) is a random 896 process with an average µ C and standard deviation σ C . When V (t) reaches a 897 threshold V th , the neuron emits a spike after which the membrane potential returns 898 to a reset voltage V reset and stays there for a duration τ ref , i.e. the refractory 899 period.

900
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. [8, Section 1], [9] for descriptions of the Fokker-Planck approach). This eventually leads to the first-passage time equation (average time for crossing the threshold), and V ss the steady state voltage The transfer function employed by Wong and Wang [5, 10], i.e. equation S6 with 901 appropriate choice of parameters, is a good approximation of equation S8 when the 902 input level is low.

903
Thus, the first passage equation S8 provides a bridge between the transfer function 904 (equation S6) and the single-cell level model (equation S7) incorporating realistic 905 biophysical parameters (Table S2). In other words, it allows one to use empirically 906 measurable quantities at the neuronal level to directly constrain the the transfer 907 function and the entire model. This is a major difference with the Wilson-Cowan 908 model [2, 3].

909
parameter interpretation value C m membrane capacitance 0.5, 0.2 (nF) g L leak conductance 25, 20 (nS) τ m membrane time constant 20, 10 (ms) τ ref refractory period 2 (ms) V L resting membrane potential -70 (mV) V th threshold for firing -50 (mV) V reset reset potential -55 (mV) Table S2: 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 S6).
According to the first passage equation S8, 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 S6).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). H p 910 matchesH p for low levels of input but flattens out eventually at r max as one would 911 expect from equation S8.  , c). A typical phase portrait from each regime is provided in Figure S2.
The local model exhibits a rich repertoire of dynamical features, including mul-918 tistability (b, c in Figure S1, S2), damped oscillation (c, d, f), and limit cycles 919 (sustained oscillation; a, b). Mathematical analysis of the local model (Section S13 in 920 Supplementary Materials) shows that nonlinearity in the dynamics can essentially be 921 controlled by two local structural properties: the strength of excitatory-to-excitatory 922 connection w EE and the strength of the excitatory-to-inhibitory connection w EI . 923 Geometrically, the two structural properties "twist" the nullclines (dashed lines in 924 Figure S2). Specifically, stronger w EE introduces a deeper twist and fold of the 925 red nullcline (compare Figure S2a and d), whereas stronger w EI introduces a more 926 vertical twist of the blue nullcline (compare Figure S2d and e). These twists are 927 the key sources of dynamic complexity-multistability and oscillation. For example, 928 when w EE is sufficiently large (equation S30), multistability becomes possible: the 929 folded red nullcline allows for multiple intersections with the blue nullcline, and 930 potentially a greater number of attractors (compare Figure S2c to a; see analytical 931 results in Section S13: Multistability). When w EI is sufficiently large (equation S38), 932 oscillatory activity becomes possible (analytical results in Section S13: Oscillation). 933 Moreover, the combination of large w EE and w EI gives rise to sustained oscillation 934 (equation S61). The characteristic frequency of such oscillation further depends on 935 the specific values of w EE and w EI . Note that the general qualitative effects of 936 these two local structural properties are consistent with those of the Wilson-Cowan 937 model, but the specific boundaries at which transitions occur are determined by the 938 biophysical constraints inherited from the reduced Wong-Wang model (see equa-939 tions S38, S61). Analytical results (Section S13) provide detailed quantification of 940 how these boundaries are shifted by different local structural properties.

941
To maintain a sufficient twist in the red nullcline (red dashed line in Figure S2) 942 and associated multistability and oscillation, inhibitory-to-excitatory feedback w IE 943 needs to be proportional to self-excitation w EE (c.f. equation S20). In the present 944 study, we simply let w IE = w EE . This equality is a simpler alternative to the 945 Feedback Inhibition Control adopted in [6] in both numerical and mathematical 946 analyses.  Figure S1. The specific parameters defining local structural connectivity are (a) w EE = 4, w EI = 1; (b) w EE = 4, w EI = 0.8; (c) w EE = 2.3, w EI = 0.75; (d) w EE = 1.5, w EI = 1; (e) w EE = 1.5, w EI = 0.5; (f) w EE = 1.5, w EI = 0.3; (g) w EE = 1.5, w EI = 0.2. 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.

S5 The effects of nonlinearity in the local model 948
Before getting into the global model, we briefly demonstrate numerically how the 949 present unified model (equation 1-2) extends the reduced Wong-Wang model [5, 6] 950 to more complex scenarios. As expected, the dynamics of the present model match 951 that of the reduced Wong-Wang model for low levels of excitation, i.e. weak local 952 excitatory connectivity (Figure S3a-b). In a regime of stronger local excitatory 953 connectivity, as explored in [7], the two models diverge (Figure S3c-d). In the present 954 model ( Figure S3c), all trajectories are well-confined within a physiologically plausible 955 range-state variables S E and S I denote the fraction of open channels, which by 956 definition are between 0 and 1. In contrast, certain trajectories of the reduced 957 Wong-Wang model ( Figure S3d) overshoot beyond the physiologically plausible 958 range. The effect of added nonlinearity in the present model manifests through the 959 curvature of the blue nullclines, which confines the flow of oscillatory activities and 960 creates extended multistability (see e.g. Figure S2b). Thus, the present model is 961 more 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).

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

967
Given a uniform global network, temporal diversity across the whole brain can be 968 induced by the diversity of local excitatory-to-excitatory connection (w EE ), as shown 969 in Figure S4a. Brain regions with relatively weak w EE (blue) have low characteristic 970 frequencies around 10 Hz (alpha range), while brain regions with strong w EE (red) 971 have higher characteristic frequencies around 30 Hz (beta/gamma range). In other 972 words, the characteristic frequency of the oscillation increases monotonically with 973 w EE (see also Figure S8a). This is expected from the behavior of isolated brain 974 regions ( Figure S1d). In addition to the expected diversity, signs of coordination 975 between regions can be seen as the wide-spread alpha peaks ( Figure S4a). In 976 contrast, regions with a higher characteristic frequency (beta/gamma range) are 977 not as influential to other regions. That is, low-frequency oscillations, rather than 978 high-frequency ones, are responsible for global coordination.

979
The above observations concern high-frequency dynamics typically measured 980 using, e.g. electroencephalography (EEG) and magnetoencephalography (MEG). For 981 low-frequency dynamics typical for functional magnetic resonance imaging (fMRI), 982 we examine the low-frequency content (0.01-0.1 Hz) of the normalized power spectra 983 of BOLD activities, derived from the same simulated neural dynamics (see Section S7 984 in Supplementary Materials for details). The result is shown in Figure S4b: there is 985 no significant dependency of low-frequency power on w EE (Spearman correlation 986 ρ = −0.029, p = 0.81). In short, we find differential effects of local structural 987 diversity on neural dynamics at the time scales typical for different neural imaging 988 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 S7 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 [11,12] 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 990 themselves to be structurally different. As shown in Figure S4c-d, locally identical 991 brain regions can behave very differently due to the topology of the large-scale 992 network (human connectome as in Section 2.2). The influence of large-scale structural 993 connectivity on temporal diversity is manifested in both the high-frequency neural 994 dynamics ( Figure S4c; Figure S9a) and the low-frequency power of the BOLD signals 995 ( Figure S4d; Figure S9b). Specifically, the low-frequency power is inversely related 996 to the degree of each node (brain region) in the large-scale network (Spearman 997 correlation ρ = −0.584, p < 10 −6 ).

998
In Section S12, we demonstrate that the above effects are robust over 200 sim-999 ulated trials of the same parameter settings. Overall, both local ( Figure S4a,b; 1000 Figure S8) and large-scale structural connectivity ( Figure S4c,d; Figure S9) con-1001 tribute to the diversification of local dynamics. The contribution of local structural 1002 differences is stonger in a higher-frequency range ( Figure S4d; Figure S8a), while 1003 the contribution of global structural connectivity is stronger in a very-low frequency 1004 range ( Figure S4d; Figure S9b). Modeling real neural dynamics requires considering 1005 both sides of the spectrum.

S7 Computation of BOLD signal and low-frequency power 1007
In the present study, we are interested in not only the high-frequency activity 1008 measurable by, for example, EEG recordings but also low-frequency fluctuations 1009 that are often a subject of investigation in fMRI studies. Therefore, we simulated 1010 the BOLD activities induced by the underlying neural dynamics and examine their 1011 low-frequency properties.

1012
BOLD (Blood-oxygen-level-dependent) activities are computed using the Balloon-Windkessel model [13][14][15][16]. The hemodynamic response of the i th brain area takes the formṡ where the interpretation and value of the parameters are given in Table S3. The 1013 initial condition is The power spectrum for each simulated BOLD time series is computed using Welch's method [17], after being subsampled at 720ms intervals (matching the TR of resting state fMRI used in the Human Connectome Project [11]). The full power spectrum P (ω) was first normalized such that where ω N is the Nyquist frequency (approximately 0.7 Hz for the chosen sampling interval). The low-frequency power is defined as S8 Discretization of regional states 1018 Although the dynamic landscape of the global model can be quite complex (Figure 3g-1019 i), each region in the globally connected network still falls into discrete states 1020 ( Figure S5) very much like in the local model (disconnected stripes in Figure 3a-c)-1021 it is the combination of regional states that produces a great variety of attractors at 1022 the global level. Discretized regional states (number on black disks in Figure S5) 1023 thus give rise to discretized attractors in the global model. Rank correlation 1024 (Spearman) between regional states across these discretized attractors are used to 1025 quantify cross-attractor coordination in the model brain ( Figure 4b in the main text). 1026 Using discretized attractors, we quantify how much two regions move up and down 1027 together across attractors without considering the distance between the attractors. 1028 The distance between attractors are considered separately as the energy cost that 1029 constrains such transitions ( Figure 5 in the main text).  The average functional connectivity is highly similar between Day 1 (Figure 4a in 1032 the main text) and Day 2 ( Figure S6a). Linear regression analysis indicates that 1033 the Day-1 matrix is highly predictive of the Day-2 matrix (red x in Figure S6b; 1034 β 1 = 0.98, t(2143)=203, p<0.001, R 2 = 0.951; only elements below the diagonal are 1035 compared due to the symmetry of the matrix). Moreover, the functional connectivity 1036 matrix ( Figure 4a in the main text) obtained using Spearman correlation is highly 1037 predictive of the corresponding Pearson correlation coefficients ( Figure S6c; β 1 = 1.03, 1038 t(2143)=897, p<0.001, R 2 = 0.997), which itself is highly consistent across two days 1039 (blue o in S6b; β 1 = 0.97, t(2143)=211, p<0.001, R 2 = 0.954).
1040 Figure S6: 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 4a (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 4a (main text). See text for statistical information. S10 Mean comparisons of human-model similarity 1041 Figure S7: Average effects of local structural connectivity on model-human similarity. Cross-attractor coordination in the model well captures human functional connectivity for both intra-(red) and inter-hemispheric interaction (blue). The model-human correlation (Spearman's ρ) is slightly worse when local excitatory connection is too strong (w EE = 2.8, w EI = 1). (* p<0.05, *** p<0.001, with Tukey HSD) S11 Results of permutation tests and bootstrapping   Figure 6) are based on 30000 random permutations. The significance levels are Bonferroni-corrected. A comparison reaches statistical significance if the mean is outside of the confidence interval.  Table S5: Permutation tests for the maximum energy gaps by three different levels of local excitatory connectivity (w EE , w EI ). The mean comparisons between different levels of w EE , w EI (bars in Figure 5b) are based on 30000 random permutations. The significance levels are Bonferroni-corrected. A comparison reaches statistical significance if the mean is outside of the confidence interval.   Figure 5c) are based on 30000 random permutations. The significance levels are Bonferroni-corrected. A comparison reaches statistical significance if the mean is outside of the confidence interval. S12 Dependency of spectral properties on local and global 1043 structural connectivity 1044 In Figure S4, we illustrate with two simulated trials how high-frequency and low-1045 frequency dynamics depend on local excitatory-to-excitatory connectivity w EE and 1046 the topology of the global network. To show that these effects are not incidental, we 1047 simulated 200 trials for each of the conditions: (1) the global network is uniform but 1048 local connectivity w EE is diverse (as in Figure S4a,b), and (2) local connectivity 1049 w EE is identical but the global network follows the human connectome (as in 1050 Figure S4c,d). We characterize the high-frequency content of a spectrum as its peak 1051 frequency, i.e. the frequency at which the spectral power is the highest (e.g. peaks 1052 in Figure S4a,c); the low-frequency content as the integral of the power between 1053 0.01 and 0.1 Hz (Section S7). The dependency of these features on local (w EE ) and 1054 global structural properties (node degree) is quantified using Spearman correlation. 1055 The distributions of the correlation coefficients (ρ) and corresponding p-values are 1056 shown in Figure S8 and Figure S9 for condition 1 and 2 respectively. Figure S8 1057 shows that local structural connectivity w EE strongly affects the peak frequency of 1058 the brain region (a,c) but not so much the low-frequency power (b,d). The stronger 1059 the local connectivity, the higher the peak frequency. Figure S9 shows that the node 1060 degree of the global network has a strong and negative effect on the low-frequency 1061 power, and a weak and positive effect on the peak frequency. Figure S4 illustrates 1062 such dependencies using typical trials (median correlation coefficients) from the 1063 distributions ( Figure S8-S9).
1064 Figure S8: 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 S4a,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 S4a, 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 S4b shows this lack of dependency in an example trial that corresponds to the median of the distribution (b). Figure S9: 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 S4c,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 S4c, 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 S4d illustrates this dependency with an example trial that corresponds to the median of the distribution (b).
nullclines are confined within the interval Within this interval S I = f (S E ) (equation S20; red nullcline), overall, goes down 1066 from +∞ to −∞, while S E = g(S I ) (equation S21; blue nullcline) goes up from −∞ 1067 to +∞. This results from the dominant effect of H −1 p for a very large or very small 1068 input.

1069
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 S20 becomes In this simplified case, increasing self-excitation w EE reduces the influence of H −1 E 1070 such that the slope of middle branch approaches 1.

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

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

1094
This result highlights the importance of self-excitation w EE in equation S20 1095 and equation S23-multistability can only occur when w EE is sufficiently large. 1096 Correspondingly in the numerical result ( Figure S1), the region of multistability 1097 appears only for larger w EE 's.

1098
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 S8). 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 Instead of finding the minimum of equation S25, we aim to find a representative 1099 point S * E such that equation S25 is relatively small.

1100
One option is to use the minimum of the numerator. The minimum of the numerator . By design, H E reaches its maximum 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 1101 S * E is very close to one, which further results in a small denominator in equation S25 1102 and a large slope for H −1 E . Thus, the inflection point of H E (v) is not a very good 1103 choice.

1104
To avoid the small denominator problem for equation S25, 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 S26, we need and Now for the nullcline S I = f (S E ) to have a positive slope at S * E , one simply needs 1105 w EE > h E .
Here h E is approximately 0.2 based on the present parameter choices, inherited from 1106 Wong and Wang's initial derivation [5]. This result is confirmed numerically by the 1107 bifurcation diagrams (Figure 3a-c vs. Figure S10a) of the local model-multistability 1108 exists for some level of input I E when w EE > 0.2.

1109
Oscillation. Now we look for the conditions for oscillation to emerge. Here we 1110 are mainly concerned with the oscillation occurring on the ascending segment of 1111 S I = f (S E ) (red nullcline). Following a similar argument as Wilson and Cowan 1112 [2], one notice that for the flow around a fixed point-an intersection between the 1113 nullclines-to have consistent rotation, the nullcline g −1 (S E ) (blue) must have a 1114 greater slope than f (S E ) (red nullcline). Qualitatively, one would expect oscillation 1115 to be induced by increasing w EI , which twists g(S I ) (blue nullcline) counterclockwise. 1116 This expectation is confirmed by the numerical results in Figure S1a-d: oscillation 1117 emerges for sufficiently large w EI for fixed points on the ascending branch of 1118 S I = f (S E ) ( Figure S2a-d). 1119 Quantitatively, we consider the derivative of the two nullclines at a respective representative point. First, we extend the results in equation S28-S29 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 1120 w EI > w II + h I .
Given the present parameter choices, we need w EI > 0.45 to induce oscillation for 1121 some level of input I E and I I . This is in line with the numerical results in Figure S1. 1122 For h E > 0, as assumed here, lowering w EE also lowers the threshold for oscillation. 1123 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 S20-S21) intersect. Conveniently, we let them intersect at their respective representative points (S * E , f (S * E )) and (g(S * I ), S * I ) (see equation S28 and equation S31), 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 1124 parameters used in the present study, we have α E ≈ 14, α I ≈ 111, β E ≈ 71, and 1125 β I ≈ 276. 1126 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 re-1127 quires w EE > 2. Correspondingly in the numerical results ( Figure S1), equation S61 1128 provides an estimate of the lower bound of the Hopf bifurcation (gray dashed line). 1129 Note that stronger inhibitory-to-inhibitory connection w II increases the minimal 1130 w EE required to induce sustained oscillation. Overall, these analyses show that 1131 sustained oscillation requires both strong self-excitation and a sufficiently active 1132 inhibitory population.

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

1156
Proof. Assume there are at least two distinct fixed points of the system: S * and S * * , 1157 where S = ( S E , S I ) and S p = (S (1) p , · · · S (i) p , · · · S (N ) Figure S10: 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 3a, 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 3g, h, i). See text for further discussion.
What we have not addressed in the above analyses is to what extent the global 1182 system is multistable-what is the number of stable states, or the size of the memory 1183 capacity-and what are the contributions from local self-excitation and global 1184 network connectivity. An analytical approach to this problem is difficult; thus, it 1185 is mainly addressed numerically (c.f. Figure S10 and Figure 3). Nevertheless, we 1186 provide an intuitive argument below as to how local and global connectivity affects 1187 the relevant geometrical properties of the dynamical system.

1188
Local origin of geometrical complexity. At an intuitive level, the number of 1189 intersections between these hypersurfaces (equation S65-S66) is likely to increase with 1190 the number of folds of each surface. In the present case, the folding of hypersurfaces 1191 entails the temporary reversal of the sign of its partial derivative along a certain 1192 direction. Observe equation S65 and see that global coupling cannot create any 1193 folding of the surfaces. Thus, the geometrical complexity of the nullclines purely 1194 depends on the local properties of each node, in particular, the folding effect of 1195 self-excitation w way dependent on the structure connectivity C ij . This may remove or introduce new 1199 intersections between the surfaces without changing the geometrical complexity of 1200 these surfaces. Thus, global coupling allows system-level multistability to be created 1201 synergistically, given appropriate structural connectivity C ij .

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