Normalized Connectomes Show Increased Synchronizability with Age Through Their Second Largest Eigenvalue

The synchronization of different brain regions is widely observed under both normal and pathological conditions such as epilepsy. However, the relationship between the dynamics of these brain regions, the connectivity between them, and the ability to synchronize remains an open question. We investigated the problem of inter-region synchronization in networks of Wilson-Cowan/Neural field equations with homeostatic plasticity, each of which acts as a model for an isolated brain region. We considered arbitrary connection profiles with only one constraint: the rows of the connection matrices are all identically normalized. We found that these systems often synchronize to the solution obtained from a single, self-coupled neural region. We analyze the stability of this solution through a straightforward modification of the Master Stability Function (MSF) approach and found that synchronized solutions lose stability for connectivity matrices when the second largest positive eigenvalue is sufficiently large, for values of the global coupling parameter that are not too large. This result was numerically confirmed for ring systems and lattices and was also robust to small amounts of heterogeneity in the homeostatic set points in each node. Finally, we tested this result on connectomes obtained from 196 subjects over a broad age range (4-85 years) from the Human Connectome Project. We found that the second largest eigenvalue tended to decrease with age, indicating an increase in synchronizability that may be related to the increased prevalence of epilepsy with old age.


Introduction
Neuroscience has gained the ability to non-invasively map the structural and functional connectivity of the brain in a diverse group of subjects with the emergence of large scale connectomics initiatives [Van Essen et al., 2013, Behrens andSporns, 2012]. This capability has created a cornucopia of neuroimaging data in the form of connection matrices or connectomes. These connectomes determine how brain regions are coupled either structurally via white matter tracts, or functionally, through correlations. However, a natural problem emerges with this glut of data: How do we use it?
The initial approach in analyzing these connectomes originates from graph theory with the nodes of the network corresponding to the vertices of a mathematical graph and the elements of the connection matrix determining the edges of the graph [Bondy et al., 1976]. The connectomes are then analyzed for graph-based statistics such as the connectivity, mean in/out degrees, average shortest path length, average motif strengths and frequencies, or community structures, to name a few [Van den Heuvel et al., 2016]. Typically, these statistics are average properties of the connections between individual nodes or small populations of nodes over the entire graph, although some may be more global in nature. These graph-based statistics have proven useful for summarizing and elucidating the differences between connectomes of healthy individuals, and those with pathological conditions [Van den Heuvel et al., 2016].
However, the dynamics of a large network are not solely determined by graph based statistics. Indeed, dynamical systems theory tells us that the eigenvalues of the connection matrix are often the critical determining factor of the stability of any particular solution [Schaub et al., 2015]. This statement was clarified further in [Pecora and Carroll, 1998] through the derivation of the so-called Master Stability Function (MSF). The MSF determines under what conditions a dynamical system synchronizes and only depends on the eigenvalues of a connection matrix, rather than any explicit graph-based measure.
Given the power of the MSF approach, computational neuroscience has recently started analyzing connectome eigenvalues rather than graph-based statistics [Tang et al., 2017]. This is typically done by taking the raw structural connectome and computing the Laplacian matrix. This is a critical step required to apply MSF analysis. Owing to its construction, the Laplacian effectively assumes that the coupling between nodes is diffusive with strong recurrent coupling within a node and weak coupling between nodes [Tang et al., 2017, Nishikawa andMotter, 2010]. Unfortunately, the evidence for diffusive coupling is scant [Papo and Buldú, 2019].
Here, we consider an alternate assumption on the Diffusion Tensor Imaging (DTI) derived structural connectome that may have some empirical support. Rather than assuming diffusive connectivity, we instead assume that each neural region has the same level of incoming activity. Biologically, this could be achieved by the growth or pruning of dendritic spines at synapses, which would adjust the synaptic weights [Segal andAndersen, 2000, Bonhoeffer andYuste, 2002]. In our theoretical work, we achieve this through a normalization process, where the incoming weights for each node are divided by the total weight sum for that node. Like the diffusive coupling assumption, this normalization leads to an analytically tractable criterion for synchronization through a MSF approach. We test our approach on a network of Wilson-Cowan nodes with intranode inhibitory homeostatic plasticity [Hellyer et al., 2016, Nicola et al., 2018, Vogels et al., 2011] and various connectivity matrices. We find that the magnitude of the second largest eigenvalue of the normalized connectivity matrix is the critical determining factor for the stability of the synchronized state, while the largest eigenvalue is identical for all connectomes due to normalization. We then apply this approach to DTI derived structural connectome data from the human connectome project for n = 196 subjects covering a large age group (4-85 years). We find that the second largest eigenvalue (SLE) of these normalized connectomes decreases with age. Extrapolating from our analysis, this implies increasing synchronizability with age. This partially mirrors the epidemiology of epilepsy where prevalence increases with old age [Beghi andGiussani, 2018, Leppik et al., 2006], indicating a potential diagnostic use of the eigenvalues of the DTI derived structural connectomes.

