Propagation of transient explosive synchronization in a mesoscale mouse brain network model of epilepsy

Epileptic seizures, characterized by bursts of excessive neuronal synchronization, which usually self-initiate and self-terminate, are considered as a dynamical disease of brain networks (Da Silva et al., 2003). Seizures can be classified into distinct subtypes, broadly including those that are confined to a circumscribed area (focal) and those that involve larger sections of the brain (generalized) (Sarmast et al., 2020). A wide range of microscopic mechanisms contribute to this limited repertoire of seizure types (Richardson, 2012). Although seizures can originate from different brain regions in patients, they may still manifest similar macroscopic features, as seen in EEG recordings (Bromfield et al., 2006; Gleichgerrcht et al., 2015; Richardson, 2012). Furthermore, studies have suggested that seizures with similar microscopic mechanisms can present as either focal or generalized depending on the macroscopic network structure (Terry et al., 2012). This decoupling between microscopic and macroscopic dynamics underscores the importance of directly modeling emergent macroscopic properties of seizures and highlights the significance of adopting a network-level approach to studying epilepsy, which has also been recognized by the International League Against Epilepsy (Berg et al., 2010).

Experimental evidence shows that network structure alone is not sufficient, but the dynamics supported by it also play an important role in seizure generation and propagation in a brain network (Chowdhury et al., 2014; Petkov et al., 2014). Seizures have been hypothesized to exist in a bistable regime of dynamical networks that exhibit multiple stable states: normal (unsynchronized) and abnormal (hypersynchronized). In such a system, random fluctuations (noise) or resource availability can transition the network between the different states, giving rise to transient hypersynchronized activity seen during epileptic seizures (Da Silva et al., 2003; Jirsa et al., 2014).

The phenomenon of explosive synchronization (ES) that is widely studied in complex systems and network science can provide a general framework to understand the role of network structure in facilitating seizure dynamics (Gerster et al., 2020; Gomez-Gardenes et al., 2011; Leyva et al., 2013; Zhang et al., 2015). ES is characterized by first-order, discontinuous, and irreversible transitions between globally coherent and incoherent states. These features are highly relevant to seizure dynamics, which also show signatures of critical transitions at both onset and termination across multiple spatial scales (Kramer & Cash, 2012; McSharry et al., 2003). Consequently, ES models have been employed to study abrupt transitions in brain networks (Kim et al., 2022; Wang et al., 2017). Moreover, complete brain coverage recordings, which can elucidate the dynamics of generalized seizures, are lacking in both patient and animal models. Thus, integrating ES with biological networks enables the study of seizure-like synchronization dynamics at the whole-brain level.

Although ES in complex networks has been successfully modeled using two common microscopic mechanisms—the presence of microscopic correlation features, such as frequency-degree coupling (FDC) and adaptive coupling (Gomez-Gardenes et al., 2011; Zhang et al., 2015)—it does not explain the transient and recurrent nature of seizures. Experimental and computational studies have linked this transient nature of seizures with the dynamics of energy metabolism (Bromfield et al., 2006; Jirsa et al., 2014). Several indirect measures of energy metabolism (ATP use), such as oxygen, intracellular NADH levels, and extracellular ions, have been known to change slowly across the timescale of seizure, and they have consequently been used to model onset and termination of epileptic seizures (Jirsa et al., 2014). Consistent with this observation, a recent study has shown the occurrence of transient ES (tES) in a resource-constrained scale-free network with FDC (Frolov & Hramov, 2021), where the time-varying nature of resource consumption is shown to cause the transient behavior.

While resource-constrained networks with FDC exhibit tES for certain scale-free networks, there are several other frequently occurring families of network structures (Figure 1, panels A–C) across which this mechanism does not appear to generalize. Especially in the context of neural dynamics, the network structure of brains (Figure 1D) is often found to show characteristics of small-world networks (SWNs) as well as scale-free networks (SFNs) (Figure 1, panels E and F). Moreover, adaptive coupling schemes have been shown to exhibit ES more generally across several network topologies, and in fact ES has been suggested to be a generic property of networks with adaptive coupling (Zhang et al., 2015). Adaptive coupling is biologically plausible and has often been used to model interaction in biological systems (Lopes et al., 2023). Specifically, it models spike-timing-dependent plasticity (STDP) for populations of neurons, where connectivity changes are influenced by the level of synchrony between populations (Duchet et al., 2023). Therefore in this study, we combine resource constraint with adaptive coupling in a model (Figure 1G) that can manifest tES across several types of network structures, including classic small-world and scale-free networks as well as a real biological network (a mesoscale mouse brain network, MBN, obtained from the Allen Institute public dataset; Oh et al., 2014).

