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

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 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][59][60][61][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][48][49][50] to a comprehensive characterization of the emerging spatiotemporal behavior.

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 matrixW ¼ ðw ij Þ where i, j = 1, . . ., N. The connection weightsw ij are non-negative quantities that indicate the synaptic strength (average number of connections) from region j to 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 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 matrixW establishes that all network nodes receive the same amount of input, i.e., P N j¼1wij ¼ 1 for 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 p in 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 p small (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 p only 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 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 þ 2 supercritical 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) 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  (18) obtained by varying p and 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 p and with fixed � = 50. (e,f) Two-parameter bifurcation diagram of Eq (18) depending on the external input p and 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 (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: • Oscillatory activity at the theta or alpha ranges (or bistability between both), mostly delimited by the HB þ 2 and 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 matrixW . 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 ofW in 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 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 Λ α , 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).  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. In order to validate the emergence of heterogeneous dynamics in the system we perform numerical simulations of system (10) for different values of p and �, starting from initial conditions close to a homogeneous state. At each time step we compute the standard deviation across network nodes, 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 p and �. 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.

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 propagation, whereas it has been argued that travelling waves in the scalp are mediated by an intra-cortical mechanism [74]. In fact, travelling waves in the cortex are usually observed not at the whole-brain level but localized in smaller regions [5,75,76]. A reconciling approach could consist on applying the same analysis in a model closer to experimental conditions, based on short-range coupled NMMs in a small section of the cortex.
Vector fields in Fig 3(c) and 3(d) show coherent spatial propagation patterns visible to the naked eye. Nonetheless, in the model, the spatial support is discrete and the coupling among nodes is set by a complex network, thus local irregularities or outliers should be expected. We quantify the local directional agreement among neighbouring propagation vectors by means of a measure of local polarization a j (see Eq (34) in Methods). Fig 3(f) shows the average (black squares) and individual (small circles) polarizations for simulations A and B. In both cases the different velocity vectors display a good agreement with the propagation directions of its neighbouring regions, as indicated by a j being close to 1. Nonetheless, simulation B presents a few outliers, which we have checked to correspond to parietal brain regions close to the source of the wave, where the pattern spreads in many different directions.
The well-structured wave patterns in Fig 3(a) resemble some of the eigenmodes given by the diagonalization ofW depicted in Fig 2(a). Indeed, we expect the eigenmodes associated to positive Floquet exponents in the MSF to have a larger role on the final state of the system. Nonetheless, the random perturbations used to initialize the system might be more localized along a specific direction, making that direction more prominent among other unstable modes. Fig 3(g) shows this specific contribution for simulations A and B as approximated by Eq (29) (Methods). The perturbation used in A is mostly being expanded through the second eigenvector (see Fig 2(a)), which shows a left-to-right division in close agreement to the observed pattern. Instead, perturbation B is mostly being expanded according to the third structural eigenvector, leading to the back-to-front wave propagation.