Results
In order to apply the MSF approach, an explicit model of the dynamics of the nodes and coupling between nodes is necessary. Here, we consider a modification of the Wilson-Cowan system [Wilson and Cowan, 1972] first developed in [Hellyer et al., 2016, Nicola et al., 2018. This system is described by three dynamical variables: where E k is the activity of the excitatory population of neurons within the kth node, I k is the activity of the inhibitory population in the kth node ( Figure 1A), and W EI k is the homeostatically adjusted inhibitory weight. The population activities are variables that are confined to the interval (0, 1). This loosely corresponds to the proportion of neurons active in node k, and are thus "neural-field" or "mean-field" approximations to a large population of neurons in node k commonly. These are commonly used in computational neuroscience to investigate diverse network level phenomenon [Kilpatrick and Ermentrout, 2013, Ermentrout, 1998, Folias and Ermentrout, 2012, Coombes, 2010, Breakspear, 2017, Bressloff, 2019, Pinto et al., 1996, Park and Ermentrout, 2018. The inhibitory homeostatic weight W EI drives the excitatory population activity E k towards the homeostatic set-point p. Further, we assume that inhibition is always a local interaction within a node while excitation is global and dependent on the coupling matrix W EE .
All synaptic weights satisfy the requirement that where the excitatory/inhibitory character of the weight is carried by the signs in (1)- (3) The sigmoidal transfer function φ(x) satisfies the following properties: In numerical simulations, we will primarily consider the logistic transfer function To prevent individual nodes from saturating due to excessively high inputs, or becoming quiescent due to insufficient input, we assume that the input weights of each neural region are normalized ( Figure 1B). Mathematically, this assumption becomes the following: The two conditions (4)-(5) immediately imply that the weight matrix W EE has identical L 1 normalized rowsums. Recall that the L 1 norm of a vector x ∈ R n is defined as x 1 = n i=1 |x i |. Further, equation (4) also implies that the largest eigenvalue of W EE is W E by the Perron-Frobenius theorem.
Before we proceed further we note that it is convenient to write the coupling matrix W EE as follows:

Local Analysis of Equilibria
Before we consider the synchrony of solutions, we must first investigate when non-equilibrium solutions to (1)-(3) actually exist. For any connectome, the stable equilibrium point determined by the homeostatic set point destabilizes for a sufficiently large W E . This is independent of the connectome and only depends on the dynamics of the nodes. This results in a Hopf bifurcation curve W E = g Hopf (W IE ) where g Hopf (x) can be analytically determined (see Appendix, [Nicola et al., 2018]). This bifurcation curve has no dependence on any L 1 normalized connectome whatsoever.
This Hopf bifurcation coincides with the bifurcation curve of the single recurrently coupled node. In fact, all invariant sets (equilibria, limit cycles, etc.) for the single recurrently coupled node exist for any L 1 normalized connectome. For the single recurrently coupled node, these invariant sets are attractor states which starting with low amplitude oscillations, chaos, mixed-mode chaos, and mixed-mode oscillations ( Figure 1C-E) for progressively large W E [Nicola et al., 2018]. As we will see from the Master Stability Function Analysis below, networks with L 1 normalized connectomes can synchronize to these states depending on the eigenvalue spectrum of the L 1 normalized weight matrix, and in particular, the magnitude of the Second Largest Eigenvalue (SLE) of the normalized connectome.

Master Stability Function Analysis
To investigate this system more generally than in [Nicola et al., 2018] and determine the stability of synchronized solutions, we conducted large scale simulations of a variety of connectome types ( Figure 2). First, as in [Nicola et al., 2018], we found that networks of homeostatically coupled Wilson-Cowan nodes readily synchronized to attractor states of the single recurrently coupled node system with a self coupling strength of W E ( We refer to the non-equilibrium point attractor solutions to (7)-(9) as (E s (t), I s (t), W EI s (t)). This result was consistent across a variety of common coupling types such as small ring networks, lattices, small world networks. However, some connectomes could display non-synchronized solutions for a given value of W E such as larger rings (N > 9), weak coupling networks, and larger lattices ( Figure 2E). Thus, some connectomes have the capability to destabilize the synchronized solution arising from the single recurrently couple node.
To investigate the loss of stability as a function of the underlying weight matrices, we slightly modified the traditional MSF approach (See Appendix) ( Figure 3A). By linearizing around the synchronized solution (E s (t), I s (t), W EI s (t)) and assuming that W EE is diagonalizable, one can derive the following diagonalized variational equation for perturbations to the synchronized solution: wherer k is an eigenvalue of W EE . Note that for all subsequent simulations and plots, we will consider r k =r k W E , which corresponds to an eigenvalue of the L EE , which has a maximal eigenvalue of 1, rather than the eigenvalues of W EE While equation (10)- (12), may not seem particularly useful at first glance, it can compute the stability of the synchronized state for any connectome with a simple numerical algorithm [Pecora and Carroll, 1998] by decoupling the computation of 3N eigenvalue problems into N separate 3 dimensional systems. One exploits this by treating the equations (10)-(12) as a general form and using them (in conjunction with (7)-(9)) to compute the Lyapunov exponents over a mesh in the (Re(r), Imag(r)) space ( Figure 3B). Then, the eigenvalues can be computed for any connectome W EE with the Lyapunov exponents of the synchronized solution immediately determined simply by a look-up operation on the pre-computed mesh ( Figure 3B).
We applied this procedure and found that for values of W E near the Hopf bifurcation of the single recurrently coupled node, the maximum Lyapunov exponent is positive for sufficiently large eigenvalue magnitudes, |r|. However, for a suitably large W E , the synchronized solution is stable for all |r| as all Lyapunov Exponents are non-positive for |r| < 1.
To test this result, we considered two diagonalizable connectomes, the ring and lattice solutions, and determined when the synchronized state loses instability as a function of the network size, N . Our MSF analysis shows that the eigenvalues of the ring cross into the positive Lyapunov exponent region for N = 9 ( Figure 3C). Indeed, we find that simulations of rings for N = 8, 9, 10 shows that the synchronized solution is stable for N = 8, and unstable for larger N . However, the attractor state for N > 8 node ring systems has periods of synchrony interspersed with bursts of asynchrony ( Figure 3D,E). For lattice-based connectomes, we find that for an √ N × √ N lattice, the MSF analysis predicts the instability of the synchronized solution initiating at √ N = 16 ( Figure 3G), which was also confirmed by our numerical simulations ( Figure 3H). Thus, for the Wilson-Cowan system considered here, the stability of the synchronized solution is largely determined by the SLE of the connectome. The larger the magnitude of the SLE, the greater the potential for desynchronization, while the smaller the magnitude of the SLE, the greater the stability of the synchronous solution.

Network Simulations with Heterogeneous Nodes
While MSF analysis can determine when synchronized solutions lose stability, it is limited by the assumption that all nodes are homogeneous [Papo and Buldú, 2019]. However, recent analysis of alternate diffusively coupled neural systems has indicated that so long as the heterogeneity is suitably small, MSF analysis is still a reasonable predictor of synchronization [Pereira et al., 2014].
To investigate if the stability of the synchronized solution was robust to heterogeneity in the nodes, we simulated the ring networks with varying levels of heterogeneity ( Figure 4). In particular, we considered the case where the homeostatic set point p for the different nodes was chosen from a uniform distribution over small intervals around the default value of p = 0.2: The plastic weight W EI now tries to drive the excitatory nodes to their distinct set points: E k → p k ( Figure  4A).
To quantify the magnitude of synchrony, we computed the Kuramoto order parameter (Materials and Methods) where φ k is the phase of the kth node, as measured by E k (t) (Materials and Methods). Fully synchronized states correspond to |R(t)| = 1 while states with a uniform (or uniformly clustered) phase distribution over the unit circle correspond to |R(t)| = 0 ( Figure 4B). For small amounts of heterogeneity ( = 10 −2 ), the stability of the synchronized solution is largely similar to the homogeneous case ( Figure 4C), with the destabilization of synchrony at a ring size of N = 9. However, larger amounts of heterogeneity between nodes ( = 10 −1 ) can completely destroy the synchronized attractor state, including restoring the stability of the equilibrium point ( Figure 4D). Thus, the numerical simulations indicate that the synchronized solution in a mildly heterogeneous system may retain the stability characteristics of the homogeneous system, as in the network considered by [Pereira et al., 2014].

The Second Largest Eigenvalue of Normalized Structural Connectomes
The analysis of the model above indicates that the Second Largest Eigenvalue (SLE) of an L 1 normalized connectome determines the stability of the synchronized solution. Larger values lead to instability while smaller values indicate greater stability of synchronized solutions. This is markedly different from stability criteria invoked by the Laplacian, where the spread of the eigenvalues is often used as the determinant of stability [Tang et al., 2017, Nishikawa andMotter, 2010] (Materials and Methods). Indeed, the eigenvalues of the Laplacian and an L 1 normalized matrix need not coincide (Appendix), thus they generate different metrics for synchronizability that need not be related.
Given the lack of correspondence between these synchronization criteria, we investigated whether or not the SLE differs across real connectomes ( Figure 5). Thus, we utilized publicly available data sets from the UCLA multimodal connectivity database (NKI RS sample) consisting of DTI connectomes from 196 subjects ranging from 4-85 years of age [Brown et al., 2012]. These connectomes consist of 188 nodes corresponding to different brain regions ( Figure 5A). For the purposes of comparison, we considered both the L 1 normalization and the Laplacian transformation of the raw-connectomes ( Figure 5B).
After normalization, we found that the SLE of the connectomes exhibited a linear decrease with age ( Figure  5C-D) as measured by a Pearson-Correlation coefficient (ρ = −0.2812, p = 6.54·10 −5 ). This decreasing strength in the SLE indicates a broad increase in the synchronizability associated with old age. We found that this was also present in the synchronizability metric of the Laplacian, albeit with a weaker correlation (ρ = 0.21, p = 2.6 · 10 −3 ).
Our analysis has shown that the SLE is a potential predictor for synchronizability. However, it is inherently global and does not implicate any particular node in the network as contributing towards synchronizability. To investigate the impacts of individual nodes further, we applied a node-deletion protocol. Nodes were individually deleted for each of the n = 196, then the resulting weights renormalized and the SLE recomputed ( Figure 5F-I). We found that across the subjects, many nodes consistently increased/decreased the SLE for the resultant matrix ( Figure 5G-H), with the most impactful nodes increasing or decreasing the SLE by no more than 1 − 2%. We found that overall, the deletions were slightly more likely to result in SLE increases (56.46% of deletions), rather than decreases (43.54% of deletions), which was statistically significant (Wilcoxon Sign-Rank test, n = 36847, p 10 − 4).

Discussion
With the glut of neuroimaging data, a natural question emerges as to how to best make use of structural connectomes. We considered a dynamical systems model in the form of a network of Wilson-Cowan nodes with intranode inhibitory homeostatic plasticity, with the coupling between nodes determined by arbitrary L 1 normalized connectomes. We considered a variety of connection matrix types, both synthetic (rings, weak coupling, lattices, etc.) and experimentally obtained from Diffusion Tensor Imaging data. Through modification of the Master Stability Function (MSF) approach, we found that networks with a sufficiently large second largest eigenvalue (SLE) could destabilize the synchronized state. This result persisted when we considered small amounts of heterogeneity in the homeostatic set points of the individual nodes, but perished for a sufficiently large distribution of heterogeneity. Finally, we tested the SLE of real-connectomes from the NKI Rockland data set and found a decrease in the SLE of the normalized-connectomes with age, indicating an increase in the network synchronizability.
Our results indicate a general increase in synchronizability with age. This, at face value, seems to contradict the decrease in synchronizability recently reported by [Tang et al., 2017]. However, there is one critical difference between the data sets we consider: the age of the subjects. In [Tang et al., 2017], the authors considered subjects in the 8-22 age group and found a general decrease in synchronizability with age in 882 developing subjects. Our result, on the other hand, consists of 196 subjects over a much broader age range (7-85), and we find an increase, rather than a decrease in synchronizability with old age. One possible parsimonious explanation is that both results are valid. In developing children, there is a general decrease in synchronizability which reaches a plateau. However, as adults enter into old age, synchronizability begins to increase again due to white matter degeneration. Thus, these combined results would indicate that for a sufficiently large sample covering a broad age range, synchronizability is a unimodal function with a minimum at adulthood. This hypothesis is supported by the epidemiology of epilepsy. The prevalence of epilepsy with age is a u-shaped curve with increased synchronizability in children and the elderly [Beghi andGiussani, 2018, Leppik et al., 2006].
Recently, the applicability of MSF analysis to brain synchronizability has been criticized [Papo and Buldú, 2019]. Two of the major criticisms levied at MSF analysis in [Papo and Buldú, 2019] are 1) The functional form of the MSF is dependent on the dynamics of the system, thus considering the eigenvalues alone can lead to erroneous results, 2) the homogeneity and diffusive coupling assumption in the MSF approach make it difficult to apply. All of the criticisms raised by [Papo and Buldú, 2019] are indeed valid. Here, we find that by considering a Wilson-Cowan system with homeostatically regulated inhibitory plasticity yielded one of two behaviours: Stability of the synchronized state regardless of the eigenvalues of the normalized matrix, or instability of the synchronized solution for weight matrices with a sufficiently large SLE. Unlike prior work ( [Tang et al., 2017]), we considered the normalized coupling case, yielding alternate MSF profiles and biologically plausible modelling for the neural regions. Further, we did find that the MSF analysis was somewhat permissive of heterogeneity in the nodes. While the criticisms in [Papo and Buldú, 2019] are valid, the analyses conducted here and in [Tang et al., 2017] do yield diagnostic measures that may be worth examining in real subjects for comparison across subjects.    The activity for the excitatory population E(t) for a single recurrently coupled node for W E = 2.115, in the chaotic dynamics regime considered by [Nicola et al., 2018]. (B) An Erdős-Rényi coupled network with N = 100 nodes, with an L 1 normalized connectome with an identical W E as (A). The nodes all synchronize onto a chaotic attractor. (C) The steady state attractors for the single-node system in (A) and the 100 node Erdős-Rényi system in (B). (D) The return maps for the excitatory populations of both networks, where the (n − 1)th peak of E(t) is plotted against the nth peak. As the return maps broadly overlap, the two attractors are identical. (E) Steady state attractors for 1) a single self coupled node at W E = 2.05, W E = 2.115, and W E = 2.25, two nodes with reciprocal coupling, a ring consisting of N = 7 nodes, A Watts-Strogatz smallworld network with N = 200 nodes, an Erdős-Rényi Network with N = 100 nodes, a 10 × 10 lattice network, a 20 × 20 lattice network, a ring with N = 10 nodes, a system with N = 50 nodes that have strong self-coupling, but random weak-coupling to other nodes. Some systems display a robust synchronization phenomenon across states, while other systems desynchronize. To derive the master stability function for the synchronized solution (SS) to any connectome, the original dynamical system is first linearized and diagonalized into a general form involving the independent perturbations to the dynamics , and the eigenvalues of the matrix W EE . The master stability function can be numerically approximated for any weight matrix by resolving the Lyapunov exponents for the general block over a mesh in r. The largest positive eigenvalue yields the master stability function. (B) The numerically derived master-stability functions for the WC system (1)