Depending on the resource availability, our model exhibits desynchronized activity, bistability of desynchronized and hypersynchronized activity (i.e., tES), as well as steady-state hypersynchronized activity in SFNs, SWNs, and the MBN. Furthermore, during the sudden transition to the synchronized state, we observe a wave-like propagation of synchronization across subnetworks within the MBN, beginning with cortico-thalamic subnetworks, followed by subcortical and deeper subnetworks. We also develop a novel algorithm to analyze how the synchronization propagates across individual nodes (brain areas) in the MBN and identify key brain areas that may be responsible for initiation, sustenance, and propagation of the hypersynchronized state. Our results agree with several observations from experimental studies, suggesting that a few key parameters can successfully capture the network-level phenomenology of seizure dynamics. Finally, the model allows us to study the relationship between the hypersynchronized state and the resource consumption to recovery rate ratio. Specifically, the model predicts an optimal intermediate ratio for which the likelihood of tES, that is, epileptic attacks, is minimal. This and related predictions of our model should be directly testable in experiments.

Our model consists of N sinusoidal oscillators that form the nodes of a network. Following the Kuramoto model, the interactions between connected oscillators depend on their phase difference. The interaction strength is determined by the structural weight of the connection and is further modulated by both the synchronization levels of neighboring nodes and the availability of excitability resources (Figure 1G).

We begin by investigating the properties of our model in an SWN. For this, we generate an SWN comprising 400 nodes while maintaining parameters such as average degree (〈k〉 = 40), average path length (L_p = 1.97), and clustering coefficient (C_p = 0.37) similar to the MBN for later comparison (Figure 1A).

We first characterize the resource dependence of the network dynamics in the absence of resource dynamics (Figure 2A). Thus, in Equation 1, λ_i is replaced by Λ, a fixed resource available to each node at all times. We simulate this model with varying Λ and observe the steady-state behavior. For this, we begin with Λ = 0, adiabatically increase (decrease) Λ with increment (decrement) of ΔΛ = 0.003, simulate the model for 1,000 time steps, and compute the stationary value of global synchrony (R) (see Methods section), going up to Λ = 0.12. For very small Λ the network exhibits normal activity, and for very high values of Λ, it goes into the hypersynchronized state.

Interestingly for intermediate values, as Λ is slowly varied, we observe an abrupt first-order irreversible transition, with the presence of a hysteresis region. Depending on the direction of change of Λ, or in other words, depending on the current state of the dynamics, the network goes into either the normal or the hypersynchronized state (Figure 2A).

Such existence of hysteresis has been shown to give rise to tES when resource constraint is imposed (Equation 3) (Frolov & Hramov, 2021). Mean resource level, 〈λ〉, in resource-constrained models, functions similarly to fixed resources (Λ) in an unconstrained model, except that it is a dynamic quantity. Its value can oscillate around the hysteresis region, causing the network to transition between incoherent and coherent states (Figure 2C). When sufficient resources become available, that is, when the average level of resources in the network, denoted as 〈λ〉, crosses a certain threshold, called the forward tipping point (Λ_c(fw); Figure 2A, red), the network shifts into a state of very high synchronization. This change is marked by a sudden drop in available resources because of increased consumption (Figure 2, panels E and F). Depending on the resource availability, the network can stay in this highly synchronized state for only a limited time before it runs out of resources and falls back to a less synchronized state. This backward shift happens when the average resource level 〈λ〉 drops below another threshold, known as the backward tipping point (Λ_c(bw); refer to Figure 2A, blue). While in this less synchronized (or incoherent) state, the system replenishes its resources. Once there are enough resources again to reach the forward tipping threshold, the cycle can repeat, with the network moving back into the highly synchronized state.