Chaotic wave dynamics
Numerical analysis displayed in Fig 2(f) and 2(g) shows a large region of chaotic dynamics. Here we analyze an instance of such states. Lyapunov exponents are a common tool to characterize the complexity of a chaotic regime (see Methods) [77]. We set parameter values to p = 170 and � = 20, for which at least 5 Lyapunov exponents are positive, an indication of highdimensional chaos. Single nodes still oscillate at the alpha range (�9Hz), however, the phase differences between brain areas change over time. Fig 4(a) shows consecutive time snapshots of the difference between the phase of each region ϕ i and the collective phase C, with the corresponding propagation direction vector field depicted in Fig 4(b). Such phase gradient preserves certain levels of locally coherent structure, although the scenario looks much more complex than the periodic travelling waves presented in the previous section (see also S3 Movie).
In order to understand the main features of these chaotic patterns, we first analyze the level of phase synchrony in the network. By observing the time evolution of the phase differences (see Fig 4(c)), we can appreciate that the dynamics alternates between periods in which all nodes evolve with similar phases (all phase differences are distributed around zero), with periods in which a group of nodes (but not all) lose synchrony. Such type of breathing dynamics can be monitored by means of the Kuramoto order parameter R (see Methods for a definition). This can be seen in Fig 4(d), which shows the irregular oscillatory activity of R. These results indicate that the level of collective synchronization in the network itself displays chaotic behavior. The regularity of both the phase differences and R can be quantitatively analyzed by means of their corresponding autocorrelation functions (Fig 4(e)). The quick decay of both autocorrelations illustrates the short-time scale of the dynamics, which agrees with that indicated by the inverse of the largest Lyapunov exponent (1/λ 1 ' 0.93s). Hence, overall, the synchronization patterns in the network, either pairwise or collective, are mediated by the spontaneous appearance of a slow modulation of the alpha activity.
We next characterize this chaotic wave propagation using different measures. Fig 4(f) shows the average and standard deviation of the propagation speed for each node. The velocity range is similar to that of the periodic waves (1-8m/s), with a few outliers with larger velocities. These outliers are not present on the time-averaged distribution (right panel in Fig 4(f)), thus indicating that they correspond to sporadic fluctuations. The local polarization (Fig 4(g)) exhibits a similar behavior: Instantaneous snapshots of a j display a large variability on the phase coherence across regions, but with a time average concentrated around ha j i � 0.5 for all nodes. Finally, we quantify the change of the propagation direction over time using the mean resultant length ρ j = khζ j /z j ik (see Methods), which measures the directional coherence of single nodes, being 0 for completely random directions, and 1 for quenched directionality (as in periodic travelling waves). The low values of ρ j in Fig 4(h) indicate that most nodes propagate their dynamics to statistically almost any direction over the entire simulation. Altogether, this is a manifestation of short-lived and locally coherent wave dynamics, which travel through the network without a well-identified pattern.

Onset and multistability of irregular dynamics
In this section we characterize the emergence the spatiotemporally chaotic dynamics displayed by our large-scale brain model, as well the effect of transverse instabilities on the frequency of the oscillations. We perform different sets of numerical simulations for fixed � = 50 and varying p, starting from random initial conditions. Fig 5(a) shows the 5 largest Lyapunov exponents of the system as p is varied (colored thin lines, left axis). The transverse instability takes place at p ' 265.5 (right-most vertical dashed line in the figure), and the first instances of chaos emerge for values of p smaller than p ' 184.5 (first vertical dashed line in Fig 5(a)). Between these two values, the behavior of the system is either periodic (single zero LE) or quasiperiodic (two LEs equal to zero). Chaotic dynamics appear rather abruptly from the quasiperiodic states as p is decreased.
We therefore conclude that the route to chaos is that of a torus breakdown [78,79]. Examples of periodic, quasiperiodic, and chaotic dynamics are displayed in Fig 5(b). In this panel, each plot shows the time series of the sample mean voltage,  then breaks into faster irregular modulations of the alpha rhythm at the chaotic state (bottom row of Fig 5(b)).
After the onset of chaotic dynamics at p ' 265.5, further decrease of p leads to consecutive Lyapunov exponents becoming positive, thus increasing the dimensionality of the chaotic attractor. In fact, although we only display 5 Lyapunov exponents, many more directions become unstable. We estimate the (fractal) dimension of the chaotic attractor in terms of the Kaplan-Yorke dimension, D KY [77,80] (see Methods), which quantifies the effective number of degrees of freedom that characterizes the irregular dynamics (black thick line in Fig 5(a), right axis). This analysis reveals the high-dimensionality of the chaotic regime, which can reach up to 100 for low p values. For even lower p the complex oscillatory state suddenly vanishes, and the dynamics converge to the stable homogeneous fixed point.
Does the onset of heterogeneous states modify the oscillatory frequency of the system? To answer this question we analyze the leading frequency of each network node upon varying p. The result is shown in Fig 5(c), where dots indicate the peak frequency of each brain area, as determined by the maximum of the power spectral density in each case. For each value of p, results from 5 different simulations are shown in gray, with results from a single randomly chosen simulation depicted in red. Before the synchronous alpha activity loses stability (right-most vertical dashed line in Fig 5(c), p > 265.5), all nodes in all simulations evolve according to the frequency given by the analysis of the self-coupled system (18), obtained in Fig 1(d) and indicated by a black curve. Between the transverse instability of the homogeneous state and the emergence of chaos (184.5 � p � 265.5, area between the two vertical dashed lines in the figure), all brain areas in a given simulation oscillate with the same leading frequency, even though different sets of simulations fall in different attractors with slightly different oscillation frequencies (due to multistability between different types of periodic states, as at the case analyzed in Fig 3). In the chaotic regime (p � 184.5) this situation changes, and within each simulation different nodes might oscillate at different peak frequencies, which cover a broader range as the chaotic state increases dimensionality. In any case, it is worth highlighting that oscillatory activity remains at the alpha range, around 9Hz, for most values of p. This frequency range differs from that of the underlying homogeneous state, which transitions to the theta regime and finally vanishes (see black curve).
As shown in Fig 5(c), vestiges of theta activity exist only for values of p close to the steadystate regime (e.g., p ' 70). In those cases, brain regions oscillate intermittently between the alpha and theta frequency bands. Fig 5(d) shows an example of such multifrequency dynamics for p = 70 and � = 50, with the corresponding power spectra displayed in Fig 5(e). Similar multifrequency dynamics are usually associated to the role of stochastic input inducing jumps between the two oscillatory states of the uncoupled Jansen model [81,82]. Here instead, it is the collective deterministic chaos what causes a two-peak spectrum for single nodes.

Non-normalized connectivity
All results shown so far stem from the row-normalized connectivity matrixW given by Eq (14), which allows for a semi-analytical treatment of the model. Next, we turn to a numerical analysis of the model using the non-normalized network connectivity W directly as obtained from the data. In this case one cannot ensure the existence of a homogeneous manifold, and thus the fixed points and their stability should be assessed by solving a system of 540 nonlinear equations and analyzing the corresponding Jacobian. Here we avoid such a full stability analysis of the system, since direct simulations are enough to compare with our previous results.
We analyze numerically Eq (10) using the SC matrix W. Therefore, we replace the input Eq (11) in Eq (10) by where hsi is the average node strength, i.e., The scaling factor hsi −1 ensures that the distribution of incoming connections is centered at 1, so that the effect of the coupling strength � in the non-normalized model is comparable to that the row-normalized case studied so far. Fig 6(a) and 6(b) show the classification of dynamical states as obtained from the two largest Lyapunov exponents in numerical simulations. The first figure shows results obtained initializing the system close to a homogeneous state, whereas the second corresponds to random initial conditions. Continuous lines show the bifurcations of Eq (18) as in Fig 2. The diagrams reveal a scenario very similar to that of the normalized topology, with a central large region displaying different types of oscillatory states bounded by two white regions corresponding to a lowactivity (left) and high-activity (right) fixed points. Also, the presence of hyperchaos (red region in the plots) is ubiquitous, and shares a small island of bistability with periodic oscillatory states which contrasts with the large bistable areas present in the normalized system ( Fig  2). Additionally, the non-normalized system displays a large region of quasiperiodicity and a small island of low-dimensional chaos not present in  Despite the quantitative dissimilarities between the dynamical landscapes of the row-normalized and non-normalized systems described above, there is a fundamental difference between using W orW with regard the intrinsic dynamics of the system: In the non-normalized system all states are always heterogeneous. Fig 7(c), for instance, shows a simulation for parameter values close to the onset of oscillations (p = 320 and � = 50).
This oscillatory state arises from an underlying unstable heterogeneous fixed point that we obtain by solving the nonlinear fixed point equations of the system (see Methods). Indeed, Fig  6(d) shows that the time averaged mean membrane potential hv i (t)i, is practically identical to the underlying steady state v ð0Þ i for large enough p (blue circles). In a chaotic state (p = 170, � = 50, red squares in Fig 6(d)) the center of the oscillations shifts from the fixed point, but the two quantities still show a high degree of correlation. Remarkably, the heterogeneous unstable fixed points remain distributed around the homogeneous steady state of the model with normalized connectivity for same parameter values (vertical dashed lines in Fig 6(d)). Hence, the use of the matrix W does not seem to fundamentally alter the nature of the different steady states.
Next we study to what extent these oscillatory states depend on the distribution of inputs s i . Fig 6(e), 6(f) and 6(g) show the focus, amplitude, and phase difference of each node plotted against the node strength for the periodic (blue circles) and chaotic (red squares) cases. First, the focus of the oscillations follows a direct relation with the input of each node (Fig 6(e)), for 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)

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][85][86][87][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 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][59][60][61][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][92][93][94][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 functional and resting brain activity [98][99][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][102][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][105][106][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.
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]).

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 W was 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 h e ðtÞ ¼ Aate À at HðtÞ where H is the Heaviside step function, A and B are the PSP amplitudes, and a and b quantify 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 differential equation Analogously, a inhibitory firing rate generates a inhibitory post-synaptic potential (iPSP) according to On the other hand, a population with mean membrane potential v(t) generates a firing rate according to the sigmoid function 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: 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 N nodes, 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 i read [35] _ y 0;i ðtÞ ¼ y 3;i ðtÞ _ y 1;i ðtÞ ¼ y 4;i ðtÞ _ y 2;i ðtÞ ¼ y 5;i ðtÞ _ y 3;i ðtÞ ¼ Aa Sigm½y 1;i ðtÞ À y 2;i ðtÞ� À 2ay 3;i ðtÞ À a 2 y 0;i ðtÞ _ y 4;i ðtÞ ¼ AafI i ðtÞ þ C 2 Sigm½C 1 y 0;i ðtÞ�g À 2ay 4;i ðtÞ À a 2 y 1;i ðtÞ _ y 5;i ðtÞ ¼ BbC 4 Sigm½C 3 y 0;i ðtÞ� À 2by 5;i ðtÞ À b 2 y 2;i ðtÞ : 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 common constant firing rate p. Altogether, I i (t) takes the form of the sum of two independent terms: where � is the coupling strength, andW ¼ ðw ij Þ is the row-normalized structural connectivity matrix (see next section). We study the system under the influence of the external driving firing rate p and 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 × N structural 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 × N matrix whose non-zero entries are the sum of incoming connections to each node: is the in-strength of node i. Then we consider the row-normalized SC matrix whose elements arew ij ≔ w ij =s i , i.e.,W is obtained dividing each row of W by its sum. Therefore, the sum of the elements on each row equals unity i.e., 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 2 WD À 1 2 . Hence Z is diagonalizable and has real eigenvalues. Using that Therefore, Z andW are similar (in the mathematical sense), i.e., they share the same eigenvalues. The fact that the eigenvalues oW are 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 dependent-in which all nodes of the network behave identically.

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 N and 6N dimensional vectors and matrices, and regular symbols for 6-dimensional vectors, 6-dimensional matrices, and scalar quantities. All scalar quantities except the time t and the coupling strength � have a subscript.
Let y ð0Þ ≔ ðy ð0Þ 0 ; . . . ; y ð0Þ 5 Þ T be a solution of (18). Then, the 6N-dimensional vector where each component δy k,j is the perturbation acting on the variable k of node j in 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 where J is the full 6N × 6N Jacobian of system (10), evaluated at y (0) . Notice that J can be written as where � denotes the Kronecker product, J is the 6 × 6 Jacobian of the uncoupled Jansen Eq (9), I N is the N × N identity matrix, and K = (k ij ) is a 6 × 6 matrix defined by 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 δy in an adequate coordinate system.
Let Φ ðaÞ ¼ ð� ðaÞ 1 ; . . . ; � ðaÞ N Þ T be a normalized eigenvector ofW associated with the eigenvalue Λ α for α = 1, . . ., N, so thatW The set of eigenvectors {Φ (1) , . . ., Φ (N) } constitute a basis of the vector space R N , thus we can express the perturbation δy of the homogeneous solution as a linear combination of such 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 where we have used the diagonalization of the connectivity matrix Eq (24). Now, making use of the linear independence of the eigenvectors fΦ ðaÞ g N a¼1 one obtains that the evolution of u (α) (t) becomes independent for each α = 1, . . ., N through the relation where J ðy ð0Þ ; LÞ ≔ J þ �L a K is 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, 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 J is 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 J as 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 δy we compute the quantity c a ¼ m a ku ð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).
• 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).
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 where j is the LE for which