Acknowledgements
WN is funded by an NSERC Discovery Grant, and a Hotchkiss Brain Institute start-up fund. SAC is funded by an NSERC Discovery Grant. We would like to thank Joern Davidsen for many fruitful conversations about this work.

Methods
All simulations were conducted in MATLAB2019a, with the ODE45 numerical integration sub-function used for all simulations. The tolerance parameter values (RelTol and AbsTol) were increased to 10 −14 to ensure accurate integration from the default 10 −6 values.

Network Sub-type Parameters
For all networks considered, the parameters used were as in Table 1.
100 (Small), 400 (Large) N (Weak Coupling) 50 Table 1: Table of parameters for simulations, unless otherwise specified by a figure caption.

Ring Networks
Rings were generated by having each node connected to the next with a weight of W E :

Small World Networks
Small world networks are generated by creating a bidirectional ring network initially where neuron is connected to its k previous and next nearest neighbours with the weight W E 2k . Then, a parameter 0 ≤ β ≤ 1 is used to set the probability of random rewirings where node (i, j) is randomly permuted with (i, k) for each i = 1, 2, . . . N , and i < j < i + k/2. The values β = 0.7 and k = 20 were taken in all simulations.

Erdős Rényi Networks
Erdős Rényi networks were generated by first creating a random weight matrix Q EE ij = q ij where q ij was a uniform random variable on [0, 1], then normalizing rows in Q EE with W EE .