Towards testing this hypothesized mechanism for tES, we first identify the parameter range for which the full model shows bistability; that is, the network spends time in two (meta)stable states (Figure 2B). We impose resource constraint at a consumption rate β ( = 0.002) and again simulate the full model for λ_o, the size of the resource bath, varying between 0.01 and 0.3 with an increment of 0.01. For each value of λ_o, we simulate the model for 1,000 time steps, and observe the range of values taken by the global synchrony parameter over the simulated period. The resulting bifurcation diagram reveals that the SWN exhibits a globally incoherent state for λ_o < 0.095, where resources are too limited to allow hypersynchronization, and a hypersynchronized state for λ_o > 0.21. For intermediate values of λ_o the network shows a coexistence of both states, reflecting the presence of tES, as seen in the time series of the global synchrony parameter (Figure 2, panels E and F). Consistent with the proposed mechanism, the transitions to and from the hypersynchronized state occur when the mean resource availability across all nodes is close to the corresponding critical values of Λ, as revealed in the state space trajectory of the system (Figure 2C).

Since the transition back to the incoherent state occurs because of resource depletion, we expect the network to spend a longer time in the hypersynchronized state as the size of the resource bath, λ_o, increases. This prediction is supported in our simulation results (Figure 2D) as well as with a simple analytical calculation (see the Supporting Information).

Finally, we note that the phenomenon of ES or tES is not observed with SWN when correlation feature-based connectivity, such as FDC, is used (Frolov & Hramov, 2021).

Next, we use our model to study tES in a real biological neural network, using the mouse brain mesoscale connectivity data from the Allen Brain Atlas. The mesoscale atlas is constructed by injecting viral vectors to trace axonal projections across pairs of brain regions in mice (Oh et al., 2014). The dataset consists of detailed and accurate connectivity information across 426 brain areas spanning both hemispheres (see Methods section) in healthy mice. The resulting network (Figure 1D) comprises 11,000 directed edges, with weights rescaled between 0 and 1.

We repeat the analyses from SWN on MBN with an added nuance: while we assumed the SWNs to be binary undirected networks, the dataset we use allows us to define the MBN as a weighted directed network. To account for this, we use a modified version of Equations 1 and 2 (see the Supporting Information). We analyze all three variants of the MBN: binary-undirected, binary-directed, and weighted-directed. The results presented below refer to the most complete weighted-directed variant unless mentioned otherwise, while the results for the other two variants are qualitatively similar (Supporting Information Figure S4).

We again begin by characterizing the system with fixed resource availability (see Methods section for details). Surprisingly, unlike the SWN which showed hysteresis, the MBN shows bistability for intermediate values of fixed resources (Λ ∈ [2.58, 2.65]) (Figure 3A). Even binary-undirected and binary-directed versions of MBN show hysteresis but not bistability for fixed resources (Supporting Information Figures S4A and S4B). It is important to note that when examining a network with fixed resources, the hysteresis region indicates the potential for tES, but does not guarantee it. On the other hand, the presence of the bistability region directly confirms the existence of tES for the weighted-directed MBN even without the need for resource constraint. We further confirm the existence of tES with simulations of the full model that includes resource constraint (Figure 3B–D, Supporting Information Figure S3). Note that with the FDC model with or without resource constraints, we are unable to induce tES in the MBN (Supporting Information Figures S4F, S4G).

Our analysis of the MBN recapitulates the observations from SWN, including the increasing time spent in the synchronized state with increasing size of resource bath (Figure 3B, inset). An interesting deviation is that the synchronized state is only partially synchronized, as reflected in the global synchrony parameter only reaching up to 0.6 instead of 1 (Figure 3, panels A–D). This happens because certain subnetworks within the MBN never participate in the synchronization. We thus define actively participating nodes as those with a local synchrony greater than the threshold value of 0.7 during the partially synchronized state. We then compare the state space trajectory obtained by averaging the available resources over all nodes versus averaging only over the actively participating nodes (Supporting Information Figure S2). The state space trajectory shifts to the right in the former case, suggesting that inactive nodes are pushing the average available resources, 〈λ〉, to a higher value. We thus conclude that only the nodes that actively participate in the synchronization cluster increase their energy consumption during tES.

