Complex spatiotemporal oscillations emerge from transverse instabilities in large-scale brain networks

RESEARCH ARTICLE Complex spatiotemporal oscillations emerge from transverse instabilities in large-scale brain networks Pau ClusellaID 1 *, Gustavo Deco 2,3 , Morten L. Kringelbach 4,5 , Giulio Ruffini 6 , Jordi GarciaOjalvoID 1 1Department of Medicine and Life Sciences, Universitat Pompeu Fabra, Barcelona, Spain, 2Computational Neuroscience Group, Center for Brain and Cognition, Department of Information and Communication Technologies, Universitat Pompeu Fabra, Barcelona, Spain, 3Institucio ´Catalana de Recerca i Estudis Avanc¸ats (ICREA), Barcelona, Spain, 4Department of Psychiatry, University of Oxford, Oxford, United Kingdom, 5Center for Music in the Brain, Department of Clinical Medicine, Aarhus University, Aarhus, Denmark, 6Brain Modeling Department, Neuroelectrics, Barcelona, Spain *[email protected] Abstract Spatiotemporal oscillations underlie all cognitive brain functions. Large-scale brain models, constrained by neuroimaging data, aim to trace the principles underlying such macroscopic neural activity from the intricate and multi-scale structure of the brain. Despite substantial progress in the field, many aspects about the mechanisms behind the onset of spatiotemporal neural dynamics are still unknown. In this work we establish a simple framework for the emergence of complex brain dynamics, including high-dimensional chaos and travelling waves. The model consists of a complex network of 90 brain regions, whose structural connectivity is obtained from tractography data. The activity of each brain area is governed by a Jansen neural mass model and we normalize the total input received by each node so it amounts the same across all brain areas. This assumption allows for the existence of an homogeneous invariant manifold, i.e., a set of different stationary and oscillatory states in which all nodes behave identically. Stability analysis of these homogeneous solutions unveils a transverse instability of the synchronized state, which gives rise to different types of spatiotemporal dynamics, such as chaotic alpha activity. Additionally, we illustrate the ubiquity of this route towards complex spatiotemporal activity in a network of next generation neural mass models. Altogehter, our results unveil the bifurcation landscape that underlies the emergence of function from structure in the brain. Author summary Monitoring brain activity with techniques such as electroencephalogram (EEG) and functional magnetic resonance imaging (fMRI) has revealed that normal brain function is characterized by complex spatiotemporal dynamics. This behavior is well captured by large-scale brain models that incorporate structural connectivity data obtained with MRIbased tractography methods. Nonetheless, it is not yet clear how these complex dynamics PLOS COMPUTATIONAL BIOLOGY PLOS Computational Biology | https://doi.org/10.1371/journal.pcbi.1010781 April 12, 2023 1 / 34 a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 OPEN ACCESS Citation: Clusella P, Deco G, Kringelbach ML, Ruffini G, Garcia-Ojalvo J (2023) Complex spatiotemporal oscillations emerge from transverse instabilities in large-scale brain networks. PLoS Comput Biol 19(4): e1010781. https://doi.org/10.1371/journal.pcbi.1010781 Editor: Boris S. Gutkin, E ´cole Normale Supe ´rieure, College de France, CNRS, FRANCE Received: December 1, 2022 Accepted: March 24, 2023 Published: April 12, 2023 Peer Review History: PLOS recognizes the benefits of transparency in the peer review process; therefore, we enable the publication of all of the content of peer review and author responses alongside final, published articles. The editorial history of this article is available here: https://doi.org/10.1371/journal.pcbi.1010781 Copyright: ©2023 Clusella et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. Data Availability Statement: All code and data in support of this publication are publicly available at https://github.com/pclus/transverse-instabilities. emerge from the interplay of the different brain regions. In this paper we show that complex spatiotemporal dynamics, including travelling waves and high-dimensional chaos can arise in simple large-scale brain models through the destabilization of a synchronized oscillatory state. Such transverse instabilities are akin to those observed in chemical reactions and turbulence, and allow for a semi-analytical treatment that uncovers the overall dynamical landscape of the system. Overall, our work establishes and characterizes a general route towards spatiotemporal oscillations in large-scale brain models. Introduction The interplay between spiking neurons across the brain produces collective rhythmic behavior at multiple frequencies and spatial resolutions [1,2]. This oscillatory neural activity is fundamental for proper cognitive function [3,4], and is reflected in a plethora of spatiotemporal phenomena in recorded signals [5–8]. At the microscopic scale, computational and mathematical models have successfully captured some of these emergent properties by analyzing the collective behavior of networks of coupled neurons [1,9–11]. In larger scales, the task becomes increasingly challenging, as one needs to model several populations of neurons, which increases the mathematical complexity and the computational cost of the problem. Mesoscale neural-mass models (NMMs) allow to overcome this situation by capturing the neural activity of large numbers of neurons with a few equations [12–17]. Recent advances on neuroimaging allow to characterize the structural organization of the brain as a complex network [18–20]. In this representation, each node of the network corresponds to a brain region composed of densely inter-connected neurons, and edges across nodes represent pairwise interactions across distant regions. Combining NMMs with connectomics data one can create large-scale brain models whose dynamical properties reflect the principles underlying macroscopic neural activity [21–24]. This framework has been used, for instance, to unveil the nature of resting-state fluctuations [25–30], investigate the relation between structural and functional connectivity [31–35], and characterize transitions between brain states [36–39]. In clinical applications, large-scale brain models allow for studying macroscopic aspects of neurophatologies such as epilepsy or Alzheimer disease [40–43], and simulate the use non-invasive brain stimulation protocols for potential treatments [44,45]. Despite all this progress, many aspects about the mechanisms through which large-scale brain models reproduce macroscopic neural dynamics are still unknown. For instance, synaptic delays between distant regions [25,46], noise [25,29], and heterogeneities [28] are usually acknowledged as a source of dynamical complexity. Nonetheless, spatiotemporal behavior can also arise from homogeneous deterministic systems with instantaneous interactions [30,35, 44]. In particular, Forrester et. al. [35] recently investigated the dynamics of a large-scale brain model composed of weakly-coupled Jansen’s NMMs. By means of a phase-reduction approximation, they unveil phase-locked states emerging from a instability of the synchronized state (see also [33]). For arbitrary coupling values, however, bifurcation studies of whole-brain networks are rather limited to numerical investigations [30,44]. Other studies show that systems of NMMs coupled through simplified network topologies might display travelling waves and even chaotic dynamics [47–50]. Whether these results translate to irregular brain networks remains, so far, unexplored. In this paper we characterize the onset of spatiotemporal dynamics in a simple large-scale brain model without heterogeneities, noise, nor delays. In close analogy to pattern-formation mechanisms in reaction-diffusion systems [51–55], coupling among brain regions alone is PLOS COMPUTATIONAL BIOLOGY Transverse instabilities in large-scale brain networks PLOS Computational Biology | https://doi.org/10.1371/journal.pcbi.1010781 April 12, 2023 2 / 34 Funding: PC, GD, GR, and JGO have received funding from the Future and Emerging Technologies Programme (FET) of the European Union’s Horizon 2020 research and innovation programme (project NEUROTWIN, grant agreement No 101017716). JGO also acknowledges financial support from the Spanish Ministry of Science and Innovation and FEDER (grant PID2021-127311NB-I00), by the “Maria de Maeztu” Programme for Units of Excellence in R&D (grant CEX2018-000792-M), and by the Generalitat de Catalunya (ICREA Academia programme). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. Competing interests: I have read the journal’s policy and the authors of this manuscript have the following competing interests: GR is a co-founder of Neuroelectrics, a Company that manufactures tES and EEG technology. The remaining authors don’t have any conflict of interest. enough to spontaneously destabilize an underlying synchronized state, thereby generating complex oscillatory behavior. Our analysis consists of two parts: First we show that, given a common normalization on the incoming input of each region [56,57], the network possesses an invariant homogeneous manifold, i.e., a set of states in which the behavior of each node is identical across all the network. These states are described by a low-dimensional system: a selfcoupled version of the NMM used for the evolution of each brain region. Bifurcation analysis of the system reveals how different system parameters modify the onset of synchronized oscillatory states within the manifold. Second, we employ the Master Stability Function (MSF) formalism [58–62] to investigate the stability of the homogeneous states to heterogeneous perturbations. The synchronized oscillatory solution of the system turns out to be transversally unstable in a large region of parameter space, giving rise to complex spatiotemporal dynamics including travelling waves, multistability, and high-dimensional chaos. In order to illustrate the ubiquity of this mechanism, we mainly use the Jansen NMM for the single-node dynamics [15,16,63], but also briefly review the onset of chaotic dynamics using a next generation NMM [17]. Moreover, we also show that, even when the normalization condition is not fulfilled, the bifurcation diagram of the system remains similar to that of the simplified model. Our work extends previous findings on weakly-coupled NMMs [33,35] and simplified network topologies [47–50] to a comprehensive characterization of the emerging spatiotemporal behavior. Results The large-scale brain model We build a large-scale brain model using a structural connectivity (SC) network comprised of N= 90 brain regions defined by the Automated Anatomical Labeling parcellation (AAL-90) [64], which includes 76 cortical and 14 subcortical regions. Each brain region is represented by a network node, and pairwise interactions between regions are given by a row-normalized connectivity matrix ~ W¼ ð~ wijÞwhere i,j= 1, . . .,N. The connection weights ~ wij are non-negative quantities that indicate the synaptic strength (average number of connections) from region jto region i, and have been obtained from tractography data from 16 human subjects (see [65] and also Methods). Although networks obtained from DTI are usually symmetric prior rownormalization, our approach can also be applied to asymmetric networks (see Methods). In order to model the dynamics of each brain region we use, for most of the paper, Jansen’s model for a cortical column [15,16,63]. According to this model, the behavior of each brain region is given by a system of six ordinary differential equations (Eq (9) in Methods) that account for the interactions between a population of excitatory pyramidal neurons (PNs), a population of inhibitory interneurons (INs), and recurrent connections within pyramidal neurons (rPNs). Within each region, the PNs receive two sources of external input: a baseline firing rate p, which we assume constant and identical across the network, and the incoming firing rates from the PNs of other brain regions, modulated by a coupling parameter �. Hence, long range connections established by the structural connectivity matrix are all assumed excitatory, whereas inhibition acts only locally within each network node. Overall, we obtain a large-scale brain model composed of 90 ×6 = 540 ODEs (see Eq (10) in Methods). Despite its high dimensionality, this system is relatively simple to analyze, as it does not include noise nor time delays and its parameters are assumed to be identical across brain regions. Similar large-scale brain models based on the Jansen system have been analyzed in previous works, the main difference being the connectivity data used for the underlying PLOS COMPUTATIONAL BIOLOGY Transverse instabilities in large-scale brain networks PLOS Computational Biology | https://doi.org/10.1371/journal.pcbi.1010781 April 12, 2023 3 / 34 network topology [35,44]. Our work extends the results of these previous studies, providing a comprehensive bifurcation landscape of the network. Bifurcation diagram Following a customary approach in the literature [30,33,35,44,45,56,57], our definition of the structural connectome matrix ~ Westablishes that all network nodes receive the same amount of input, i.e., PN j¼1~ wij ¼1for all i= 1, . . .,N. This condition ensures that the system has homogeneous (uniform) states, i.e., states in which all network nodes evolve identically. If the system is initialized in such a state, it will remain there forever unless perturbed: mathematically speaking, these states lie on an invariant manifold. The existence of homogeneous solutions does not prevent the existence of heterogeneous states, in which the trajectories of different nodes evolve differently. If a homogeneous state proves to be unstable to arbitrary (heterogeneous) perturbations, complex spatiotemporal dynamics might arise. In this section we investigate the homogeneous states of the large-scale brain model systematically, as a means to unveil the emergence of non-trivial heterogeneous patterns. This analysis requires two steps: First, we identify the homogeneous states of the system and their stability to uniform perturbations, i.e. we study the homogeneous invariant manifold of the system. Then, we analyze the stability of these states to arbitrary perturbations using the Master Stability Function (MSF) formalism [58]. Dynamics of the homogeneous invariant manifold. If all brain regions behave identically, then each network node follows the dynamics of the self-coupled Jansen model given by Eq (18) in Methods. This system defines the homogeneous invariant manifold of the system, hence its stable states correspond to those homogeneous states of the full model (Eq 10) that are stable to uniform perturbations. In this section we analyze the dynamics of Eq (18) alone. Therefore, the term stability in this section refers to stability within the homogeneous manifold only. We uncover the different attractors of the self-coupled system with the help of the bifurcation analysis software AUTO-07p [66], using the common external input p�0 and the coupling strength ��0 as control parameters. Despite both parameters being positive quantities, we also include negative values of pin the bifurcation diagrams in order to reveal the entire bifurcation structure. We start by recalling the dynamics of the uncoupled Jansen model (for details readers can refer to [63]). Fig 1(a) shows the value of the mean membrane potential of the PN population v corresponding to different solutions of system (18) for fixed �= 0 and varying baseline input p. For small positive p, two stable steady states coexist (dark pink curves), each of them leading to a different type of stable periodic state (black curves) upon increasing the baseline input. First, the high-activity fixed point (upper pink branch) undergoes a supercritical Hopf bifurcation (HBþ 1) at p�90, which corresponds to the onset of alpha oscillatory activity (�10Hz). This periodic state persists until p�315, where it vanishes through a second supercritical Hopfbifurcation (HBþ 2) leading again to a stable high-activity steady state. Second, the low-activity stationary state for psmall (lower dark pink branch) vanishes through a saddle-node in a invariant cycle (SNIC) bifurcation for p�114, giving rise to finite amplitude oscillations at a theta range (�4Hz). This spiky oscillatory activity vanishes at p�137 through a fold (or saddle-node) bifurcation of limit-cycles (FLC 1 ). The addition of coupling (� > 0) modifies this bifurcation scenario. For instance, Fig 1(b) shows the bifurcation diagram for �= 4. The figure shows that coupling eliminates the supercritical Hopf bifurcation HBþ 1, so that for small values of ponly the low-activity fixed point is stable. As in the uncoupled case, this steady state vanishes at p�111 through a SNIC bifurcation. However, now the branch of oscillatory states arising from the SNIC connects both the PLOS COMPUTATIONAL BIOLOGY Transverse instabilities in large-scale brain networks PLOS Computational Biology | https://doi.org/10.1371/journal.pcbi.1010781 April 12, 2023 4 / 34 theta (111 ≲p≲150) and alpha (134 ≲p≲351) frequency bands, which coexist in a small region bounded by two FLCs. Hence, for small but non-zero coupling the alpha oscillatory state emerges through the fold of cycles FLC 2 instead of a Hopf bifurcation, and, as in the uncoupled case, it disappears through the HBþ 2supercritical Hopf bifurcation for large input (p�351). Further increase of �leads to a total disappearance of bistability between oscillatory states. Fig 1(c) shows an example of this situation for �= 50. Upon increasing the external input p, the low activity fixed point leads to finite-amplitude oscillations through the SNIC bifurcation (p�84.68). This oscillatory state becomes the only attractor for a wide range of p, until it vanishes at the supercritical Hopf bifurcation HBþ 2(p�330), beyond which trajectories converge to the high-activity fixed point. The frequency of this single stable oscillatory state (see Fig 1(d)) shows two distinct plateaus that dominate for most values of p, located at around the theta (3–5Hz) and alpha (7–9Hz) levels. Overall, the diagrams in Fig 1(a) and (c) reveal a change on the onset of alpha oscillatory activity, together with a loss of bistability as �increases. These changes can be better understood from the two-parameter bifurcation analysis shown in Fig 1(e) and (f). These diagrams illustrate the different regions of stability of Eq (18) for varying pand �, with panel (f) showing a zoomed version of panel (e) for small coupling. The red continuous curve in Fig 1(e) and (f) corresponds to the Hopf bifurcation HBþ 1giving raise to alpha activity for low values of �(see Fig 1(a)). This bifurcation becomes subcritical (dashed red line) through a Bautin (or generalized Hopf) codimension-2 bifurcation at p�40, ��1.8. At this point, the FLC 2 appears Fig 1. Bifurcation diagrams of the homogeneous states subject to uniform perturbations. (a-c) One-parameter bifurcation diagrams of Eq (18) obtained by varying pand with fixed �= 0 (a), �= 4 (b), and �= 50 (c). Colored curves indicate stable fixed points (dark pink), unstable fixed points (light pink), extrema of stable limit-cycles (black), and extrema of unstable limit-cycles (grey). Dashed black vertical lines indicate the relevant bifurcation points cited in the text. (d) Frequency of the stable oscillatory state by varying pand with fixed �= 50. (e,f) Two-parameter bifurcation diagram of Eq (18) depending on the external input pand the coupling strength �. Panel (f) is a zoomed version of (e). Curves indicate different bifurcation types: supercritical Hopf (HB + , continuous red), subcritical Hopf (HB − , dashed red), saddle-node (SN, brown), saddle-node in a invariant cycle (SNIC, dark blue), saddle-node of limit-cycles (FLC, green), and homoclinic (Hom., orange) The light-blue region indicates the existence of a single periodic state. The light-blue region with stripped black pattern indicates the coexistence between a limit-cycle and a fixed point. The dark-blue region indicates the coexistence of two stable periodic states. All results obtained by analyzing the system of Eq (18) using the bifurcation analysis software AUTO-07p (scripts for one-parameter bifurcations available at www.github.com/pclus/transverse-instabilities). https://doi.org/10.1371/journal.pcbi.1010781.g001 PLOS COMPUTATIONAL BIOLOGY Transverse instabilities in large-scale brain networks PLOS Computational Biology | https://doi.org/10.1371/journal.pcbi.1010781 April 12, 2023 5 / 34 (green line), and joints the FLC 1 at a cusp (p�165, ��10) at which both bifurcations vanish. Therefore, within the dark-blue region delimited by the two FLCs (green curves) and the SNIC bifurcations (dark-blue curve) two stable oscillatory states coexist. Beyond this region, and for a wide range of �values, the dynamical landscape of the system becomes simpler and coincides with that shown in Fig 1(c) (see Fig 1(e) for �= 50). This scenario changes for very high coupling, when the SNIC bifurcation turns to a SN through a saddle-node separatrix loop (SNSL) codimension-2 bifurcation [67]. From this point, a homoclinic bifurcation (Hom.) bounds the region of oscillatory dynamics, which appear for arbitrary low p. Parallel to the homoclinic line, two additional branches of FLC appear, leading to a very narrow region of oscillatory bistability. Since they have a minor effect on the overall bifurcation landscape, we do not depict them in Fig 1(e). Finally, further increase of �ceases all oscillatory activity through the HBþ 2(red continuous curve). In summary, the dynamics of the system within the homogeneous invariant manifold can be divided in three main regions: • Low activity fixed point, represented by the white region to the left of the SNIC bifurcation (dark blue curve) in Fig 1(e). • High activity fixed point, represented by the white region to the right and above of the HBþ 2 line (red curve) in Fig 1(e). • Oscillatory activity at the theta or alpha ranges (or bistability between both), mostly delimited by the HBþ 2and the SNIC bifurcation curves (light blue shadowed region) in Fig 1(e). Overall, the bifurcation analysis of Eq (18) uncovers the rich dynamical repertoire of the homogeneous manifold of the full system (Eq 10). Nonetheless, we remark again that this analysis only concerns homogeneous states, i.e., the different regions of stability and instability outlined so far only assume perturbations that act identically at each network node. In the next section we assess the stability of the homogeneous solutions to heterogeneous perturbations. Transverse stability. The analysis of system (18) discussed above reveals the loci of all homogeneous states of the network model that are stable to uniform perturbations. In order to determine the stability of these states to arbitrary non-uniform perturbations, we follow a wellestablished approach that can be applied to either fixed points or periodic states (see Methods). This technique, analogous to the one used in the study of Turing bifurcations in complex networks [54,68,69], consists on decomposing an arbitrary perturbation vector on the basis given by the eigenvectors of a suitable matrix representing the way the nodes are coupled. This provides a dispersion relation for the growth rate of the perturbations. In our case, instead of the Laplacian matrix used in diffusively coupled systems, we diagonalize the normalized structural connectivity matrix ~ W. As we will see, stable fixed points in the homogeneous manifold remain always stable to heterogeneous perturbations in the full-model. Hence, our focus will be on the stability of the limit-cycle solutions, for which this decomposition technique is known as the Master Stability Function (MSF) [58,59,61]. Fig 2(a) shows the first 5 eigenvectors of ~ Win terms of their components across the 90 brain regions. The first eigenmode is homogeneous, and its associated eigenvalue is always Λ 1 = 1 (see Methods). Perturbations along this direction are the ones already accounted by the analysis of Eq (18) in the previous section. The subsequent eigenmodes are heterogeneous and their associated eigenvalues are between -1 and 1 (see Methods). Perturbations along these directions are transverse to the homogeneous invariant manifold. Despite stemming from an irregular network topology, the eigenvectors exhibit a well defined spatial structure, which can be traced back to the well-known exponential decay of inter-region connectivity with PLOS COMPUTATIONAL BIOLOGY Transverse instabilities in large-scale brain networks PLOS Computational Biology | https://doi.org/10.1371/journal.pcbi.1010781 April 12, 2023 6 / 34 Euclidean distance observed in brain connectivity data obtained with diffusion tensor imaging (DTI) [70,71]. Following the technique explained in Methods, we found that no transverse instabilities arise from the homogeneous fixed points. Therefore next we focus on the stability of the limitcycle solutions by means of the Master Stability Function (MSF) [58,59,61]. The growth rate of an infinitesimal perturbation of a periodic state is given by the real part of the corresponding Floquet exponents [72], which we denote by μ. The MSF provides the largest growth rate μof a perturbation acting along each eigenmode Φ (α) as a function of its associated eigenvalue Λ α , ma¼MSFðLaÞ:ð1Þ This relation is analogous to dispersion relations in spatially extended systems, where the growth rate of a perturbation is given as a function of the perturbation wavenumber [55]. In the Jansen model, the MSF needs to be computed numerically (see Methods section for a detailed explanation of its derivation and numerical computation). Fig 2(b) shows the MSFs of the homogeneous limit-cycle solution for different pvalues. Each circle indicates the real part of the largest Floquet exponent corresponding to a Fig 2. Transverse instabilities of oscillatory states. (a) First 5 eigenvectors of ~ Win a spatial representation of the brain (superior view). Each network node has been colored according to their contribution to the corresponding eigenvector ΦðaÞ¼ ð�ðaÞ 1;...; �ðaÞ NÞ. (b) Master Stability Function of system (10) showing the dependence of the largest Floquet exponent μwith respect to the structural connectivity eigenvalues Λ α for three different values of the input p.�= 50 in all cases. Circles correspond to the eigenvalues of ~ W, whereas black curves are obtained by continuously tuning Λ. (c,d) Results from numerical simulations of Eq (10) with (c) p= 280 and �= 50 and (d) p= 210 and �= 50. (e) Complete bifurcation diagram of the homogeneous states. Colors of lines and regions as in Fig 1(e) and (f), with the region of transverse instability (i.e., μ 2 >0) shaded in pink and delimited with a black curve. Black circles correspond to numerical simulations in which hσi>10 −5 . Simulations initialized close to a homogeneous state. (f,g) Bifurcation diagram obtained from direct simulations of the system with initial conditions close to homogeneous (f) and random (g). Continuous curves as in Fig 1(e) and (f). Regions colored according the dynamical classification given by the two largest Lyapunov exponents (see Methods). https://doi.org/10.1371/journal.pcbi.1010781.g002 PLOS COMPUTATIONAL BIOLOGY Transverse instabilities in large-scale brain networks PLOS Computational Biology | https://doi.org/10.1371/journal.pcbi.1010781 April 12, 2023 7 / 34 perturbation applied along the αth eigenvector computed according to the MSF. Since we are considering a stable limit-cycle solution given by the analysis of system (18), the largest growth rate corresponding to the uniform eigenvector (i.e., Λ 1 = 1) is μ 1 = 0. The other exponents might be positive or negative depending on both the system parameters and Λ α . For instance, for p= 280 (green circles) the dispersion relation is negative for all the structural eigenvalues −1<Λ α <1. Therefore, in this case we expect small inhomogeneous perturbations of the homogeneous state to decay exponentially. In contrast, for p= 210 (blue circles) some of the connectivity matrix eigenvectors have positive growth rates μ α , thus small perturbations should give rise to heterogeneous patterns. Indeed, Fig 2(c) and 2(d) show the results of integrating numerically Eq (10) for p= 280 and 210, respectively. The initial conditions correspond to a uniform state plus a small random perturbation. In agreement with the results given by the stability analysis, after a short transient (not shown), the dynamics for p= 280 falls back to the homogeneous state, whereas a heterogeneous spatiotemporal pattern arises for p= 210. The MSF dispersion relations represented in Fig 2(b) show that the positive growth rates appear first through the second largest structural-connectivity eigenvalue, Λ 2 . We thus extensively analyze the loci of unstable directions by checking whether μ 2 has a positive real part for the entire region of existence of oscillatory activity. Fig 2(e) displays the region where μ 2 >0 (purple) superimposed on the bifurcation diagram of the homogeneous manifold, calculated in the previous section and shown originally in Fig 1(f). Remarkably, the region of transverse instability occupies a large portion of the parameter space. Moreover, the corresponding values of �cover a range comparable to the other coupling parameters of the system, C 1 ,. . .,C 4 2[0, 135]. In order to validate the emergence of heterogeneous dynamics in the system we perform numerical simulations of system (10) for different values of pand �, starting from initial conditions close to a homogeneous state. At each time step we compute the standard deviation across network nodes, sðtÞ≔1 NX N i¼1ðviðtÞ � vðtÞÞ2 " #1 2 where � vðtÞ ¼ 1 NX N i¼1 viðtÞ:ð2Þ This quantity vanishes in the homogeneous states, whereas it is positive in heterogeneous states. Black circles in Fig 2(b) indicate the parameter values for which hσ(t)i>10 −5 in the simulations, showing a complete overlap with the results coming from the linear stability analysis (purple region in the figure). Finally, we characterize the type of dynamical states arising in the region of transverse instability, by computing the two largest Lyapunov exponents (LE, see Methods), λ 1 and λ 2 , in independent numerical simulations with varying pand �. Using this tool we can classify the attractors of the system depending on the sign of the two exponents (see Methods). Fig 2(f) and 2(g) show the resulting numerical bifurcation diagrams corresponding to initial conditions close to a homogeneous state (Fig 2(f)) or entirely random initial conditions (Fig 2(g)). In both cases the region of transverse instability can be roughly divided in two parts: one dominated by periodic heterogeneous oscillations (large �, blue), and one displaying chaotic dynamics (small �, red) with at least two unstable directions. Additionally, simulations initialized at random (Fig 2(g)) show that the chaotic regime extends much beyond the transverse instability region, thus uncovering a coexistence region between spatiotemporal chaos and homogeneous states. In the next sections we investigate the properties of these different dynamical regimes in detail. PLOS COMPUTATIONAL BIOLOGY Transverse instabilities in large-scale brain networks PLOS Computational Biology | https://doi.org/10.1371/journal.pcbi.1010781 April 12, 2023 8 / 34 Periodic travelling waves The simplest instance of heterogeneous dynamics in system (10) is a periodic regime (blue region in Fig 2(f) and 2(g)). In this regime, the dynamics of each brain region is purely periodic, but with a non-zero phase difference between regions, i.e., there is phase-locking between nodes. Here we reveal that a multiple of such states might exist for a single choice of parameter values. Fig 3(a) and 3(b) show the difference between each node’s phase ϕ j and the collective phase C(see Methods) for two simulations with fixed p= 230 and �= 50 and different initial conditions. In both cases the initial conditions are set close to the homogeneous (unstable) limitcycle, plus a small random perturbation. The first case (simulation A, Fig 3(a)) reveals a wave pattern that travels from the left to the right hemisphere, whereas the second (simulation B, Fig 3(b)) displays a travelling wave that goes from the parietal to the frontal region, (see also S1 and S2 Movies). Following [30] we obtain the direction and speed of the wave propagation by means of constrained natural element differentiation (see Methods, and also [30,73]). The resulting vector field (Fig 3(c) and 3(d)) provides better indication on the specific type of pattern shown by each simulation, together with the wave propagation speed of each region. These propagation speeds are heterogeneous in space, as expected from the irregular brain connectome, with values ranging between 1 and 8m/s (see Fig 3(e)), in close agreement to results from non-invasive recordings, but one order of magnitude larger than those observed in invasive recordings [74]. This disagreement might be due to the fact that ours is a model of inter-cortical wave Fig 3. Multistability between different types of periodic travelling waves. (a,b) Difference between each node’s phase ϕ j and the collective phase C corresponding to simulations with p= 230 and �= 50, in the spatial representation of the brain network (from left to right: frontolateral, superior, and frontal views). Panels (a) and (b) correspond to two different initial conditions. (c,d) Propagation direction vectors ζ j /z j corresponding to the phase patterns of (a) and (b). Color indicates the propagation speed ν j . (e,f) Swarm plot of propagation speed z j (e) and local polarization a j (f). The circles correspond to individual brain regions. Black squares show the average over the entire network, and error bars indicate standard deviation. (g) Contribution of each structural connectivity eigenmode to the growth rate of the perturbations in simulations A and B (see Eq (29) in Methods). https://doi.org/10.1371/journal.pcbi.1010781.g003 PLOS COMPUTATIONAL BIOLOGY Transverse instabilities in large-scale brain networks PLOS Computational Biology | https://doi.org/10.1371/journal.pcbi.1010781 April 12, 2023 9 / 34 both the periodic and chaotic dynamics. Second, the amplitude correlates mildly with the node strength at the chaotic region (red squares in Fig 6(f)), but appears independent from s i close to the onset of oscillatory activity (blue circles in Fig 6(f)) Finally, Fig 6(g) shows a negative correlation between the phase difference of each brain region and their corresponding node strength: Regions with lower input tend to oscillate first, and the hubs tend to be the last. In the periodic regime (e.g., p= 320) this relation translates to a travelling wave in which oscillations spread from outer to the inner region of the brain (see S4 Movie). In the chaotic regime (p= 170) the role of the node strength on the dynamics is still prominent (blue squares in Fig 6(g)), although much blurred by the irregular behavior of the system, as shown by the large error bars (see also S5 Movie). Overall, the non-normalized network topology adds a new layer of complexity in the model, as homogeneous states no longer exist. Nonetheless, many features of the dynamics can be traced back directly to the distribution of node strengths s i . Moreover, the agreement between diagrams in Fig 6(a) and 6(b) and those of Fig 2(f) and 2(g) indicate that the onset of chaotic dynamics is mostly retained in the row-normalized network topology ~ W. Spatiotemporal chaos in a large-scale brain model with next generation NMMs We have seen so far that irregular spatiotemporal dynamics arise in networks of coupled NMMs from transverse instabilities of oscillatory states. We expect this to be a general mechanism in brain networks. In order to illustrate this ubiquity, we now analyze a large-scale brain model composed of coupled next generation NMM (NG-NMM). These models are derived from an exact mean-field theory for quadratic integrate-and-fire neurons, and therefore, the corresponding firing rate equations can be traced back to the dynamics of single neurons [17]. Here we consider a pyramidal-interneuronal network gamma (PING) setup, in which the local dynamics within each brain region produces gamma activity through the interplay between excitatory and inhibitory neurons (see Methods, and also [83]). In order to obtain a bifurcation diagram, we proceed as we have done above for the Jansen system: first we investigate the homogeneous manifold of the system, and then apply the MSF formalism. The manifold of homogeneous trajectories corresponds again to a self-coupled version of the model (Eq (38) in Methods). Fig 7(a) shows the bifurcation diagram corresponding to homogeneous trajectories of the system, using as control parameters the external-to-excitatory baseline input η e and the coupling strength �. For weak coupling (��1) and low η e the system remains in a single low-activity fixed point (white area in the plot) corresponding to a state of asynchronous dynamics. Upon increasing the constant external input η e , such steady state gives raise to fast oscillatory activity (blue-shaded area) through a supercritical Hopf bifurcation (red continuous line). For large values of η e the fixed point recovers stability through a subcritical Hopf (dashed red curve), giving raise to a bistable state between gamma activity and asynchrony (orange-shaded area). For even larger values of η e gamma activity finally vanishes through a saddle-node of limit cycles (outside figure range). As it happened with the Jansen system, increasing �leads to a change on the onset of oscillatory activity. Following a typical scenario in oscillatory systems (see, e.g., [84–88]), in a tiny region of the parameter space (see panel (b)) three codimension-2 bifurcations coexist: a Bogdanov-Takens (BT), a saddle-node separatrix loop (SNSL), and a cusp of saddle nodes. The Hopf bifurcation vanishes at a Bogdanov-Takens (BT) point, and a saddle-node separatrix loop (SNSL) gives raise to a SNIC branch (dark blue line). Therefore, for most values of �, gamma activity arises through a infinite-period (SNIC) bifurcation. Also, an increase of the coupling causes the region of bistability between the oscillatory states and the fixed-point to PLOS COMPUTATIONAL BIOLOGY Transverse instabilities in large-scale brain networks PLOS Computational Biology | https://doi.org/10.1371/journal.pcbi.1010781 April 12, 2023 16 / 34 increase (orange-shaded area). Overall, homogeneous oscillatory activity –with or without bistability– dominates the bifurcation diagram. An example of such homogeneous gamma activity is displayed in Fig 7(c), corresponding to η e = 5 and �= 30 (black triangle in Fig 7(a)). By applying the MSF formalism on these oscillatory states, we unveil a region of transverse instability (pink-shaded region in Fig 7(a)), which emerges for large values of �. As a result, simulations initialized close to the homogeneous state within this parameter region exhibit irregular spatiotemporal patterns, as the one depicted in Fig 7(d) (η e = 5, �= 40 i.e, black square in Fig 7(a)). Notice that the frequency of the oscillations becomes almost twice as that of the (unstable) homogeneous state, which in this case is approximately 46Hz. This is a differentiating feature with respect to the Jansen model, in which heterogeneous activity is faster than that of the underlying homogeneous state, but always close to the typical values displayed by the model. Overall, the large-scale brain model with NG-NMMs displays a mechanism towards the onset of irregular states analogous to the Jansen case, despite the different nature of the two models [89]. For instance, we did not include synaptic dynamics in Eq (36), i.e., neurons receive instantaneous delta-like pulses. The onset of spatiotemporal chaos shown here supports thus the generality of transverse instabilitites in the oscillatory dynamics of coupled NMMs. Discussion The idea that systems composed of simple deterministic subunits can give rise to complex spatiotemporal behavior traces back to the seminal work of Alan Turing [51]. That work showed that a homogeneous equilibrium in a system of diffusively coupled units might lose stability due to the coupling between neighbouring sites, thereby producing spatial patterns. Decades of research have extended this simple yet powerful framework to cover a wide range of possibilities, including instabilities arising from uniform oscillatory states [52,53,90] and pattern formation in complex networks [54]. Both of these extensions are embodied in the field of collective synchronization, in which the Master Stability Function provides a proper formalism to analyze the stability of homogeneous oscillatory states [58–62]. In brain networks, the stability of synchronized states has been analyzed only in simplified network topologies [48,50] or by means of phase-reduction approximations that only apply for weak coupling [33,35]. In this paper we have studied a general scenario, that reveals transverse instabilities of homogeneous states as an ubiquitous mechanism for the onset of travelling waves and high-dimensional chaos in large-scale brain models. In computational neuroscience, the spontaneous emergence of patterns through instabilities of a uniform state has been emphasized mostly in the context of neural fields [91–95]. Neural fields are models defined in a continuous spatial support, where the synaptic coupling is a smooth function of the distance between regions. Hence, these type of models cannot capture the fine macroscopic organization of the brain connectome represented by complex networks, and are thus more adequate to model local intra-cortical dynamics [96,97]. Nonetheless, the assessment of transverse instabilities is general enough to apply to both continuous spatial support with simplified interaction rules and neural mass models interacting through complex networks. The main difference between the two cases lays on the decomposition of the perturbation vector. In neural fields and regular network topologies [50], as in the Turing framework, stability analysis of homogeneous states is attained by decomposing a spatial perturbation in Fourier space. Instead, in complex networks composed of coupled NMM, the MSF requires the diagonalization of the structural connectivity matrix. Interestingly, some studies show that spectral analysis of whole-brain networks enables the characterization of PLOS COMPUTATIONAL BIOLOGY Transverse instabilities in large-scale brain networks PLOS Computational Biology | https://doi.org/10.1371/journal.pcbi.1010781 April 12, 2023 17 / 34 functional and resting brain activity [98–100]. Our work provides thus a mathematical framework to further explore the relation between structural connectivity (in terms e.g. of graph spectra) and functional activity data. An important difference between large-scale neural dynamics and classical pattern formation lies in the nature of coupling. Classical pattern forming systems usually involve diffusive coupling, for which the homogeneous manifold of the system coincides with the dynamics of the uncoupled units [55]. In contrast, the long-range brain connections represented by structural connectivity matrices correspond to myelinated nerve tracts across brain regions. Hence, the interactions between nodes are mediated by chemical synapses driven by the firing rate of pre-synaptic regions. As a result, the homogeneous manifold of the system is given by a selfcoupled version of the single-node model, and therefore the coupling strength modifies the uniform states of the system in a non-trivial manner. We have shown, for instance, that in both Jansen and next generation NMMs, there is a change on the onset of synchronized activity from a Hopf to a SNIC bifurcation upon increasing the coupling parameter �. Additionally, coupling eliminates all forms of bistability in the homogeneous states of the Jansen model, whereas it enlarges the region of coexistence between homogeneously attracting limit-cycles and steady states in the NG-NMM. Most of the theory of pattern formation is based on studying instabilities from fixed points. Although we also analyzed this case, we did not find instabilities arising from stationary states. Instead, the linear stability of oscillatory solutions revealed the ubiquity of transverse instabilities as a route for the emergence of spatiotemporal chaos in large-scale brain models. This mechanism explains the onset of deterministic chaos outlined in previous works by means of direct numerical simulations [44,49]. Moreover, the numerically computed Master Stability Functions (Fig 2(b)) show that the onset of unstable modes occurs through eigenvalues arbitrarily close to the zero eigenmode. This scenario is very close to that of the Benjamin-Feir instability in the Ginzburg-Landau system, which was studied by Kuramoto as a main route to turbulence in chemical reaction systems [52,53,60]. The emergence of high-dimensional chaos in large-scale brain networks is in line with recent studies that highlight the prominence of turbulent dynamics as a signature of healthy brain activity [8,71,101–103]. These complex fluctuations are observed in blood-oxygenlevel-dependent (BOLD) signals obtained from fMRI scans, which are characterized by slow (below 1Hz) fluctuations. These are in the same time scale as the chaotic slow modulations of alpha rhythms that we observed in our brain network composed of Jansen NMMs. Remarkably, experimental studies have also found a negative correlation between the amplitude of alpha rhythms and BOLD activity [104–107]. Thus, we conjecture that the chaotic dynamics generated by transverse instabilities might capture the spatiotemporal fluctuations of BOLD signals observed in recordings. Another important feature of the chaotic dynamics in the Jansen model is the coexistence of two frequency ranges, theta (�5Hz) and alpha (�10Hz), in some regions of the parameter space (see Fig 4(c)–4(e)). These alternating bursts of activity between the two rhythms are consistent with the multifrequency and transient behavior of oscillations in EEG recordings. Large-scale brain models aim to trace the basic principles behind neural activity. Therefore, one might need to account for other sources of complexity not included in this work. For instance, despite being a common choice in the literature [30,33,35,44,45,56,57], in some cases it would be preferable to use non-normalized topologies. We have shown that regardless of this simplification, complex spatiotemporal patterns and high-dimensional chaos exist already in the row-normalized system. Moreover, the numerically-derived bifurcation diagram of the non-normalized system shown in Fig 6(a) and 6(b) qualitatively matches the dynamical landscape uncovered by the analysis of the homogeneous states of the simplified model. PLOS COMPUTATIONAL BIOLOGY Transverse instabilities in large-scale brain networks PLOS Computational Biology | https://doi.org/10.1371/journal.pcbi.1010781 April 12, 2023 18 / 34 Additionally, here we have considered long-range excitatory coupling targeting only pyramidal neurons. The MSF formalism can also be applied in the case where long-range excitation targets both excitatory and inhibitory units. A simple exploration of this setup using Jansen NMMs indicated that no transverse instabilities emerged in this situation (results not shown). This is in agreement with [50], which shows that self-sustained traveling waves emerge in a ring of Wilson-Cowan units when cross-excitation only affects excitatory populations but vanish if interneurons are also targeted. Finally, the assumption of homogeneous dynamics across brain regions might be challenged by previous findings that indicate hierarchical heterogeneities in synaptic strengths and time scales [108,109]. Nonetheless, a deeper understanding of transverse instabilities in simplified systems provides a solid ground to analyze the emergence and propagation of spatiotemporal neural activity in comprehensive models including such heterogeneities. To this end, future work should consider quantitative comparisons between the dynamics emerging from transverse instabilities and dynamical data from EEG or fMRI recordings. This might help to validate some of the aforementioned modeling choices as well as providing relevant estimations for �, usually given by fitting the model to available data (see, e.g. [8,29]). Materials and methods Structural connectivity data The structural connectivity of our NMM network has been obtained from diffusion tensor imaging (DTI) data of 16 healthy subjects, collected and analyzed in a previous study [65]. The different brain regions were defined using the AAL90 parcellation [64]. A 90 ×90 structural connectivity matrix Wwas obtained averaging the connectome of the 16 subjects. For more details about data collection and prepossessing we refer the readers to [65]. For details on the further normalization used in most of this paper, see Section entitled “Normalized connectivity” below. The Jansen NMM model Single region. Jansen’s NMM (also known as Jansen-Rit model) describes the dynamical activity of a cortical column of neurons [15,16,63]. Following principles derived from earlier works [13,110], the model assumes populations of pyramidal neurons (PN) and inhibitory interneurons (IN), Recurrent connections are only present in the PN population, whereas inhibitory neurons solely receive inputs from pyramidal neurons For simplicity the model assumes that recurrent connections within the PN population are mediated through neurons that do not receive direct input from inhibitory neurons, and can therefore be interpreted as a third independent population, which we call recurrent pyramidal neurons (rPN). Finally, pyramidal neurons also receive excitatory external stimuli from other brain regions, modelled through a firing rate variable I(t). In the Jansen model the excitatory and inhibitory post-synaptic potentials (PSP) are given by heðtÞ ¼ AateatHðtÞ hiðtÞ ¼ BbtebtHðtÞð5Þ where His the Heaviside step function, Aand Bare the PSP amplitudes, and aand bquantify the synaptic time scales. As a result, a neural population receiving an excitatory firing rate of r(t) generates an excitatory post-synaptic potential (ePSP) described by the second-order PLOS COMPUTATIONAL BIOLOGY Transverse instabilities in large-scale brain networks PLOS Computational Biology | https://doi.org/10.1371/journal.pcbi.1010781 April 12, 2023 19 / 34 differential equation € yðtÞ ¼ AarðtÞ 2a_ yðtÞ a2yðtÞ:ð6Þ Analogously, a inhibitory firing rate generates a inhibitory post-synaptic potential (iPSP) according to € yðtÞ ¼ BbrðtÞ 2b_ yðtÞ b2yðtÞ:ð7Þ On the other hand, a population with mean membrane potential v(t) generates a firing rate according to the sigmoid function SigmðvÞ≔2e0 1þerðv0vÞ:ð8Þ The model assumes that the self-dynamics of neurons within a population can be averaged out. Thus, the entire system evolves as determined by the evolution of the average PSPs of the different populations given by Eqs (6) and (7). As a result, the interactions between the three populations (PN, IN, and rPN) lead to the 6-dimensional Jansen model: _ y0ðtÞ ¼ y3ðtÞ _ y1ðtÞ ¼ y4ðtÞ _ y2ðtÞ ¼ y5ðtÞ _ y3ðtÞ ¼ Aa Sigm½y1ðtÞ y2ðtÞ� 2ay3ðtÞ a2y0ðtÞ _ y4ðtÞ ¼ AafIðtÞþC2Sigm½C1y0ðtÞ�g 2ay4ðtÞ a2y1ðtÞ _ y5ðtÞ ¼ BbC4Sigm½C3y0ðtÞ� 2by5ðtÞ b2y2ðtÞ; ð9Þ where y 0 accounts for the ePSP generated by the PNs, y 1 is the sum of the ePSP generated by the rPNs and the external excitatory inputs, and y 2 is the iPSP generated by the INs. Finally, the mean membrane potential of the pyramidal population is given by v≔y 1 −y 2 , which we use as the main observable throughout this paper. Network model. We consider a network composed of Nnodes, each representing a brain region. The dynamics of each node follows the Jansen model for a cortical column. Therefore, the 6 equations governing the dynamics of node iread [35] _ y0;iðtÞ ¼ y3;iðtÞ _ y1;iðtÞ ¼ y4;iðtÞ _ y2;iðtÞ ¼ y5;iðtÞ _ y3;iðtÞ ¼ Aa Sigm½y1;iðtÞ y2;iðtÞ� 2ay3;iðtÞ a2y0;iðtÞ _ y4;iðtÞ ¼ AafIiðtÞþC2Sigm½C1y0;iðtÞ�g 2ay4;iðtÞ a2y1;iðtÞ _ y5;iðtÞ ¼ BbC4Sigm½C3y0;iðtÞ� 2by5;iðtÞ b2y2;iðtÞ: ð10Þ The quantity I i accounts for the incoming signals from the rest of the network or other layers not represented in the model. In the brain model we consider that the different regions are coupled only through excitation, thus inhibition acts only locally. Also, all regions receive an external input from subcortical regions not represented in our model, in the form of a PLOS COMPUTATIONAL BIOLOGY Transverse instabilities in large-scale brain networks PLOS Computational Biology | https://doi.org/10.1371/journal.pcbi.1010781 April 12, 2023 20 / 34 common constant firing rate p. Altogether, I i (t) takes the form of the sum of two independent terms: IiðtÞ ¼ pþ�X N j¼1 ~ Wij Sigm½y1;jðtÞ y2;jðtÞ� ;ð11Þ where �is the coupling strength, and ~ W¼ ð~ wijÞis the row-normalized structural connectivity matrix (see next section). We study the system under the influence of the external driving firing rate pand the coupling strength �, leaving all the other parameters fixed, as defined in Table 1. As mentioned above, from the six variables that characterize the dynamics of each brain region, we monitor the mean membrane potential of the pyramidal neurons, v i ≔y 1,i −y 2,i . Normalized connectivity. Given a N×Nstructural connectivity (SC) matrix W, a detailed mathematical analysis of system (10) is generally unfeasible. In order to allow for a semi-analytical treatment, we consider a normalized version of the topology. Such normalization is usually employed in large-scale brain models [30,33,35,44,45,56,57]. Let D= (d ij ) be a diagonal N×Nmatrix whose non-zero entries are the sum of incoming connections to each node: dij ¼ sii¼j 0i6¼ j (ð12Þ where si¼X N j¼1 wij ð13Þ is the in-strength of node i. Then we consider the row-normalized SC matrix ~ W¼ ð~ wijÞ≔D1Wð14Þ whose elements are ~ wij ≔wij=si, i.e., ~ Wis obtained dividing each row of Wby its sum. Therefore, the sum of the elements on each row equals unity i.e., X N j¼1 ~ wij ¼1:ð15Þ Table 1. Parameters of the Jansen system. Parameter Meaning Value AMaximal amplitude of excitatory post-synaptic potentials 3.25mV BMaximal amplitude of inhibitory post-synaptic potentials 22mV aCharacteristic decay time for ePSP 100s −1 bCharacteristic decay time for iPSP 50s −1 C 1 ,C 2 ,C 3 ,C 4 Synaptic strength (average number of synapses) between populations 135, 108, 33.75, 33.75 e 0 Half of the maximum firing rate 2.5Hz v 0 Potential where half the maximum firing rate is achieved 6mV rNeuronal excitability 0.56mV −1 ~wij Connectivity weights from data pConstant baseline firing rate to pyramidal neurons not fixed (Hz) �Coupling strength not fixed https://doi.org/10.1371/journal.pcbi.1010781.t001 PLOS COMPUTATIONAL BIOLOGY Transverse instabilities in large-scale brain networks PLOS Computational Biology | https://doi.org/10.1371/journal.pcbi.1010781 April 12, 2023 21 / 34 Matrices with this normalization are sometimes called right stochastic matrices [111]. An important property of right stochastic matrices that we use in our analysis is that their largest eigenvalue is exactly Λ 1 = 1, which corresponds to a uniform eigenvector ϕ (1) ≔(1, . . ., 1) T . By the by the Gershgorin circle theorem [112], all other eigenvalues are bounded within the unit circle. Finally, we remark that since the original structural connectome W, is symmetric, so is the matrix Z≔D1 2WD1 2. Hence Zis diagonalizable and has real eigenvalues. Using that W¼D~ Wwe obtain Z≔D1 2~ WD1 2:ð16Þ Therefore, Zand ~ Ware similar (in the mathematical sense), i.e., they share the same eigenvalues. The fact that the eigenvalues o ~ Ware real (due to the symmetry of W) is not strictly necessary to carry our analysis based on the MSF, but it simplifies it. Homogeneous states A row-normalized connectivity matrix ensures that all the different units in the network receive the same amount of input, although distributed differently across the different nodes. Using such type of connectivities, it is always possible to find homogeneous or uniform solutions -either stationary or time dependentin which all nodes of the network behave identically. Let y 0 ,. . .,y 5 be the variables that characterize the dynamics of each brain region in a homogeneous state. Imposing then y m,i (t) = y m (t) for m= 0, . . ., 5 and i= 1, . . .,None finds that the incoming input for each brain region in (10) reads IiðtÞ ¼ pþ�X N j¼1 ~ wij Sigm½y1;jðtÞ y2;jðtÞ� ¼ pþ�Sigm½y1ðtÞ y2ðtÞ� ;ð17Þ thus it does not depend on the node index ianymore. Replacing this expression of the input in the Jansen model (10) one finds that, in a homogeneous state, the equations for the evolution of each network node read _ y0ðtÞ ¼ y3ðtÞ _ y1ðtÞ ¼ y4ðtÞ _ y2ðtÞ ¼ y5ðtÞ _ y3ðtÞ ¼ Aa Sigm½y1ðtÞ y2ðtÞ� 2ay3ðtÞ a2y0ðtÞ _ y4ðtÞ ¼ Aafpþ�Sigm½y1ðtÞ y2ðtÞ�þC2Sigm½C1y0ðtÞ�g 2ay4ðtÞ a2y1ðtÞ _ y5ðtÞ ¼ BbC4Sigm½C3y0ðtÞ� 2by5ðtÞ b2y2ðtÞ: ð18Þ Therefore, this low-dimensional system determines all homogeneous states of the coupled system (10). Moreover, this system also retains the stability of such homogeneous states subject to uniform perturbations, i.e., perturbations that act identically at each brain region, and therefore do not change the homogeneous character of the trajectories. In other words, Eq (18) define an invariant manifold of Eq (10). Nonetheless, stable states in the homogeneous invariant manifold might still be unstable to heterogeneous perturbations, i.e., perturbations transverse to the manifold. PLOS COMPUTATIONAL BIOLOGY Transverse instabilities in large-scale brain networks PLOS Computational Biology | https://doi.org/10.1371/journal.pcbi.1010781 April 12, 2023 22 / 34 Transverse stability The following stability analysis is common in the study of dynamical systems in complex networks, specially (but not only), in the context of diffusive coupling [54,68,69]. The same technique can be applied to both, homogeneous fixed points and homogeneous limit-cycles. In the later case it is known as the Master Stability Function (MSF) [58] (see also [59,61] for introductory reviews). In this section we explain such stability analysis on the Jansen model, but it can be easily extended to other systems. We use bold symbols for Nand 6Ndimensional vectors and matrices, and regular symbols for 6-dimensional vectors, 6-dimensional matrices, and scalar quantities. All scalar quantities except the time tand the coupling strength �have a subscript. Let yð0Þ≔ðyð0Þ 0;. . . ;yð0Þ 5ÞTbe a solution of (18). Then, the 6N-dimensional vector yð0Þ¼ ðyð0Þ 0;...;yð0Þ 0;. . . z}|{ N ;yð0Þ 0;. . . ;yð0Þ 5ÞTð19Þ is a homogeneous solution of (10), either stationary or periodic. Let us consider an arbitrary small perturbation δy¼ ðdy0;1;...;dy5;1;. . . ;dy0;N;. . . ;dy5;NÞTð20Þ where each component δy k,j is the perturbation acting on the variable kof node jin the system (10). Expanding the velocity fields of (10) up to first order around y (0) one obtains the linear evolution of the perturbation vector as d dt δyðtÞ ¼ Jðyð0ÞÞδyðtÞð21Þ where Jis the full 6N×6NJacobian of system (10), evaluated at y (0) . Notice that Jcan be written as Jðyð0ÞÞ ¼ Jðyð0ÞÞ�INþ�Kðyð0ÞÞ� ~ Wð22Þ where �denotes the Kronecker product, Jis the 6 ×6 Jacobian of the uncoupled Jansen Eq (9),I N is the N×Nidentity matrix, and K= (k ij ) is a 6 ×6 matrix defined by kij ¼ Aa Sigm0ðy1y2Þi¼5;j¼2 Aa Sigm0ðy1y2Þi¼5;j¼3 0otherwise : 8 > > > < > > > :ð23Þ One could, in principle, evaluate numerically the eigenvalues and eigenvectors of this Jacobian in order to obtain the stability properties of the system. Nonetheless, there is a simpler and more informative approach based on expressing the perturbation vector δyin an adequate coordinate system. Let ΦðaÞ¼ ð�ðaÞ 1;...; �ðaÞ NÞTbe a normalized eigenvector of ~ Wassociated with the eigenvalue Λ α for α= 1, . . .,N, so that ~ WΦðaÞ¼LaΦðaÞ:ð24Þ The set of eigenvectors {Φ (1) ,. . .,Φ (N) } constitute a basis of the vector space RN, thus we can express the perturbation δyof the homogeneous solution as a linear combination of such PLOS COMPUTATIONAL BIOLOGY Transverse instabilities in large-scale brain networks PLOS Computational Biology | https://doi.org/10.1371/journal.pcbi.1010781 April 12, 2023 23 / 34 vectors. In this new basis, the perturbation acting on variable kof the jth node reads dyk;jðtÞ ¼ X N a¼1 uðaÞ kðtÞ�ðaÞ jð25Þ where uðaÞðtÞ ¼ ðuðaÞ 0ðtÞ;...;uðaÞ 5ðtÞÞTare the coordinates of the perturbation vector expressed in the new basis. This can also be expressed in vectorial form using the Kronecker operator � as δyðtÞ ¼ X N a¼1 uðaÞðtÞ�ΦðaÞ:ð26Þ Now it is necessary to perform some calculations using the tensor product �. For simplicity we drop some dependences: J=J(y (0) ), K=K(y (0) ), and u (α) =u (α) (t). Applying the change of coordinates on Eq (21) and using Eq (22), we have that d dt δy¼d dt X N a¼1 uðaÞ�ΦðaÞ ¼X N a¼1 d dt uðaÞ � ��ΦðaÞ ¼ ðJ�INþ�K�~ WÞX N a¼1 uðaÞ�ΦðaÞ ¼ ðJ�INÞX N a¼1 uðaÞ�ΦðaÞ !þ�ðK�~ WÞX N a¼1 uðaÞ�ΦðaÞ ! ¼X N a¼1ðJ�INÞðuðaÞ�ΦðaÞÞþ�X N a¼1ðK�~ WÞðuðaÞ�ΦðaÞÞ ¼X N a¼1ðJuðaÞÞ�ðINΦðaÞÞþ�X N a¼1ðKuðaÞÞ�ð ~ WFðaÞÞ ¼X N a¼1ðJuðaÞÞ�ΦðaÞþ�X N a¼1ðKuðaÞÞ�ðLaΦðaÞÞ ¼X N a¼1ðJuðaÞÞ�ΦðaÞþ�X N a¼1ðLaKuðaÞÞ�ΦðaÞ ¼X N a¼1ðJþ�LaKÞuðaÞ�ΦðaÞ; ð27Þ where we have used the diagonalization of the connectivity matrix Eq (24). Now, making use of the linear independence of the eigenvectors fΦðaÞgN a¼1one obtains that the evolution of u (α) (t) becomes independent for each α= 1, . . .,Nthrough the relation _ uðaÞ¼ ðJþ�LaKÞuðaÞ;ð28Þ where Jðyð0Þ;LÞ≔Jþ�LaKis a family of 6 ×6 Jacobians that depend on the homogeneous state of the system y (0) and on the structural connectivity eigenvalues Λ α . If y (0) is a fixed point, PLOS COMPUTATIONAL BIOLOGY Transverse instabilities in large-scale brain networks PLOS Computational Biology | https://doi.org/10.1371/journal.pcbi.1010781 April 12, 2023 24 / 34 the stability of the homogeneous solution simplifies to studying the eigenvalues (and eigenvectors) of the Jacobian J, which is a function of the connectivity matrix eigenvalues Λ α . If y (0) = y (0) (t) is a periodic solution, then Jis a periodic matrix and Floquet theory applies [72]. In this context, the growth rate of a perturbation acting at the limit-cycle solution is determined by the real part of the Floquet exponents corresponding to J, which must be determined numerically. Practically speaking, for each value of Λ α , we compute the real part of the largest Floquet exponent of Jas the largest Lyapunov exponent of the linear system Eq (28). The code is available in github (www.github.com/pclus/transverse-instabilities). Finally, in order to analyze the contribution of each structural eigenmode αon the growth rate of a specific perturbation vector δywe compute the quantity ca¼makuðaÞk:ð29Þ Lyapunov exponents Lyapunov exponents (LE) provide the growth rate of small perturbations acting on a timeevolving trajectory of a dynamical system [77,113]. To compute these quantities in practice we employ a dynamical algorithm based on QR-decompositions, using Householder reflections [77]. The algorithm is embedded in the available software (www.github.com/pclus/transverseinstabilities). There are as many Lyapunov exponents as system dimensions, and they are usually sorted from largest to smallest: λ 1 �λ 2 �λ 3 �. . . �λ 6N . Using the two largest Lyapunov exponents λ 1 and λ 2 , we can classify the system state in 5 types: • Fixed points, corresponding to both LE being negative (λ 1 ,λ 2 <0). • Periodic dynamics, identified by a zero LE and the rest being negative (λ 1 = 0, λ 2 <0). • Quasiperiodic dynamics, i.e., a regime in which the system evolves with at least two incommensurate characteristic frequencies. In this case both LE are zero (λ 1 = 0, λ 2 = 0). • Chaotic dynamics, with a single positive LE (λ 1 >0, λ 2 �0). • Hyperchaos, with more than one positive LE (λ 1 ,λ 2 >0). Since the Lyapunov exponents are computed numerically, we need to impose a threshold to discern between zero and non-zero values. We found that a value of |λ k |<10 −4 was a reasonable cut-off. Another useful application of Lyapunov exponents is to compute the dimensionality of the attractor, which is fractal in chaotic states. The Kaplan-Yorke formula [80] provides an approximation of the fractal dimension of an attractor as DKY ¼jþPj i¼1li jljþ1j;ð30Þ where jis the LE for which X j i¼1 li�0and X jþ1 i¼1 li<0:ð31Þ PLOS COMPUTATIONAL BIOLOGY Transverse instabilities in large-scale brain networks PLOS Computational Biology | https://doi.org/10.1371/journal.pcbi.1010781 April 12, 2023 25 / 34 60. Nakao H. Complex Ginzburg-Landau equation on networks and its non-uniform dynamics. The European Physical Journal Special Topics. 2014; 223(12):2411–2421. https://doi.org/10.1140/epjst/e201402220-1 61. Porter M, Gleeson J. Dynamical Systems on Networks: A Tutorial. Frontiers in Applied Dynamical Systems: Reviews and Tutorials. Springer International Publishing; 2016. Available from: https://books. google.es/books?id=uzDuCwAAQBAJ. 62. Ashwin P, Coombes S, Nicks R. Mathematical Frameworks for Oscillatory Network Dynamics in Neuroscience. The Journal of Mathematical Neuroscience. 2016; 6(1):2. https://doi.org/10.1186/s13408015-0033-6 PMID: 26739133 63. Grimbert F, Faugeras O. Analysis of Jansen’s model of a single cortical column. INRIA; 2006. RR5597. Available from: https://hal.inria.fr/inria-00070410. 64. Tzourio-Mazoyer N, Landeau B, Papathanassiou D, Crivello F, Etard O, Delcroix N, et al. Automated Anatomical Labeling of Activations in SPM Using a Macroscopic Anatomical Parcellation of the MNI MRI Single-Subject Brain. NeuroImage. 2002; 15(1):273–289. https://doi.org/10.1006/nimg.2001. 0978 PMID: 11771995 65. Deco G, Cabral J, Woolrich MW, Stevner ABA, van Hartevelt TJ, Kringelbach ML. Single or multiple frequency generators in on-going brain activity: A mechanistic whole-brain model of empirical MEG data. NeuroImage. 2017; 152:538–550. https://doi.org/10.1016/j.neuroimage.2017.03.023 PMID: 28315461 66. Doedel EJ, Champneys AR, Dercole F, Fairgrieve TF, Kuznetsov YA, Oldeman B, et al. AUTO-07P: Continuation and bifurcation software for ordinary differential equations; 2007. 67. Schecter S. The Saddle-Node Separatrix-Loop Bifurcation. SIAM Journal on Mathematical Analysis. 1987; 18(4):1142–1156. https://doi.org/10.1137/0518083 68. Fernandes LD, de Aguiar MAM. Turing patterns and apparent competition in predator-prey food webs on networks. Physical Review E. 2012; 86(5):056203. https://doi.org/10.1103/PhysRevE.86.056203 PMID: 23214853 69. Asllani M, Challenger JD, Pavone FS, Sacconi L, Fanelli D. The theory of pattern formation on directed networks. Nature Communications. 2014; 5(1):4517. https://doi.org/10.1038/ncomms5517 PMID: 25077521 70. Ercsey-Ravasz M, Markov N, Lamy C, Van Essen D, Knoblauch K, Toroczkai Z, et al. A Predictive Network Model of Cerebral Cortical Connectivity Based on a Distance Rule. Neuron. 2013; 80(1):184– 197. https://doi.org/10.1016/j.neuron.2013.07.036 PMID: 24094111 71. Deco G, Sanz Perl Y, Vuust P, Tagliazucchi E, Kennedy H, Kringelbach ML. Rare long-range cortical connections enhance human information processing. Current Biology. 2021; 31(20):4436–4448.e5. https://doi.org/10.1016/j.cub.2021.07.064 PMID: 34437842 72. Grimshaw R. Nonlinear Ordinary Differential Equations. Applied mathematics and engineering science texts. Taylor & Francis; 1991. Available from: https://books.google.es/books?id=yEWlegOzWxMC. 73. Illoul L, Lorong P. On some aspects of the CNEM implementation in 3D in order to simulate high speed machining or shearing. Computers and Structures. 2011; 89(11):940–958. https://doi.org/10.1016/j. compstruc.2011.01.018 74. Hindriks R, van Putten MJAM, Deco G. Intra-cortical propagation of EEG alpha oscillations. NeuroImage. 2014; 103:444–453. https://doi.org/10.1016/j.neuroimage.2014.08.027 PMID: 25168275 75. Zhang H, Watrous AJ, Patel A, Jacobs J. Theta and Alpha Oscillations Are Traveling Waves in the Human Neocortex. Neuron. 2018; 98(6):1269–1281.e4. https://doi.org/10.1016/j.neuron.2018.05.019 PMID: 29887341 76. Halgren M, Ulbert I, Bastuji H, Fabo ´D, Erőss L, Rey M, et al. The generation and propagation of the human alpha rhythm. Proceedings of the National Academy of Sciences. 2019; 116(47):23772– 23782. https://doi.org/10.1073/pnas.1913092116 PMID: 31685634 77. Pikovsky A, Politi A. Lyapunov Exponents: A Tool to Explore Complex Dynamics. Cambridge: Cambridge University Press; 2016. Available from: https://doi.org/10.1017/CBO9781139343473. 78. Hilborn RC. Chaos and Nonlinear Dynamics. 2nd ed. Oxford: Oxford University Press; 2000. Available from: https://doi.org/10.1093/acprof:oso/9780198507239.003.0006. 79. Afraimovich VS. Torus breakdown. Scholarpedia. 2007; 2(10):1933. https://doi.org/10.4249/ scholarpedia.1933 80. Kaplan JL, Yorke JA. Chaotic behavior of multidimensional difference equations. In: Peitgen HO, Walther HO, editors. Functional Differential Equations and Approximation of Fixed Points. Berlin, Heidelberg: Springer Berlin Heidelberg; 1979. p. 204–227. 81. Wendling F, Bellanger JJ, Bartolomei F, Chauvel P. Relevance of nonlinear lumped-parameter models in the analysis of depth-EEG epileptic signals. Biological Cybernetics. 2000; 83(4):367–378. https:// doi.org/10.1007/s004220000160 PMID: 11039701 PLOS COMPUTATIONAL BIOLOGY Transverse instabilities in large-scale brain networks PLOS Computational Biology | https://doi.org/10.1371/journal.pcbi.1010781 April 12, 2023 32 / 34 82. Jedynak M, Pons AJ, Garcia-Ojalvo J, Goodfellow M. Temporally correlated fluctuations drive epileptiform dynamics. NeuroImage. 2017; 146:188–196. https://doi.org/10.1016/j.neuroimage.2016.11.034 PMID: 27865920 83. Dumont G, Gutkin B. Macroscopic phase resetting-curves determine oscillatory coherence and signal transfer in inter-coupled neural circuits. PLOS Computational Biology. 2019; 15(5):1–34. https://doi. org/10.1371/journal.pcbi.1007019 PMID: 31071085 84. Sakaguchi H, Shinomoto S, Kuramoto Y. Phase Transitions and Their Bifurcation Analysis in a Large Population of Active Rotators with Mean-Field Coupling. Progress of Theoretical Physics. 1988; 79 (3):600–607. https://doi.org/10.1143/PTP.79.600 85. Zaks MA, Neiman AB, Feistel S, Schimansky-Geier L. Noise-Controlled Oscillations and Their Bifurcations in Coupled Phase Oscillators. Physical Review E. 2003; 68(6):066206. https://doi.org/10.1103/ PhysRevE.68.066206 PMID: 14754296 86. Childs LM, Strogatz SH. Stability Diagram for the Forced Kuramoto Model. Chaos: An Interdisciplinary Journal of Nonlinear Science. 2008; 18(4):043128. https://doi.org/10.1063/1.3049136 PMID: 19123638 87. Lafuerza LF, Colet P, Toral R. Nonuniversal Results Induced by Diversity Distribution in Coupled Excitable Systems. Physical Review Letters. 2010; 105(8):084101. https://doi.org/10.1103/PhysRevLett. 105.084101 PMID: 20868099 88. Pietras B, Devalle F, Roxin A, Daffertshofer A, Montbrio ´E. Exact firing rate model reveals the differential effects of chemical versus electrical synapses in spiking networks. Phys Rev E. 2019; 100:042412. https://doi.org/10.1103/PhysRevE.100.042412 PMID: 31771022 89. Clusella P, Ko ¨ksal-Erso ¨z E, Garcia-Ojalvo J, Ruffini G. Comparison between an Exact and a Heuristic Neural Mass Model with Second-Order Synapses. Biological Cybernetics. 2022. https://doi.org/10. 1007/s00422-022-00952-7 PMID: 36454267 90. Challenger JD, Burioni R, Fanelli D. Turing-like instabilities from a limit cycle. Phys Rev E. 2015; 92:022818. https://doi.org/10.1103/PhysRevE.92.022818 PMID: 26382465 91. Breakspear M, Roberts JA, Terry JR, Rodrigues S, Mahant N, Robinson PA. A Unifying Explanation of Primary Generalized Seizures Through Nonlinear Brain Modeling and Bifurcation Analysis. Cerebral Cortex. 2005; 16(9):1296–1313. https://doi.org/10.1093/cercor/bhj072 PMID: 16280462 92. Bressloff PC. Spatiotemporal dynamics of continuum neural fields. Journal of Physics A: Mathematical and Theoretical. 2011; 45(3):033001. https://doi.org/10.1088/1751-8113/45/3/033001 93. Deco G, Jirsa VK, Robinson PA, Breakspear M, Friston K. The Dynamic Brain: From Spiking Neurons to Neural Masses and Cortical Fields. PLOS Computational Biology. 2008; 4(8):1–35. https://doi.org/ 10.1371/journal.pcbi.1000092 PMID: 18769680 94. Esnaola-Acebes JM, Roxin A, Avitabile D, Montbrio ´E. Synchrony-induced modes of oscillation of a neural field model. Phys Rev E. 2017; 96:052407. https://doi.org/10.1103/PhysRevE.96.052407 PMID: 29347806 95. Coombes S, beim Graben P, Potthast R, Wright J, editors. Neural Fields. Springer Berlin Heidelberg; 2014. Available from: https://doi.org/10.1007/978-3-642-54593-1. 96. Spiegler A, Jirsa V. Systematic approximations of neural fields through networks of neural masses in the virtual brain. NeuroImage. 2013; 83:704–725. https://doi.org/10.1016/j.neuroimage.2013.06.018 PMID: 23774395 97. Jirsa V. Large Scale Brain Networks of Neural Fields. In: Coombes S, beim Graben P, Potthast R, Wright J, editors. Neural Fields. Springer Berlin Heidelberg; 2014. p. 417–432. 98. Atasoy S, Donnelly I, Pearson J. Human brain networks function in connectome-specific harmonic waves. Nature Communications. 2016; 7(1):10340. https://doi.org/10.1038/ncomms10340 PMID: 26792267 99. Abdelnour F, Dayan M, Devinsky O, Thesen T, Raj A. Functional brain connectivity is predictable from anatomic network’s Laplacian eigen-structure. NeuroImage. 2018; 172:728–739. https://doi.org/10. 1016/j.neuroimage.2018.02.016 PMID: 29454104 100. Glomb K, Kringelbach ML, Deco G, Hagmann P, Pearson J, Atasoy S. Functional harmonics reveal multi-dimensional basis functions underlying cortical organization. Cell Reports. 2021; 36(8):109554. https://doi.org/10.1016/j.celrep.2021.109554 PMID: 34433059 101. Escrichs A, Perl YS, Uribe C, Camara E, Tu¨rker B, Pyatigorskaya N, et al. Unifying Turbulent Dynamics Framework Distinguishes Different Brain States. Communications Biology. 2022; 5(1):638. https:// doi.org/10.1038/s42003-022-03576-6 PMID: 35768641 102. De Filippi E, Uribe C, Avila-Varela DS, Martı ´nez-Molina N, Gashaj V, Pritschet L, et al. The Menstrual Cycle Modulates Whole-Brain Turbulent Dynamics. Frontiers in Neuroscience. 2021; 15:753820. https://doi.org/10.3389/fnins.2021.753820 PMID: 34955718 PLOS COMPUTATIONAL BIOLOGY Transverse instabilities in large-scale brain networks PLOS Computational Biology | https://doi.org/10.1371/journal.pcbi.1010781 April 12, 2023 33 / 34 103. Cruzat J, Perl YS, Escrichs A, Vohryzek J, Timmermann C, Roseman L, et al. Effects of Classic Psychedelic Drugs on Turbulent Signatures in Brain Dynamics. Network Neuroscience. 2022; p. 1–21. 104. Goldman RI, Stern JM, Engel J, Cohen MS. Simultaneous EEG and fMRI of the Alpha Rhythm:. NeuroReport. 2002; 13(18):2487–2492. https://doi.org/10.1097/01.wnr.0000047685.08940.d0 PMID: 12499854 105. Moosmann M, Ritter P, Krastel I, Brink A, Thees S, Blankenburg F, et al. Correlates of Alpha Rhythm in Functional Magnetic Resonance Imaging and near Infrared Spectroscopy. NeuroImage. 2003; 20 (1):145–158. https://doi.org/10.1016/S1053-8119(03)00344-6 PMID: 14527577 106. Laufs H, Kleinschmidt A, Beyerle A, Eger E, Salek-Haddadi A, Preibisch C, et al. EEG-correlated fMRI of Human Alpha Activity. NeuroImage. 2003; 19(4):1463–1476. https://doi.org/10.1016/S1053-8119 (03)00286-6 PMID: 12948703 107. Feige B, Scheffler K, Esposito F, Di Salle F, Hennig J, Seifritz E. Cortical and Subcortical Correlates of Electroencephalographic Alpha Rhythm Modulation. Journal of Neurophysiology. 2005; 93(5):2864– 2872. https://doi.org/10.1152/jn.00721.2004 PMID: 15601739 108. Wang XJ. Macroscopic Gradients of Synaptic Excitation and Inhibition in the Neocortex. Nature Reviews Neuroscience. 2020; 21(3):169–178. https://doi.org/10.1038/s41583-020-0262-x PMID: 32029928 109. Wang XJ. Theory of the Multiregional Neocortex: Large-Scale Neural Dynamics and Distributed Cognition. Annual Review of Neuroscience. 2022; 45(1):533–560. https://doi.org/10.1146/annurev-neuro110920-035434 PMID: 35803587 110. van Rotterdam A, Lopes da Silva FH, van den Ende J, Viergever MA, Hermans AJ. A model of the spatial-temporal characteristics of the alpha rhythm. Bulletin of Mathematical Biology. 1982; 44(2):283– 305. https://doi.org/10.1007/BF02463252 PMID: 7074252 111. Seneta E. Non-negative Matrices and Markov Chains. Springer Series in Statistics. Springer New York; 2006. Available from: https://books.google.es/books?id=J3bsjqQBCZUC. 112. Wilkinson JH. The Algebraic Eigenvalue Problem. Monographs on numerical analysis. Clarendon Press; 1967. 113. Politi A. Lyapunov exponent. Scholarpedia. 2013; 8(3):2722. https://doi.org/10.4249/scholarpedia. 2722 114. Pikovsky AS, Rosenblum MG, Kurths J. Synchronization, a Universal Concept in Nonlinear Sciences. Cambridge: Cambridge University Press; 2001. 115. Clusella P, Montbrio ´E. Regular and sparse neuronal synchronization are described by identical mean field dynamics; 2022. Available from: https://arxiv.org/abs/2208.05515. 116. Galassi Mea. GNU Scientific Library Reference Manual; 2018. Available from: https://www.gnu.org/ software/gsl/. PLOS COMPUTATIONAL BIOLOGY Transverse instabilities in large-scale brain networks PLOS Computational Biology | https://doi.org/10.1371/journal.pcbi.1010781 April 12, 2023 34 / 34