Lattice Networks
Lattices are generated by setting the nodes on a √ N × √ N square lattice. Node i, j receives input weights of strength W E 4 from node (i + 1, j), (i − 1, j), (i, j + 1) and (i, j − 1). Finally, periodic boundary conditions are applied at the edges i = N, 1 and j = N, 1.

Weak Coupling Networks
Weak coupling networks are generated by first creating a random matrix: where q ij is a uniformly distributed random number, and δ ij is the Kronecker delta. Then, W EE is given by normalizing Q EE with W E : The parameter λ controls how weak the non-self-coupling components of the weight matrix were and was set to λ = 10 −3 .

Numerically Computing Lyapunov Exponents
The Lyapunov exponents were numerically computed by the algorithm described in [Wolf et al., 1985] (see Appendix).

Numerically Computing the Kuramoto Order Parameter
The Kuramoto order parameter is computed directly by first using the findpeaks function in MATLAB2019a to detect the peaks of the chaotic or periodic solutions. These peaks are taken to have a phase of 2π at time t * , the location of the peak (φ k (t * ) = 2π) with the phase value reset to 0 in the next time step t * + ∆t. The phase is linearly interpolated for all other time points between [0, 2π]. While this is a crude approximation to the actual phase, we found that other numerical methods (e.g. the Hilbert Transform) were yielding poor estimates to the phase due to the mixed mode nature of the periodic or chaotic solutions.