Synchronization and wave-propagation analysis
We quantify the synchronization of the network by extracting the phase variables ϕ j (t) from each node's dynamics v j (t) using the Hilbert transform [114], i.e., where v ðHÞ j ðtÞ is the Hilbert transform of v j (t), v ðHÞ j ðtÞ ¼ 1 p P:V: with P.V. indicating that the integral corresponds to the Cauchy principal value. Then we can compute the complex quantity in which R 2 [0, 1] is the usual Kuramoto order parameter, and C is the collective phase variable [53,114]. From a given phase gradient in space, we extract the propagation vector of each node ζ j by following the same procedure as in [30] (see also [5]). Briefly, the method consists on numerically finding the temporal and spatial derivatives of ϕ and then solving for ζ. The numerical differentiation of the phase is obtained through a constrained natural element method [73]. The instantaneous speed of each node is then given by the modulus of the velocity vector, z ≔ kζk.
In order to measure the degree of alignment of a certain direction vector with neighbouring nodes, we use a measure of local polarization. Let N j be the set of nodes that are located at less than 40mm than node j. Then the local polarization of j is defined as A value close to 1 indicates a perfect alignment of node j with its neighbours, whereas a value close to zero indicates a nearly isotropic propagation at j. Finally, in order to assess how much the direction of propagation changes over time we compute the directional coherence as the modulus of the time-average normalized vectors Notice that ρ j corresponds to the mean-resultant length of the temporal distribution of directions ζ j /z j .