Another deviation from SWNs is that the average resource level (even after restricting to participating nodes) at the time of transitions is observed to be higher than the corresponding resource level in the fixed resource model (Figure 3, panels A and C). We hypothesize that even within the synchronized cluster, there may be a core subcluster that drives tES, for which the average energies at transition may be lower, but the peripheral nodes in the synchronized cluster, with higher energy availability, push the apparent average transition energy levels higher.

This does not seem to be the case, however. We conclude that the higher complexity of biological networks, perhaps owing to their modular, hierarchical, and heterogeneous nature compared with the SWN, makes the relationship between the fixed resource and resource-constrained dynamics less predictable.

In the following sections we investigate the dynamics of propagation of synchronization in the MBN.

We group the 426 nodes in MBN into distinct communities purely based on the network structure, using the Louvain algorithm (see Methods section) (Blondel et al., 2008). We hypothesize that such structure-based partitioning of the nodes will group together nodes that are highly likely to form a synchronization cluster. The process naturally partitions the MBN into eight communities that can be clubbed into three broad classes: cortico-thalamic (four communities), subcortical (three communities), and hindbrain (Figure 4A).

We then generate a set of 52 transitions from the incoherent to the hypersynchronized state by running four long simulations, with fixed λ_o = 2.87, but with different initial conditions (see Methods section). Each transition is characterized by the presence of an “onset index,” a point where nearly complete desynchronization occurs just before the abrupt transition is about to begin (Figure 4B, inset), subsequent to which synchronization rapidly increases. To study the spatial propagation of synchronization within the short temporal window in which it occurs, we analyze the transient dynamics over a 50-time-unit window following the onset index.

We first study the dynamics of synchronization during the transient window within each community, by computing the intramodular synchrony (a measure of phase alignment across nodes within the community; see Methods section) averaged over all 52 transients. Based on the temporal evolution of the intramodular synchrony, we can group the eight communities into three distinct cohorts that perfectly overlap with the structural classes defined earlier: cortico-thalamic, subcortical, and hindbrain (Figure 4C). The steady-state intramodular synchrony is significantly higher for the cortico-thalamic communities, followed by subcortical, and finally hindbrain, which shows no internal synchronization (p < 10⁻⁴ for all the three pairs).

At the start of the transient process, the orbital, olfactory, and sensorimotor areas exhibit significantly higher intramodular synchrony than the rest of the communities (p < 0.05 for each pair; Supporting Information Figure S8), and also show a steeper rate of increase (Supporting Information Figure S8), indicating their potential role in initiating and driving the abrupt transition. As the transient progresses, orbital and olfactory areas consistently maintain a significantly higher level of synchronization compared with the rest (p < 0.05; Supporting Information Figure S8). These results suggest that, while the onset sites of tES can be the orbital, olfactory, or sensorimotor areas, it is likely that the orbital and olfactory areas play a dominant role in driving the transition throughout the entire duration of the transient process. We note that, while intramodular synchrony in the visual area does not seem to significantly differ from the orbital area (Supporting Information Figure S8), it shows much higher trial-to-trial variation in terms of its participation. This reinforces previous observations in the literature that the visual area is not critical for the propagation of synchronization (Nersesyan et al., 2004).

The dynamics of synchronization propagation can be understood by analyzing the intermodular synchrony between community pairs (see Methods section). Consistent with what we found earlier, ES initiates in the orbital and sensorimotor areas (Figure 4D, t = 2.2), and then spreads across all cortical regions (Figure 4D, t = 4.4). Notably, the hippocampal/midbrain areas synchronize with cortical areas before synchronizing among themselves (Figure 4D, t = 4.4, 6, 50), indicating that the cortical areas are driving their synchronization. At steady state (t = 50), all cortical, thalamic, and subcortical regions achieve synchrony, while the hindbrain exhibits minimal participation throughout. Overall, these findings suggest a hierarchical synchronization process, with cortical areas potentially driving synchronization across subcortical regions.