Synchronization Metric of the Laplacian
Following [Tang et al., 2017, Nishikawa andMotter, 2010], we also computed the synchronization metric (σ −2 ) of the Laplacian as: where L is the Laplacian of DTI structural connectivity matrix and λ i are the eigenvalues of the L. Note that as the row-sum of the Laplacian is identically 0, only the non-zero eigenvalues (hence the N − 1 terms in the sums) are used to compute this metric Appendix A: Master Stability Function Derivation 5.5 Appendix A1: Laplacian/Diffusive Connectivity The Master Stability Function (MSF) approach introduced originally by [Pecora and Carroll, 1998] allows one to analyze the stability of synchronized solution x i (t) = x s (t), i = 1, 2, . . . N , of the coupled dynamical network: by knowing little more than the eigenvalues of of the matrix A. To apply the MSF formulation, one requires that N j=1 A ij = 0, ∀i, and thatẋ These two criteria force the synchronized solution to be an invariant set of the coupled network. The stability analysis is accomplished by first linearizing the subsequent dynamics around the invariant set x s (t) : where DF and DG denote the Jacobians of F and G. This non-autonomous dynamical system can resolve the Lyapunov exponents of the synchronized solution x s (t) for any A. However, for diagonalizable matrices A, i.e., matrices such that P −1 AP is diagonal, one can consider the substitution: which yields the following: which is a decoupled block diagonal system. This reduces the problem of determining N p Lyapunov exponents in (18) to solving separate N separate Lyapunov exponents in p dimensional dynamical systems in (20). The final step in MSF analysis is to note that (20) is a generic system of the eigenvalue r i which can be considered as a free parameter, r, thus yielding the master stability function λ max (r) as the maximum Lyapunov exponent of (20) for a generic r. The requirement that N j=1 A ij = 0 is seldom satisfied by real-world connectomes, and thus users of MSFs often assume so-called "diffusive coupling": where L ij = δ ij · N j=1 A ij − A ij . This assumption of the structural form of the underlying network model allows one to use an arbitrary connectome A ij with MSF analysis, albeit with some rigid assumptions on the nature of the coupling. Thus, in the diffusive coupling case, one considers an arbitrary connectome and analyzes the eigenvalues of the Laplacian, L, to determine the stability of the synchronized state.