Next generation neural mass model
The firing rate equations for populations of quadratic integrate-and-fire neurons are based on an exact mean-field theory derived by Montbrió et. al. [17]. This theory has been recently shown to apply also for neurons subject to Cauchy white noise [115]. Here we consider a PING mechanism in which the oscillatory activity of each brain region is in the gamma range, similarly to the scenario studied in [83]. Hence, at each network node j we consider an excitatory and a inhibitory populations of QIF neurons, each characterized by its firing rate, r, and mean membrane potential, v. We also consider delta pulses for the PSP, as a means to illustrate that the shape of the PSP is not crucial for the transverse instability of the synchronized state. Hence, the neural activity of brain region j is given by t_ r e;j ¼ D tp þ 2r e;j v e;j t _ v e;j ¼ Z e À ðpr e;j tÞ 2 þ v 2 e;j þ tðJ ee r e;j À J ie r i;j þ I j Þ where r e,j and v e,j (resp. r e,j , v e,j ) are the firing rate and mean-membrane potential of the excitatory (resp. inhibitory) population of region j. The explanation and value of the different system parameters are given in Table 2. This setup corresponds to the PING mechanism for gamma activity, and shows that a single uncoupled brain area displays activity at a frequency around 40Hz [83]. We have verified that the results are robust upon considering other parameter choices and, in particular, when including recurrent inhibition (J ii > 0). As in the network of Jansen NMMs, we consider that different regions interact only through excitation, i.e., Following the same argument as for the Jansen system, if all nodes of the network evolve identically, then the dynamics of each network node follows t_ r e ¼ D tp þ 2r e v e t _ v e ¼ Z e À ðpr e tÞ 2 þ v 2 e þ tððJ ee þ �Þr e À J ie r i Þ

Numerical simulations and root-finding algorithm
Numerical simulations of the systems have been performed with a fourth-order Runge-Kutta algorithm with time step Δt = 10 −3 . The C code to simulate the model composed of Jansen NMMs is available at www.github.com/pclus/transverse-instabilities. Most of the simulations have been run for 1000s after discarding an extra 1000s of relaxation time. The only exception is Fig 5(a), for which we used a relaxation time of 10 4 s in order to achieve a good convergence of the zero Lyapunov-exponents needed to compute the Kaplan-Yorke dimension. Two different types of initial conditions have been considered: • Close to homogeneous: we select a point in the trajectory of the homogeneous dynamics of the system using the equations of the homogeneous manifold, Eqs (18) and (38), and adding an independent random perturbation to each variable, uniformly distributed between −10 −3 to 10 −3 .
• Random: all variables of the system randomly distributed between −1 and 1.
In order to obtain the stationary states of the unnormalized network topology W in Section entitled "Normalized connectivity" above, we use the adaptive Newton-Raphson multidimensional root-finding algorithm provided by the GNU Scientific Library [116].