Going beyond the community level, below we assess synchronization propagation at a single node level. Assessing the synchronization between individual nodes typically involves computing the correlation of activity of the pair over short time windows (Schmidt et al., 2015; Zhang et al., 2014). However, given the extremely short-lived transition window, this technique cannot provide sufficient temporal resolution for studying abrupt transitions. We therefore came up with a novel algorithm that we call the synchronization cluster tracking algorithm (SCTA) (see Methods section) for this analysis. The SCTA employs a two-step procedure for quantifying the participation of individual nodes in the synchronization process: (1) generating synchronization clusters by starting from a seed node and expanding the cluster until the local synchrony drops below a fixed threshold; and (2) tracking the synchronization clusters progressively through time to get a “cluster lineage” for each node (see Methods section). We apply SCTA to all 52 transients with a local synchrony threshold of 0.5 for the subsequent analysis, although the choice of threshold does not qualitatively change the results (Supporting Information Figure S5).

The transition typically begins with multiple small synchronization clusters (cluster size < 10) during the early phase (Figure 5A, t ≈ 17). As the transition progresses, we observe the emergence of a “main synchronization cluster” that quickly comes to dominate in size (Figure 5A) throughout the rest of the transition. The size of the main cluster varies as the small peripheral clusters or individual nodes continue to join or leave it. This indicates that while the specific onset sites may vary across different trials initially, once the abrupt transition begins, it rapidly spreads from that initial onset site to encompass a broad set of core nodes that hold it together.

To test the “core” nodes hypothesis, we quantify the time spent by individual nodes in the main synchronization cluster across the 52 transients. We find that although a majority of the nodes exhibit participation in the main synchronization cluster, certain nodes consistently spend a significant amount of time in the main cluster (median time spent > 35, SD < 10, Figure 5B). This core spans across multiple brain areas including the cortex (19 nodes), hippocampus (6), striatum (9), thalamus (10), and midbrain (3). In contrast, nodes with high variability likely represent peripheral nodes that frequently attach and detach from the main cluster, contributing to the observed variation in cluster size.

This analysis allows us to hypothesize the existence of driver nodes for hypersynchronization, as nodes with consistent, high participation, and a high out-degree. By spending longer times in the main synchronization cluster and influencing several downstream nodes, they are likely to play a key role in the propagation of synchronization (Figure 5C). Numerous other nodes with lower out-degrees also consistently spend more time in the main cluster, indicating their higher susceptibility to the influence of the drivers, rather than themselves influencing other nodes. These driver areas include the perirhinal (PERI), entorhinal (ENTl), orbital (ORBI), reticular nucleus (RE), basolateral amygdala (BLA), piriform (PIR), and agranular insular area (Ai). In contrast, we observe some nodes that have very high out-degree, yet spend very little time in the synchronization cluster, and are thus likely not driver nodes. Notably, a majority of these “driver nodes” are located in cortical areas. The hindbrain spends the least time in the main cluster, consistent with our earlier findings.

Since we model resource dynamics explicitly, our model allows us to investigate the impact of resource dynamics on the propensity of tES for the network. In particular, we study the impact of the resource recovery-to-consumption rate ratio (α/β) by fixing the recovery rate (α = 0.01) and varying the consumption coefficient (β ∈ {0.01, 0.005, 0.002, 0.001, 0.0005, 0.00025, 0.0002}). For this set of parameters, we identify the range of resource bath sizes (λ_o) that support tES. We find that as the recovery-to-consumption rate ratio increases, both SWN and MBN exhibit a shift in the range boundaries to lower values, indicating a higher propensity for tES at low resource levels (Figure 6, panels A and B). The exponential decrease in range boundaries with an increase in the ratio suggests that even a slight change in the recovery-to-consumption ratio can significantly alter the network’s propensity to generate tES.

In contrast, as the recovery-to-consumption rate ratio increases, the width of the bistability region decreases. Moreover, within the bistable region also, the duration of time spent in the synchronized state decreases, resulting in a decreased probability of tES occurrence (Figure 6C).

Together, these results suggest an optimal intermediate ratio for a healthy brain as follows (Figure 6C): Assume that the brain always operates at a resource bath size such that it is near the bistable region (criticality hypothesis; Beggs, 2022; Kim et al., 2022). Moreover, the resource bath size also undergoes fluctuations. In this scenario, a very low ratio would mean a high probability of tES, that is, a high likelihood of seizures. If the ratio is very high, the width of the bistable region is so small that fluctuations in the resource bath size push the brain into the monostable hypersynchronized state, again increasing the likelihood of seizures. But an intermediate value of the ratio ensures a balance between these two extremes and should be observable in a healthy brain.

