Neuronal Oscillations on Evolving Networks: Dynamics, Damage, Degradation, Decline, Dementia, and Death
NNeuronal Oscillations on Evolving Networks:Dynamics, Damage, Degradation, Decline, Dementia, and Death
Alain Goriely , Ellen Kuhl , and Christian Bick , , Mathematical Institute, University of Oxford, Oxford, OX2 6GG, UK Living Matter Laboratory, Stanford University, Stanford, CA 94305, USA Department of Mathematics, University of Exeter, Exeter, EX4 4QF, UK Institute for Advanced Study, Technische Universit¨at M¨unchen, 85748 Garching, Germany
Neurodegenerative diseases, such as Alzheimer’s or Parkinson’s disease, show characteristic degradation ofstructural brain networks. This degradation eventually leads to changes in the network dynamics and degra-dation of cognitive functions. Here, we model the progression in terms of coupled physical processes: Theaccumulation of toxic proteins, given by a nonlinear reaction-diffusion transport process, yields an evolvingbrain connectome characterized by weighted edges on which a neuronal-mass model evolves. The progressionof the brain functions can be tested by simulating the resting-state activity on the evolving brain network. Weshow that while the evolution of edge weights plays a minor role in the overall progression of the disease,dynamic biomarkers predict a transition over a period of 10 years associated with strong cognitive decline.
PACS numbers: 87.10.Ed, 87.15.hj, 87.16.Ac, 87.19.L-, 87.19.lp, 87.19.xr
Introduction.—
Neurodegenerative diseases are not only ma-jor health and societal problems [1], they are also formidablescientific challenges. Their neuropathology, characterized bydiseased brain tissue and cortical atrophy, is linked to the ac-cumulation of toxic proteins. These structural modifications change the way neurons interact [2] and lead to cognitive de-cline and neurobehavioral symptoms [3]. Specifically, axonaldeath has a direct effect on the collective brain network dy-namics , including synchronization [4], dynamics [5, 6] andconnectivity [7]. There are three interconnected physical andcognitive processes at work: disease progression through thebrain, structural damage created by the disease, and dynamicchanges from damage with the associated functional loss.Here, we build a model that predict both the spatio-temporal evolution of the disease but also how it affects ba-sic cognitive functions. Our approach combines dynamicson multiple temporal scales (years for the disease and sec-onds for the resting-state dynamics) with multiphysics at var-ious levels (transport, aggregation, damage, and oscillations).Specifically, we look at interacting dynamical processes on anevolving network structure: Disease progression changes thestructural properties of the brain connectome, which resultsin changes to characteristic dynamics of the brain networkdynamics. To probe cognitive functions of a given connec-tome we simulate whole-brain resting states against functionalmagnetic resonance imaging (fMRI) resting-state data [8].Since Gamma activity emerge in neural populations [9] andis related to both hippocampal memory formation [10] andAlzheimer’s disease [11, 12], we use a minimal neural-massmodel with intrinsic frequency in the Gamma range (definedby
Γ = [30 , Hz) [13].
Disease progression.—
We follow the prion-likeparadigm [14–16], which proposes that degeneration iscaused by the invasion and conformational autocatalyticconversion of misfolded proteins transported along axonalpathways. [17]. The prion-like idea has served as an impor-tant unifying concept through which many features can be frontalparietaloccipitaltemporallimbicbasal gangliabrain stem
FIG. 1. The connectome G with N = 83 nodes and M = 1654 edges (592 shown). understood such as staging [18], biomarker evolution [19],and neural atrophy. Basic models based in this notion recovermost of these observations [20, 21]. The basic features ofthese diseases can be obtained by restricting all physicalquantities on the structural connectome, a brain networkrepresenting the connections between different regions ofinterest [22–24]. Our model [25, 26], with parameters takenfrom [27], combines network diffusion encoded by a graphLaplacian and a reaction term characterizing the populationamplification due to the conversion of healthy proteins.We model the connectome as an evolving undirectedweighted graph G T at time T ≥ . The initial graph G is extracted from from the tractography of diffusion tensormagnetic resonance images of 418 healthy subjects of the Bu-dapest Reference Connectome v3.0 [28, 29]. Each node k isassociated with a particular brain region R s , s = { , . . . , } ,corresponding to the frontal, parietal, temporal, occipitallobes, the limbic area, the basal ganglia, and the brain stem;cf. Fig. 1. The graph has M edges with weight w kj ( T ) , be-tween nodes k and j , defined as the ratio of the number offibers between the nodes and the mean fiber length, whichleads to the symmetric matrix W = ( w kj ) , k, j ∈ { , . . . , N } a r X i v : . [ n li n . AO ] S e p and the weighted graph Laplacian L = ρ ( D − W ) where D = diag ( (cid:80) Nj =1 w kj ) and ρ is an overall velocity constantdefining the time-scale of transport.Assuming the concentration of healthy proteins remainsmostly unchanged [21], the protein concentration c k ( t ) atnode k obeys a network discretization of a Fisher–KPP equa-tion [25] given by ˙ c k = − N (cid:88) j =1 L kj c j + αc k (1 − c k ) , k = 1 , . . . , N, (1)where α characterizes the conversion from healthy to toxicproteins. Since L depends on the weights W , we require anequation for the evolution of the weights in time. Network damage and evolution.—
The accumulation oftoxic proteins influences the network properties due to its ef-fect on synapses, plasticity, and eventual cell death [2, 30, 31].We quantify the damage at each node by a variable q k ∈ [0 , ( healthy, maximal damage), for which we assume a first-order rate model: ˙ q k = βc k (1 − q k ) , q k (0) = 0 , k = 1 , . . . , N, (2)where β characterizes the protein toxicity. Since transportaway from a node depends on the node’s health, the damage ata node affects the connectivity to other nodes. We assume thatthe relative change of an edge weight depend on the damagesat the nodes with a rate γ according to ˙ w kj = − γw kj ( q k + q j ) , k, j = 1 , . . . , N, (3)and the systems (1–3) form a closed system of N + M or-dinary differential equations. Initially, the concentration oftoxic proteins vanishes at all nodes except at seeding nodesthat is disease-dependent. For illustrative purpose, we use thepropagation of tau proteins as the main source of toxic proteinand seed the system in the entorhinal region [32] by setting c k (0) = 0 for all k except for c (0) = c (0) = 0 . . Ini-tially, the connectome is healthy, with q k (0) = 0 for all k , andthe weights w kj (0) are given by G . Damage does not slow down disease progression.—
Beforeconsidering the resting-state dynamics on an evolving con-nectome, we study the effect of damage on the propagationof the disease. To quantify disease progression, we computethree key structural biomarkers evaluated at a sampled time T :(a) the average concentration C ( T ) = N (cid:80) Nj =1 c j ( T ) , (b) theaverage damage Q ( T ) = N (cid:80) Nj =1 q j ( T ) , and (c) the scaledaverage connection weight W ( T ) = (cid:107) W ( T ) (cid:107) / (cid:107) W (0) (cid:107) ,where (cid:107) W (cid:107) = N − (cid:80) Nk,j =1 w kj . We also compute the av-erage damage Q s in region s by only summing over indicesin R s and normalizing accordingly.While one may expect a slowing down of disease progres-sion as the transport network is affected, the actual overalleffect is negligible as shown in Fig. 2. We first compare C in the absence of damage (solid curves), with the case of se-vere damage over a period of 30 years (dotted curve) and see WW QQ C
Time T (years) FIG. 2. Evolution of averaged toxic concentration and damage. Theblack dotted curve (superimposed with the blue solid curve) is the av-erage concentration in the absence of damage ( β = γ = 0 ) whereasthe solid curve is the case of severe damage ( β = (1 / /year, γ = (1 / /year) leading to a reduction of the connection weightof 50% after 20 years. Even for unrealistic values of the parameters( β = 4 , γ = 2 -dashed curve) with a reduction of 99% of all weightsafter 15 years, the delay in the evolution of the concentration is onlya year ( α = (3 / /year, ρ = 1 / mm/year in all simulations). a negligible difference. Even for unrealistic values (dashedcurves) leading to the destruction of the network within a pe-riod of 15 years, the delay in invasion is only about a year.Hence, network damage does not slow down significantly theinvasion of the disease even in extreme cases. Indeed, the re-duction of diffusion associated with damage mostly affects re-gions where the concentration is high. In these regions, mostof the tissue is already damaged and only very little transportis taking place. In terms of front dynamics, the front velocityin Fisher–KPP depends entirely on the asymptotic zero state.Once nodes have been seeded, the local increase in concentra-tion takes place even in the absence of a network. This non-linear effect due to the conversion of healthy to toxic agents isfundamentally different than the diffusion process where thetoxic protein must be carried from its source.Disease staging can be obtained by computing damage ineach region (Fig. 3) [33]. Physical damage first appears in thelimbic region where the disease originates but moves rapidlyto the temporal and parietal lobes. This early invasion canbe understood by looking at the topological properties of thelinearized system [26] (that can be solved explicitly [34] asshown in the Supplemental Material). Note that increase indamage in the limbic region is slower than in other regions andthat by year 13 the total damage for example in the temporallobe is larger than in the limbic region. Eventually, the diseaseinvades all cortical areas. Resting state brain dynamics.—
To test the declining cog-nitive functions of the brain, we focus on resting-state braindynamics. The time scales involved in the process are of theorder of months for the disease and of the order of secondsfor the rest-state activity. Therefore, the disease dynamics isquasi-stationary and at time t = T we consider the connec-tome G T to be constant when probing resting-state activity.As a proof of principle, we consider a simple neural-mass Time T (years) Q s frontalparietaloccipitaltemporallimbicbasal gangliabrain stem FIG. 3. Evolution of damage in different brain regions. Physicaldamage first manifests itself in the limbic region then moved to thetemporal lobe, the basal ganglia, and the parietal and occipital lobesbefore invading all cortical areas. model on each node representing large interacting excitatoryand inhibitory neural populations to approximate a Wilson–Cowan type model [13, 35]. In the absence of coupling, theintrinsic node dynamics are given by a supercritical Hopf bi-furcation. The state of node k is given by z k ∈ C , and, apartfrom an offset, the real part of z k encodes the activity of theexcitatory population and the imaginary part the activity ofthe inhibitory population. For a given network G T with asso-ciated weights W = W ( T ) the neural populations are cou-pled through the amplitude of the excitatory population andmodulated by a sigmoidal function S ( x ) = 1 / (1 + exp( − x )) through the delay differential equation ˙ z k = F ( z k ) + κS Re N (cid:88) j =1 w kj z j ( t − τ kj ) , (4)where F ( z k ) = z (cid:0) λ + iω k − (cid:12)(cid:12) z k (cid:12)(cid:12) (cid:1) with decay λ = − . ,intrinsic frequencies ω k = ω + δ k = 40Hz + δ k whose devi-ations δ k are sampled from a normal distribution (mean zero,variance . ), coupling gain κ = 10 , and delays τ kj pro-portional to the distance between node k and j from the con-nectome data with transmission speed of . / s (discretizedto have a maximum of 40 distinct delays); these model pa-rameters were chosen to approximate the neural-mass modelin [13] validated against resting-state fMRI. For the initialgraph G , we observe collective oscillations. In the absence ofcoupling, all amplitudes decay exponentially with a frequencyof around 40Hz. Global cognitive decline after physical damage—
As in-dicators of cognitive processes, we consider three dynamicbiomarkers obtained from t sim = 10s of resting-state dynam-ics (4): (d) the overall power in the Gamma-range P ( T ) = (cid:82) Γ PSD ( (cid:104) z (cid:105) )(Ω) d Ω , where PSD is the power spectral den-sity of the signal (cid:104) z (cid:105) = N (cid:80) Nj =1 Re( z j ) , (e) the average os-cillatory activity A ( T ) = N − (cid:80) Nj =1 t − sim (cid:82) t sim | z j ( t ) | d t , and(f) the metastability index B ( T ) = N − (cid:80) Nj =1 σ t ( | z j ( t ) | ) ,where σ t is the variance of the signal over the time interval FIG. 4. Dynamic biomarkers are stable up to year 10 followed by arapid decline from year 16 onward. (a) power in the Γ band; (b) oscil-lation mean amplitude; there are no oscillations in the absence of net-work coupling. (c) mean variability which indicates non-stationary(and potentially metastable) brain dynamics. Mean and standard de-viation for 12 realizations with different intrinsic frequencies and ini-tial conditions (solid gray line is C ( T ) for comparison). [0 , t sim ] . We define the corresponding measures P s , A s , B s for R s by summing over the corresponding nodes and normal-izing accordingly. These dynamic biomarkers have been as-sociated with cognitive processes and neurodegenerative dis-eases [4]. The average amplitude is a measure of the generalactivity and the metastability index is associated with infor-mation processing [36–39].To evaluate how the disease affects the dynamics, wesolved (4) on the evolving brain connectome for differenttimes T ; see Supplemental Material for details. Fig. 4 showsthe mean and standard deviations of the dynamical biomark-ers scaled with respect to the healthy response at T = 0 . Wesee that all dynamical indicators remain fairly unchanged upto year 20 as the disease progresses. At that time, the brain hasalready suffered significant physical damage even if this dam-age cannot be easily assessed from the dynamics (see Fig. 2).It suggests that the brain state is structurally stable againstdamages for an extended time. However, after 20 years, thedynamics undergoes a clear transition when nodes are unableto sustain oscillatory activities.The local dynamical biomarkers show a more differentiatedpicture of how the spreading disease changes the dynamics.The transition of the global dynamics is preceded by a de-cay of oscillatory activity in the temporal lobe. In terms ofthe oscillatory activity, the markers for all other regions areclose to the global mean (Fig. 4b). Given that damage first FIG. 5. Oscillation amplitude changes with the homeostasis pa-rameter. Increasing homeostasis does not lead to a shift in the onsetof decline but only changes the shape of the decline. The globalmean oscillation amplitude (gray lines) shows how the slope variesfor different homeostasis parameters; in the extreme case ξ = 1 , theoscillation amplitude remains constant. By contrast, the mean am-plitude in the temporal lobe (blue lines) decreases in any case. Forcomparison, C ( T ) is given by the solid gray line. accumulates in the limbic system (Fig. 3), this observationmay be counterintuitive. However, it indicates that both speedof damage accumulation and network structure determine thecritical threshold for the decay of oscillatory dynamics: Theregion where one first sees significant structural damage canbe distinct from the region that first undergoes a dynamicaltransition. This is consistent with experimental findings thatindicate that temporal lobe activity is a precursor for the dis-ease onset [40]. Adaptation slows cognitive decline.—
Connections betweenneurons are typically not static but adjust in response to thedynamics [41, 42]. Thus, in addition to damage-inducedchanges to the network, the neural dynamics itself affects thenetwork to maintain homeostasis [43]. To capture this effect,we implement a minimal model of homeostatic adaptation thataims to keep the mean connectome coupling constant overtime. For an adaptation parameter ξ ∈ [0 , , we define thescaled matrix with W (0) = W (0) and W ( T ) = (cid:18) (1 − ξ ) + ξ (cid:18) (cid:107) W ( T − (cid:107)(cid:107) W ( T ) (cid:107) (cid:19)(cid:19) W ( T ) (5)for T = 1 , , . . . . This homeostatic adaptation creates a newweighted connectome for the resting-state dynamics (4) forevery year T of disease progression. Through (5), home-ostasis acts by globally rescaling the coupling weights. For ξ = 0 , there is no homeostatic adaptation and the networkstructure changes solely through the disease progression. For ξ = 1 , there is complete homeostatic adaptation in the sensethat we have a constant mean coupling weight (cid:107) W ( T ) (cid:107) = 1 for all T . This does not imply W ( T ) = W (0) . Rather, dueto the global nature of the homeostatic adaptation, disease-induced changes in coupling strength in one brain region willyield hyper-excitability in another brain region to balance theoverall decay in coupling strength.Homeostatic adaptation modulates the disease progression.We simulated the resting-state dynamics (4) subject to thehomeostatic adaptation (5). Fig. 5 shows the average ampli-tude A ( T ) as disease progresses for different values of theadaptation parameter ξ . An increase in homeostasis does not change the actual onset of loss of oscillatory activity (aroundyear 13) but yields a slow-down of the degeneration. How-ever, assuming that cognitive decline sets in once a certainthreshold is reached, the slow-down of disease progressionwill alter the onset of the overall transition. Fig. 5 also showsmean oscillation amplitude for the temporal lobe. The tem-poral lobe shows a decay in oscillatory activity for all val-ues of the adaptation parameters preceding the overall decayof oscillatory activity. This happens even for full adaptation, ξ = 1 , where the overall oscillatory activity stays constant inthe time-window considered here. This decline implies thatother brain areas have to be upregulated in order to keep theoverall activity almost constant. However, it also means thatfor our simple model of adaptation, the variation of oscilla-tion amplitude is a precursor for the overall transition of thedynamics independently of the adaptation parameter. Discussion.—
Compared to previous studies that focused oneither static properties of the declining network [44], synchro-nization [4] or activity-based decline [45], we focused on theimportance of the underlying interacting physical processes.The physical damage of the brain and neural dynamics—twoindependent yet coupled processes—interact on the same con-nectome. Neurodegeneration, modeled as an invasion processdue to the accumulation of toxic proteins, provides a natu-ral evolution of the connectome on long time scales that canbe probed dynamically on short time scales. Our observa-tions are compatible with the activity-dependent spreading hy-pothesis [45]: toxic proteins will spread predominantly alonghighly connected nodes which also show high activity due tothe amount of input they receive.While the actual brain connectome is not an undirectedgraph, there is very little information on directionality forthe human brain. However, an analysis of a mesocale mouseconnectome [46] that has been used for toxic protein diffu-sion [24, 47] reveals that asymmetry delays significantly theonset of the disease but preserves its main characteristics (seeSupplemental Material).Our setup provides a unified framework that combines thebiophysics of disease spreading with whole-brain dynamics togive mechanistic insights into the dynamics of neurodegener-ative diseases and its associated cognitive decline. First, weshowed that damage does not slow the disease propagation asdamage is delayed with respect to seeding. Very low levelof toxic proteins diffuse and seed new regions. Then, evenin the absence of transport, there is a local autocatalytic in-crease of toxic protein. Second, we gained insight into the dy-namical transitions appearing in some brain regions and con-firms the prediction [40] that the temporal lobe is one of thefirst to see alteration in brain dynamics, hence showing cog-nitive deficiencies related to that region for Alzheimer’s dis-ease as shown in Fig. 4. Third, our results elucidate the in-terplay between network adaptation and spreading. We foundthat incorporating a simple model of global homeostasis doesnot change the onset of dynamical changes, but how fast theyevolve as the disease progresses. Interestingly, we saw that, inthe particular case of Alzheimer disease, a decline of oscilla-tions in the temporal lobe is a universal indicator independentof the adaptation parameter.Further insights are needed on how the multiphysics ofdisease propagation interact with brain network dynamics toidentify reliable noninvasive biomarkers to assess disease pro-gression as early as possible and develop an integrated ap-proach for treatment. Other models for whole-brain dynam-ics are available that focus on different features of brain dy-namics [8] and relate to microscopic neural properties [48].When combined with data from a comprehensive longitudi-nal study and suitable dynamical models, our approach willbe able to explore new treatment approaches and interventionscenarios on the combined level of network structure and dy-namics [2, 49, 50].This work was supported by the Engineering and PhysicalSciences Research Council grant EP/R020205/1 to AG and bythe National Science Foundation grant CMMI 1727268 to EK. [1] Y. H. El-Hayek, R. E. Wiley, C. P. Khoury, R. P. Daya, C. Bal-lard, A. R. Evans, M. Karran, J. L. Molinuevo, M. Norton, andA. Atri, J Alzheimers Dis , 323 (2019).[2] J. J. Palop and L. Mucke, Nat Neurosci , 812 (2010).[3] K. Koch, N. E. Myers, J. G¨ottler, L. Pasquini, T. Grimmer,S. F¨orster, A. Manoliu, J. Neitzel, A. Kurz, H. F¨orstl, V. Riedl,A. M. Wohlschl¨ager, A. Drzezga, and C. Sorg, Cereb Cortex , 4678 (2015).[4] P. J. Uhlhaas and W. Singer, Neuron , 155 (2006).[5] J. Zimmermann, A. Perry, M. Breakspear, M. Schirner,P. Sachdev, W. Wen, N. A. Kochan, M. Mapstone, P. Ritter,A. R. McIntosh, and A. Solodkin, NeuroImage-Clin , 240(2018).[6] R. C. Budzinski, B. R. R. Boaretto, T. L. Prado, and S. R.Lopes, Phys Rev E , 022402 (2019).[7] K. Wang, M. Liang, L. Wang, L. Tian, X. Zhang, K. Li, andT. Jiang, Hum Brain Mapp , 967 (2007).[8] M. Breakspear, Nat Neurosci , 340 (2017).[9] A. Fisahn, F. G. Pike, E. H. Buhl, and O. Paulsen, Nature ,186 (1998).[10] J. Chrobak and G. Buzs´aki, J Neurosci , 388 (1998).[11] H. F. Iaccarino, A. C. Singer, A. J. Martorell, A. Rudenko,F. Gao, T. Z. Gillingham, H. Mathys, J. Seo, O. Kritskiy, F. Ab-durrob, et al. , Nature , 230 (2016).[12] B. McDermott, E. Porter, D. Hughes, B. McGinley, M. Lang,M. O’Halloran, and M. Jones, J Alzheimers Dis , 363 (2018).[13] G. Deco, V. Jirsa, A. R. McIntosh, O. Sporns, and R. K¨otter,Proc Nat Acad Sci USA , 10302 (2009).[14] M. Jucker and L. C. Walker, Ann Neurol , 532 (2011).[15] L. C. Walker and M. Jucker, Annu Rev Neurosci , 87 (2015).[16] M. Goedert, Science , 1255555 (2015).[17] S. B. Prusiner, Proc Nat Acad Sci USA , 13363 (1998).[18] H. Braak and E. Braak, Acta Neuropathologica , 239 (1991).[19] C. R. Jack and D. M. Holtzman, Neuron , 1347 (2013).[20] J. Weickenmeier, E. Kuhl, and A. Goriely, Phys Rev Lett ,158101 (2018).[21] J. Weickenmeier, M. Jucker, A. Goriely, and E. Kuhl, J Mech Phys Solids , 264 (2019).[22] A. Raj, A. Kuceyeski, and M. Weiner, Neuron , 1204 (2012).[23] Y. Iturria-Medina, R. C. Sotero, P. J. Toussaint, A. C. Evans,and A. D. N. Initiative, PLoS Comput Biol , e1003956(2014).[24] M. X. Henderson, E. J. Cornblath, A. Darwich, B. Zhang,H. Brown, R. J. Gathagan, R. M. Sandler, D. S. Bassett, J. Q.Trojanowski, and V. M. Lee, Nat Neurosci , 1248 (2019).[25] S. Fornari, A. Sch¨afer, A. Goriely, and E. Kuhl, J R Soc Inter-face , 20190356 (2019).[26] S. Fornari, A. Sch¨afer, E. Kuhl, and A. Goriely, J Theor Biol , 110102 (2020).[27] F. Kundel, L. Hong, B. Falcon, W. A. McEwan, T. C. Michaels,G. Meisl, N. Esteras, A. Y. Abramov, T. J. Knowles, M. Goed-ert, et al. , ACS chemical neuroscience , 1276 (2018).[28] J. A. McNab, B. L. Edlow, T. Witzel, S. Y. Huang, H. Bhat,K. Heberlein, T. Feiweier, K. Liu, B. Keil, J. Cohen-Adad, et al. , NeuroImage , 234 (2013).[29] B. Szalkai, C. Kerepesi, B. Varga, and V. Grolmusz, CognitiveNeurodynamics , 113 (2017).[30] J. Collinge, M. A. Whittington, K. C. L. Sidle, C. J. Smith,M. S. Palmer, A. R. Clarke, and J. G. R. Jefferys, Nature ,295 (1994).[31] J. J. Palop and L. Mucke, Nat Rev Neurosci , 777 (2016).[32] A. De Calignon, M. Polydoro, M. Su´arez-Calvet, C. William,D. H. Adamowicz, K. J. Kopeikina, R. Pitstick, N. Sahara, K. H.Ashe, G. A. Carlson, et al. , Neuron , 685 (2012).[33] J. P. Kesslak, O. Nalcioglu, and C. W. Cotman, Neurology ,51 (1991).[34] A. Goriely, Integrability and Nonintegrability of DynamicalSystems (World Scientific Publishing Company, 2001).[35] G. Deco, J. Cabral, M. W. Woolrich, A. B. Stevner, T. J.van Hartevelt, and M. L. Kringelbach, NeuroImage , 538(2017).[36] M. I. Rabinovich, V. S. Afraimovich, C. Bick, and P. Varona,Phys Life Rev , 51 (2012).[37] E. Tognoli and J. A. S. Kelso, Neuron , 35 (2014).[38] S. Heitmann and M. Breakspear, Net Neurosci , 150 (2018).[39] M. Shanahan, Chaos , 013108 (2010).[40] L. L. Beason-Held, J. O. Goh, Y. An, M. A. Kraut, R. J.O’Brien, L. Ferrucci, and S. M. Resnick, J Neurosci , 18008(2013).[41] P. J. Sj¨ostr¨om, E. A. Rancz, A. Roth, and M. H¨ausser, PhysiolRev , 769 (2008).[42] H. Markram, W. Gerstner, and P. J. Sj¨ostr¨om, Front Syn Neu-rosci , 4 (2011).[43] S. Frere and I. Slutsky, Neuron , 32 (2018).[44] C. J. Stam, B. Jones, G. Nolte, M. Breakspear, and P. Scheltens,Cereb Cortex , 92 (2006).[45] W. de Haan, K. Mott, E. C. van Straaten, P. Scheltens, and C. J.Stam, PLoS Comput Biol , e1002582 (2012).[46] S. W. Oh, J. A. Harris, L. Ng, B. Winslow, N. Cain, S. Mihalas,Q. Wang, C. Lau, L. Kuan, A. M. Henry, et al. , Nature , 207(2014).[47] M. X. Henderson, S. Sedor, I. McGeary, E. J. Cornblath,C. Peng, D. M. Riddle, H. L. Li, B. Zhang, H. J. Brown, M. F.Olufemi, et al. , Neuron , 822 (2020).[48] C. Bick, M. Goodfellow, C. R. Laing, and E. A. Martens, JMath Neurosci , 9 (2020).[49] W. de Haan, E. C. W. van Straaten, A. A. Gouw, and C. J. Stam,PLoS Comput Biol , e1005707 (2017).[50] R. G. Canter, J. Penney, and L.-h. Tsai, Nature539