Appendix A2: Local Stability Analysis of Equilibria
Somewhat surprisingly, the local analysis of equilibria is largely similar for arbitrary N , and L 1 normalized L EE [Nicola et al., 2018]. In particular, we have the following: 1. All nodes share the common equilibrium point E k = p, I k = φ(θp), W EI k = W E p−φ −1 (p) φ(W IE p) , k = 1, 2, . . . N . This is the only equilibrium point of the N node system. For convenience, we will refer to this equilibrium point as z 2. The characteristic polynomial for the system (1)-(3) always decomposes into a product of N cubic polynomials, each cubic polynomial is of the same general form, and whose coefficients depend on the eigenvalues r 1 , r 2 , . . . r N of W EE as the following: 3. The cubic polynomial associated with the maximum eigenvalue of L EE induces a Hopf bifurcation of the equilibrium point. While the general form for the bifurcation curve is complicated, it is typically of the form where g(x) is a nonmonotonic, unimodal function ( Figure 1C) [Nicola et al., 2018].

Appendix A3: Normalized Connectivity
Here, we will apply the MSF approach to the Wilson Cowan system in (1)-(3). Rather than deriving an alternate form of the MSF for a general model such as (16), we will proceed directly with the derivation for system (1)-(3). The extension to other models is straightforward. First we note that with the L 1 normalization condition (4), the synchronized solution(E s (t), I s (t), W I s (t)) is a solution to the model for a single recurrently coupled node: and thus is an invariant set of the full N node system for all connectomes with normalization constant W E . Then, consider the variational equations generated by Then, we have the following: where H.O.T. denotes Higher Order Terms. This system written in matrix form becomes: As in the traditional MSF approach, we consider the case where the matrix W EE is diagonalizable: W EE = P DP −1 , η = P −1 , η i = P −1 i, η ω = P −1 ω Then we have: τ 1 η = P −1 −P η + φ (W E E s − W I s I s ) P Dη − W I s · P η i − I s · P η ω (35) Which yield the following N independent, 3-dimensional systems: The important point here is that for some attractor state for fixed W E , (E s (t), I s (t), W I s (t)), we can assess the stability for any weight matrix W EE by analyzing the stability of the equilibrium point of the decoupled blocks z above. These blocks only differ based on the dependence of the eigenvaluesr k of W EE . In fact, for fixed W E we need only do the analysis forr k such that |r k | ≤ W E , as required by the L 1 normalization constraint on W EE . This system in conjunction with the system (27)-(29) is used to compute the Lyapunov spectrum for any possible normalized connectome by varyingr k over a 2D mesh, and computing the maximum Lyapunov exponent over this mesh.