In this study, we have analyzed transient explosive synchronization dynamics in an abstract, though biologically inspired, and constrained model. Our model incorporates two key features that lead to transient explosive synchronization (tES) across diverse network topologies: resource constraint and adaptive coupling. In the model, energy resources form a pool with a maximum capacity, analogous to metabolic resources like ATP (Chhabria & Chakravarthy, 2016), and regenerate at a constant rate. Each node represents a mesoscopic brain region consisting of a neural population. While past models linked energy consumption to mean firing rate (Liu & Ching, 2017), our model ties it to local synchrony of the population, which in turn is correlated with the mean firing rate (Chawla et al., 1999). Additionally, adaptive coupling in our model is based on spike-timing-dependent plasticity (STDP) for populations of neurons, where connectivity changes are influenced by the level of synchrony between populations (Duchet et al., 2023).

With this model, we first show in an idealized SWN (and other topologies such as SFN) how adaptive coupling can give rise to resource-level-dependent hysteresis. Upon the addition of resource dynamics, this gives rise to tES for intermediate sizes of the resource bath. Recent research comparing diffusive and adaptive coupling, common modeling choices in networks of neural masses, has demonstrated a higher likelihood of networks with adaptive coupling to generate seizures (Lopes et al., 2023). Our findings reinforce this preference for adaptive coupling in exhibiting a higher tendency for seizure generation.

Our results hold well qualitatively when we apply the same model to a real biological neural network—the mesoscale mouse brain network from the Allen Brain Atlas (Figure 3). Although the structural network comes from healthy rather than epileptic mice, our results demonstrate the ability of the model to generate seizure-like dynamics in a biologically realistic network. The framework can be used to further study how perturbations to this network can increase their susceptibility to seizures, thus understanding the specific potential structural elements in diseased mice (or humans) that lead to epilepsy. The choice of a mouse brain as the model system is made because of the greater precision and completeness with which the anatomical connectivity can be measured, compared with noninvasive methods employed in humans.

We observe some very interesting deviations for the MBN compared with SWNs. The MBN reaches only a partially synchronized state, with specifically the hindbrain subnetwork never participating in the synchronization (Figure 3). Even within the synchronized cluster, we hypothesize the presence of a core that becomes fully synchronized, and drives the tES event, and a periphery that does not necessarily reach full synchronization. The energy levels of the nodes at the time of transition may provide a means to identify the core and the periphery. Moreover, the constitution of this core likely depends on the size of the resource bath, so that for different bath sizes, we observe different average transition energies (Figure 3C). We speculate that this added complexity is a result of the weighted, hierarchical, and modular network structure in real brains compared with our idealized SWN. These hypotheses and speculations are areas for further study to understand how the network structure affects its susceptibility to tES.

We then study the dynamics of synchronization propagation across the network at the level of communities (also known as subnetworks) and individual nodes. A salient feature we observe is that preceding each abrupt transition, there is a point of near-complete desynchronization across the network (Figure 4B, inset). This is accompanied by the formation of several small clusters that later merge into the main synchronization cluster. These phenomena are consistent with results obtained through mean field analyses (Zhang et al., 2014, 2015), as well as experimental observations at micro- and macroscopic levels (Aarabi et al., 2008; Jiruska et al., 2013; Mormann et al., 2003; Wendling et al., 2003).

At the intra-community level, we find that the cortico-thalamic networks (particularly the orbital, olfactory, and sensorimotor areas) exhibit a higher starting synchrony, a faster increase of synchrony, and a higher steady-state synchrony during the transitions, compared with subcortical communities (Figure 4C). This suggests their role in initiating and propagating the hypersynchronized state, consistent with extensive observations and predictions in literature (Holmes et al., 2004; H. Meeren et al., 2005; H. K. M. Meeren et al., 2002; Miao et al., 2014; Szaflarski et al., 2010; Zhong et al., 2022). Additionally, we find certain cortical networks (orbital, sensorimotor) to be more critical for synchronization propagation than others (visual) (Nersesyan et al., 2004).

Quantification of intermodular synchronization reveals that the synchronization expands in a hierarchical manner, as a propagating wave from the cortical to subcortical regions (Figure 4D), so that the subcortical areas synchronize with cortical areas before they synchronize among themselves. Although whole-brain recordings during generalized epilepsy are lacking, this would be an interesting hypothesis to test in model organisms with invasive electrophysiology.

We develop a novel algorithm to track the synchronization cluster lineages for individual nodes, which reveals the existence of a single large synchronization cluster during the transition, with several small clusters that dynamically join or leave it. Based on the consistency of time spent by the nodes in the main cluster, and their out-degrees, we find a set of driver nodes that hold the cluster together, irrespective of the initiating site. This driver set includes perirhinal, entorhinal, orbital, reticular nucleus, basolateral amygdala, piriform, and agranular insular areas. These predictions are supported by several experimental findings (Chen et al., 2003; Hamasaki et al., 2014): For instance, the entorhinal, perirhinal, and piriform cortex form a highly interconnected network with other limbic structures and have been shown to possess characteristics that make them susceptible to the initiation and spread of epileptic seizures (Vismer et al., 2015).

According to theoretical analysis (Zhang et al., 2015) and our cluster tracking results, abrupt transitions during tES are preceded by the formation of numerous small synchronization clusters. This is consistently preceded by almost complete desynchronization. The more of these clusters, the more abrupt the transition (Zhang et al., 2015). A similar phenomenon of synchronization cluster formation and interictal/preictal desynchronization is observed before critical transitions during seizures (Jiruska et al., 2010, 2013; Mormann et al., 2003). Experimental evidence shows that these individual clusters exhibit high-frequency oscillations (HFOs) of 80–500 Hz (Jiruska et al., 2010). These observations suggest that the preictal/interictal dynamics of HFO may vary depending on the seizure class that exhibits preictal desynchronization. Testing this hypothesis is intriguing, as it could emphasize the importance of considering seizure type when using HFO as a biomarker of ictogenesis (Saggio et al., 2020).

Lastly, our mechanistic model highlights the importance of an intermediate resource recovery-to-consumption ratio, effectively balancing the heightened tES likelihood and the occurrence of monostable hypersynchronous activity. This implies an optimal recovery-to-consumption range where seizures are infrequent (Protachevicz et al., 2019), and a constant hypersynchronous state is improbable. Deviations from this range may trigger abnormal brain states, suggesting a testable hypothesis for the susceptibility to epileptic attacks in relation to ATP demand and oxygen consumption rates observed during ictal and interictal epileptiform activity (Muddapu et al., 2020; Schoknecht et al., 2017).

To summarize, our mesoscale network model for generalized epilepsy applied to a real biological brain network makes several predictions that are consistent with experimental data and more biologically realistic and complex models. The simplicity, coupled with the generality, of the model holds significant value for two key reasons: First, its simplicity allows for the simulation of large-scale brain networks without significant concerns about computational load. Second, it potentially enables the study of seizure dynamics in a wide range of whole-brain networks and could have applicability from a translational perspective. By identifying the propagation pattern during seizures, we can potentially identify strategies to halt the propagation. Therefore, the model and techniques developed here can be applied to connectome data from actual epileptic brains, with the hope of identifying the seizure onset site and its progression.

To evaluate the dynamics of the model with different network structures (small-world, scale-free, mouse brain), we perform two types of simulations: adiabatic progression and bifurcation diagram construction. In the adiabatic progression, we systematically increase or decrease the fixed resource Λ to observe the global synchrony parameter for the conventional adaptive coupling model with fixed resources. This allows us to determine the hysteresis region of the system. For the bifurcation diagram construction, we increase the resource bath size parameter λ_o and measure the global order at each time point. For all simulations, we use α = 0.01, β = 0.002 (unless otherwise specified) and the initial phases θ_i are distributed uniformly in the range [0, 2π). Equations are simulated using the Euler method with a step size of 0.05.

To construct hysteresis (bifurcation) diagram in MBN, unlike SWN, we run simulations for a duration of 2000 time units for each value of Λ / λ_o through adiabatic progression with ΔΛ = 0.02 / Δλ_o = 0.02.

To study the progression of ES in a weighted MBN, we conduct four separate runs with a fixed value of λ_o = 2.87, each using distinct initial conditions for (phase) θ and (frequency) ω. The simulations spanned a duration of 20000 time units. From these simulations, we extract a total of 52 transients by selecting segments of length 50-time units, starting from the “onset index” just before the abrupt transition. We run simulations with different initial conditions to account for the sensitivity of spatial spread of tES to initial conditions in complex networks (Supporting Information Figure S9).

The Synchronization cluster tracking algorithm (SCTA) performs two major tasks: i) finds the synchronization clusters at each time unit (Supporting Information Figure SM1), ii) tracks the temporal evolution of identified clusters across different time units (Supporting Information Figure SM2).

The SCTA aims to expand synchronization clusters within the network based on a given synchronization threshold. It follows three key steps:

1. Expansion around Central Nodes: The algorithm begins by expanding the synchronization clusters around the central nodes, which are selected from the previous time step. These central nodes act as trackers for the clusters across different time steps. Nodes with local synchronization exceeding a predetermined threshold become part of the cluster, which terminates with nodes that fall below the threshold (Supporting Information Figure SM1).

2. Expansion for Unassigned Nodes: Next, the algorithm expands the synchronization clusters for any nodes that have not yet been assigned to a cluster. This step ensures that all nodes are considered and included in appropriate clusters based on the synchronization threshold (Supporting Information Figure SM2).

3. Update Central Node List: Finally, once all nodes have been traversed or become part of some cluster, the algorithm updates the list of central nodes for each cluster. Central nodes are the top 5 nodes in each cluster with highest local synchrony. These central nodes are important reference points for the clusters, preserving their cluster membership over time, and acting as seed nodes for the expansion of clusters in the next time step (Supporting Information Figure SM2).

By executing the SCTA over a 50-time unit window for a specific transient, we obtain the cluster sizes as a function of time. The largest cluster is defined as the main synchronization cluster, and we measure the time spent by each node as part of the main cluster within the 50 time step window.

We are grateful to Dr. Anton Arkhipov from the Allen Institute, Dr. Tapan Kumar Gandhi from IIT Delhi, and Dr. Srinivasa Chakravarthy from IIT Madras for helpful discussions on the topic.

Supporting information for this article is available at https://doi.org/10.1162/netn_a_00379.

Avinash Kumar Ranjan: Conceptualization; Formal analysis; Investigation; Methodology; Validation; Visualization; Writing – original draft; Writing – review & editing. Saurabh R. Gandhi: Conceptualization; Methodology; Project administration; Resources; Supervision; Writing – review & editing.

Avinash Kumar Ranjan, IITD Institute Assistantship. Saurabh R. Gandhi, IITD Young Faculty Incentive Fellowship. Saurabh R. Gandhi, Science and Engineering Research Board (SERB), Award ID: SRG/2023/000595.

The code and parameters that have provided the results presented here are available at GitHub (https://github.com/csndl-iitd/tES_mesoscale_connectivity_model.git; Ranjan & Gandhi, 2024).

Explosive synchronization:

Rapid and abrupt transition of a complex network from a state of disorder to a state of synchronized activity.

First-order transitions:

Sudden, irreversible changes in network behavior, like explosive synchronization.

Frequency-degree coupling:

The rate at which something (node) is oscillating in the network is directly related to how many other things it is connected to (degree).

Adaptive coupling:

Network connection strength grows stronger for pairs of nodes in a network with larger phase coherence, mimicking neural plasticity.

Scale-free network:

Networks with power law degree distributions, where a few nodes have many connections.

Small-world network:

Networks with high local connections and short paths between nodes, such that any two nodes can reach each other in just a few steps.

Spike-timing-dependent plasticity (STDP):

Synaptic plasticity mechanism where connection strength changes based on temporal correlation of neural spikes.

Mesoscale mouse brain network:

Connectome of mouse brain at intermediate scale, between the macroscale (large brain regions) and the microscale (individual neurons and synapses).

Bifurcation diagram:

Graph illustrating changes in system behavior as its parameters vary, often showing sudden transitions.

Interictal:

Period between seizures, characterized by absence of seizure activity.

Preictal:

Period preceding a seizure, marked by changes in brain activity signaling an impending seizure event.

Ictogenesis:

Process leading to the onset of a seizure, involving the buildup of abnormal electrical activity in the brain.

Connectome:

A map of neural connections in the brain, detailing how different neural regions are connected to each other.

https://github.com/csndl-iitd/tES_mesoscale_connectivity_model.git