Symmetry’s Edge in Cortical Dynamics
Multiscale Dynamics of Ensemble Excitation and Inhibition

Abstract

Creating a quantitative theory for the cortex poses several challenges and raises numerous questions. For instance, what are the significant scales of the system? Are they micro, meso or macroscopic? What are the relevant interactions? Are they pairwise, higher order or mean-field? And what are the control parameters? Are they noisy, dissipative or emergent?

To tackle these issues, we suggest using an approach akin to what has transformed our understanding of the state of matter. This includes identifying invariances in the ensemble dynamics of various neuron functional classes, searching for order parameters that connect important degrees of freedom and distinguish macroscopic system states, and identifying broken symmetries in the order parameter space to comprehend the emerging laws when many neurons interact and coordinate their activation.

By utilizing multielectrode and multiscale neural recordings, we measure the scale-invariant balance between excitatory and inhibitory neurons at a population level, referred to as ensemble E/I balance. This differs from the input E/I balance typically studied at the single-neuron level, focusing instead on the collective behavior of large neural populations. We investigate a set of parameters that can assist us in differentiating between various functional system states (such as the wake/sleep cycle) and pinpointing broken symmetries that serve different information processing and memory functions. Furthermore, we identify broken symmetries that result in pathological states like seizures.

This study provides new insights into the multiscale dynamics of excitation and inhibition in cortical networks, advancing our understanding of the underlying principles governing neural computation and dysfunction.

This research provides a novel perspective on understanding cortical excitability, drawing parallels between the established principles of state of matter and collective neural dynamics. We delve into relevant interactions and control parameters of excitatory and inhibitory neuronal assemblies, using multielectrode neural recordings to capture a multiscale portrait of neural ensemble activities. This work presents innovative analytical tools such as “collapse curve”, “partition curve” methods, along with the introduction of multiscale descriptors, providing comprehensive insights into the functional system states. This study underscores the impact of broken symmetries on information processing, memory consolidation, and association creation, while outlining their role in pathological states like seizures. Through a careful examination of cortical dynamics, we advance a quantitative theory for the cortex, underlining the importance of symmetry and scale-invariance in maintaining ensemble excitatory-inhibitory balance. We draw attention to the apparent dichotomy of simple underlying organizing principles and the resultant complex neural activity patterns. This study represents a significant stride towards understanding the cortical operational principles and dynamics, promising new insights and directions in the field of neurodynamics.

Introduction

Over the past two decades, advances in microelectronic fabrication techniques have made it possible to record the electrical activity of large populations of neurons both in vitro and in vivo. With the increasing size of these populations, researchers have shifted their focus from studying the responses of individual neurons to stimuli (such as in the seminal works of [1, 2]) to examining population dynamics and correlations between neurons [3, 4, 5], pairwise and/or higher-order correlations [6, 7, 8, 9, 10].

Despite the general agreement that information is represented and processed through correlated activity patterns in neuron populations, the existing descriptive statistical models often fall short in providing deep theoretical insights into the collective behavior of these populations. As a result, attempts to establish principles that recognize high variability in spiking patterns and spontaneous activity as fundamental theoretical constraints, typically adopt a top-down compartmentalization approach [11]. This approach follows Marr’s tri-level hypothesis [12, 13] and separates the problem into a computational (behavioral) level, an algorithmic level, and a biological implementation level. However, recent critiques argue that this framework is not always useful in the context of neuroscience [14], while others defend its relevance [15].

Alternatively, we need a quantitative computational framework that starts from the tri-level approach and adds further requirements such as hierarchical representation and learning [16]. However, while valuable, these approaches often fail to fully encompass the complexity of the brain’s functioning due to their assumption of separation of scales, which contradicts the interconnected and dynamic nature of complex systems [17, 18]. To accurately describe hierarchical representation in a way that reflects the cortex as an adaptive physical system, a successful theory must embrace the interaction between sub-assemblies and focus on the transfer of dynamics between observation scales without getting lost in the details [18, 19].

Given the challenges in modeling the activity of large neural populations, it becomes imperative to seek out identifiable signatures that can direct us in developing an effective theory. This theory should capture the macroscopic characteristics of the system, despite the high degrees of freedom at the microscopic scale of individual neurons. As the number of measured neuronal activities increases linearly, the number of interactions among neurons grows exponentially, complicating the task of obtaining a representative sample of the full distribution of activity patterns.

Traditional methods like correlation measures or dimensionality reduction of population dynamics often prove inadequate in this context. This is where the insights from physics become invaluable, particularly when we focus on identifying patterns and symmetries. These can reveal order parameters and broken symmetries, offering clues to the underlying laws governing the system [20, 21]. Instead of limiting our focus to cognition or animal behavior as metrics for describing the system, a more productive approach involves identifying invariances in ensemble dynamics, seeking out order parameters that link critical degrees of freedom to differentiate macroscopic states of the system, and exploring broken symmetries in the order parameter space to uncover the emergent laws when many neurons interact and orchestrate their activation.

In the context of neural dynamics, it is crucial to distinguish between different types of excitation-inhibition (E/I) balance. Some previous theoretical studies [22, 23], have focused on the “input” E/I balance, which pertains to the balance of excitatory and inhibitory inputs received by individual neurons within a network. This type of balance is essential for maintaining stable firing rates and preventing runaway excitation or inhibition at the single-neuron level. However, our study focuses on a different aspect: the “ensemble” E/I balance. This refers to the collective behavior of large populations of excitatory and inhibitory neurons. By aggregating the spiking activity of these neurons and normalizing by their respective numbers, we examine the macroscopic balance of excitation and inhibition across multiple temporal scales.

In our recent work, we have highlighted the invariance of ensemble E/I balance across different scales [24]. Our further investigations, however, also revealed that while a maximum entropy model of weak pairwise interactions (Ising model) effectively describes inhibitory neurons, it falls short in accurately capturing the spiking patterns of excitatory neurons [25]. The crucial role of diverse computing node classes is underscored when a disruption in ensemble E/I balance leads to seizures, a distinct pathological state characterized by diminished information processing capacity and heightened levels of macroscopic neuronal synchronicity [24]. Furthermore, traditional thermodynamic-based measures of population activity, such as heat capacity, have proven inadequate in differentiating macroscopic states like awake, light sleep, and deep sleep [25]. Here, we aim to expand on these findings by delving deeper into the multiscale properties, further investigating the invariance of ensemble E/I balance, searching for order parameters that distinguish different states of the system, and examining the interplay between microscopic and macroscopic aspects of pathological symmetry breaking during seizures.

Figure 1: Data. A. The implant location (left) and multielectrode array (right). The 10x10 array covers an area with the span of 4mm x 4mm. B. Sample recording from the multielectrode array. C. Multiscale representation of balanced Excitation and Inhibition Ensemble spiking Excitation-Inhibition balance. Sample ensemble spiking with different levels of coarse-graining are shown in the four color panels (Note, these panels are example subsets from the rastergram and LFP shown in B. Preservation of Excitation-Inhibition balance across scales shows mirrored activity for different lengths of times (scales). The negative sign is only used conventionally to represent the opposing nature of Inhibition (red) versus Excitation (blue). In each panel, the ensemble firing of a given cell category is normalized by the total firing power in the same category during the shown epoch. Each panel has 656 bins and the duration of each bin is shown in the middle panel. For example, the top right panel represents scale 30 and a duration of 1.0938 hr, showing 656 bin’s duration of 6003(sec). The mean (green) and standard deviations (magenta) of each epoch are shown with dashed horizontal lines, respectively. Vertical bars show the time window from which another coarser scale 656-bin epoch was chosen. The middle panel shows the even logarithmic spacing of scales, which was chosen for computational efficiency.

results

Understanding cortical neural activity presents the challenge of dealing with the multiple spatial and temporal scales inherent in neural dynamics. This multiscale neural activity exhibits behaviors that resemble state transitions in other physical systems with dynamics that transcend many different length scales [26]. In this study, we use the term ‘renormalization’ to describe the process of analyzing neural population activity at multiple scales by aggregating spiking data and normalizing by the number of neurons ($N$). While this approach differs from traditional renormalization in statistical physics, which involves rescaling interactions at different spatial scales, our method serves a similar purpose of identifying invariant properties across scales. By aggregating spikes from all excitatory neurons into one time series and normalizing by the number of excitatory neurons, and similarly for inhibitory neurons, we capture the multiscale dynamics of excitation and inhibition within the cortical network, providing insights into how these dynamics are preserved or altered across different temporal scales.

The renormalization group provides a powerful framework for studying symmetries and changes in a physical system as viewed at different scales [27, 28]. In the context of phase transitions, renormalization group theory is based on the idea that only slow modes are responsible for the phase transition, and phase transitions are closely tied to long-range behavior [29, 30]. The defining property of criticality, according to renormalization group theory, is that the characteristic length scale of the structure of the physical system, also known as the correlation length, becomes infinite [31]. In neuroscience, particularly in studying large-scale neural networks, analogous principles can be applied. Similarly, large ensembles of neural events provide data for describing long-range coupling in the temporal domain, which is analogous to spatial coupling in phase transitions. In neural systems, information is conveyed not just by the neural packets but also by their temporal order, leading to emergence of various collective states as the functional network undergoes changes. Thus, our study relies on multiscale dynamical features of ensemble E:I balance across different states of the wake/sleep cycle to better understand collective neural activity (see Fig.1.B,C).

.1 Multiscale Signatures of Balance Invariance

.1.1 Excitation-Inhibition Balance Across Temporal Scales

Recently, we have shown that across all different states of the wake-sleep cycle, excitatory and inhibitory ensembles are well balanced, and co-fluctuate with slight instantaneous deviations from perfect balance, mostly in slow-wave sleep. Remarkably, these correlated fluctuations are seen for many different temporal scales. The similarity of these computational features with a network model of self-generated balanced states suggests that such balanced activity is essentially generated by recurrent activity in the local network and is not due to external inputs [24].

To further examine invariance, we tested balance across multiple temporal scales ranging from 1 millisecond to 10 seconds (see Fig.1.C inset). When normalized, the excitatory and inhibitory ensemble activities reveal a scaled mirroring of each other across these scales. This invariance was observed in both awake and different sleep states (i.e. slow-wave sleep, light sleep, and rapid-eye movement sleep). For example, a detailed analysis across scales (Fig.1.C) demonstrates that when normalized to total cell-type ensemble firing power, Excitation and Inhibition qualitatively track each other’s fluctuations. To investigate whether the reduction in variance is a straightforward consequence of the central limit theorem, we compared the observed scaling of variance reduction to the expected 1/(bin size) scaling. Our analysis revealed that while the variance reduction partially aligns with the central limit theorem, deviations from the expected scaling indicate additional underlying dynamics. Specifically, the observed scaling exponent suggests that factors beyond simple averaging contribute to the multiscale balance of excitatory and inhibitory activities.

In a prior study [24], to ensure the robustness of our findings and address the concern that the observed variance reduction might be a straightforward consequence of the central limit theorem, we compared our results to those obtained from a network model of interconnected excitatory and inhibitory neurons displaying self-generated balanced activity states (COBA model). This model, as detailed in [24], consists of a conductance-based network and has been shown to exhibit balanced mirrored activity across multiple scales. The COBA model demonstrates that the overall balance is preserved across scales, with instantaneous deviations from perfect balance, similar to the experimental data. These findings indicate that the multiscale balance observed in our data is not merely a result of central limit effects but is consistent with the dynamics of an intrinsic network model that exhibits self-sustained activity. Specifically, the COBA model’s behavior shows that such balanced activity across scales arises from recurrent interactions within the network rather than from simple averaging processes. Therefore, the observations here align with the behavior observed in the COBA model, reinforcing that the observed variance reduction in our experimental data reflects the complex dynamics of cortical networks rather than a trivial consequence of averaging.

Importantly, however, it is crucial to acknowledge that, given the limitations of current recording technologies, we are significantly subsampling the neural landscape. Therefore, our estimates of balance are approximations, reliant on the data available from a limited number of neurons. This caveat suggests that a complete understanding of the true scaling of balance in neural ensembles remains an aspirational goal for future technological advancements. A more extensive sampling could potentially reveal further complexities and nuances in the multiscale dynamics of excitation and inhibition balance.

.1.2 Scale-Invariance in Correlated Fluctuations of Excitation and Inhibition in the Cerebral Cortex

A common pattern we observe is that the ensemble E-I distribution’s standard deviation (and variance) converges around the fixed mean as we transition from finer to coarser scales (see Fig.1.C panels and Fig.2.A1). This convergence suggests a systematic reduction in the variability of the E-I balance across different temporal scales. Examining the variance around the mean (i.e., coefficient of variation) for the entire recordings, we observe a symmetric decrease of coefficient of variation across scales (Table1). This symmetry underlines the robustness of the E-I balance irrespective of the temporal scale, and is disrupted when we randomize and break the temporal relation between excitatory and inhibitory ensembles (Fig.2.A2,A3). The dispersion from the diagonal line of balance in Fig.2.A1 indicates the instantaneous predominance of either excitatory or inhibitory systems, maintaining a dynamic balance (the diagonal in Fig.2.A1 represents the perfect line of balance).

To further clarify, the Mean Absolute Deviation (MAD) reflects the average distance between each data point and the mean, indicating overall variability. The Skewness measures the symmetry of the distribution, with values close to zero indicating a symmetric distribution. The Slope of the Coefficient of Variation ($S_{CV}$) measures how the ratio of the standard deviation to the mean changes across scales, with values close to one indicating a stable variability relative to the mean across scales.

Table 1: Indices of dispersion. The Mean Absolute Deviation (MAD) indicates overall variability, Skewness measures the symmetry of the distribution, and the Slope of the Coefficient of Variation ($S_{CV}$) reflects the stability of variability across scales. These measures demonstrate the consistency of the E/I balance across different temporal scales.

	MAD	Skewness	$S_{CV}$
AWAKE	$0.0023\pm 0.0030$	$-0.013\pm 0.046$
REM	$0.0028\pm 0.0032$	$0.052\pm 0.015$	$|S_{CV}-1|<\epsilon$
LS	$0.0021\pm 0.0024\$	$-0.010\pm 0.015$	$\epsilon<10^{-2}$
SWS	$0.0027\pm 0.0034$	$-0.012\pm 0.087$

Across sleep-wake cycle, this dynamic balance is consistent, as evidenced by the Mean absolute value deviation (MAD: $\frac{1}{n}\sum_{i=1}^{n}|x_{i}-\overline{X}|$,) for different states (awake, slow-wave sleep, light sleep, and REM) showing similar structure and statistical dispersion (see Table1). We measured the dispersion, representing the extent of dominance between Excitation and Inhibition, in a scatter plot of their activity ratios for each neuron pair. The dispersion decreased across all states as we increased the timescale of our analysis, indicating a more pronounced E:I balance at larger scales. However, the scatter plot’s structure remained consistent, even at higher timescales, suggesting the retention of E:I activity characteristics at coarser temporal scales.

In both slow-wave sleep and wakefulness, the distributions of the ensemble fractions of excitation and inhibition are symmetric in magnitude and frequency. The cross-subject average skewness for different states across the wake-sleep cycle were close to zero, indicating a near-perfect symmetry (see Table.1).This symmetry across scales reinforces the notion that the E:I balance is a fundamental aspect of cortical computation. The balance allows for a wide parameter space where input can be computationally processed while preserving balance, as suggested before [32]. Furthermore, this balance aligns with previous findings that dynamic interactions between excitatory and inhibitory conductances finely tune spike timing [33, 34]. Simulations have suggested that the (proximal vs distal) spatial distribution of inhibitory synaptic activity on pyramidal cells maintains the system in a state of inhibition-dominant excitability, awaiting excitation from the thalamus [35].

.1.3 Disrupting Temporal and Spatial Structure of Ensemble Spike Patterns through Randomization

We examined how randomizing the spike patterns of the neurons affected the multiscale balance of excitatory and inhibitory activity. Our approach involved two distinct randomization techniques: one that permuted the inter-spike intervals in the ensemble series, and one that circularly shifted the inhibitory and excitatory ensembles (for details, see the methods section). Both methods effectively disrupted the alignment of excitatory-inhibitory pairs along the diagonal line, which is essential for maintaining the balance between the two systems. Instead, randomizations transformed the data to a circular cloud, indicating a loss of the original structure and a shift to noise-dominated patterns. This resultant “sphering” is akin to a whitening” transformation because it converts the input vector into a white noise vector, with no remaining correlation or structure, as shown in Fig.2.A2,A3. It’s important to note that, unlike linear whitening transforms (such as Cholesky, PCA, or ZCA) which assume a multivariate normal distribution, nonlinear whitening transforms (implemented through neural networks, kernel whitening or ICA, or like those used in our study) can handle non-Gaussian distributions and remove higher-order dependencies through nonlinear transformation.

Each of these randomizations were designed to achieve a specific decorrelation objective. The “random permutation” method involved calculating the inter-spike interval (ISI) for the pooled ensemble series of each functional category (excitatory or inhibitory), then randomly shuffling the ISI for that category, and reconstructing a new temporal order of ensemble spikes by cumulatively summing the shuffled ISI. This process changes the arrangement of spikes within the series while preserving their number and distribution. The “circular shift of spike ensemble” method involved calculating the inter-spike interval (ISI) of each unit’s spike series, then shifting the spikes of each unit by a random value between 1 and the maximum ISI of that unit. The randomized neurons were then aggregated to create the randomized ensemble series of each functional category. This process preserves the internal spike timing of each unit, but disrupts the timing between units within each category.

Figure 2: Randomization of ensemble activity. A1: Data. Scatterplot of normalized ensemble excitation and inhibition. Blue, red, green and black represent 4 sample scales. Paired E-I scatter as elipsoids with their major axis of variation along the diagonal (perfect symmetry). A2:Random permutation of ISI in the ensemble series. A3:Fixed-ISI circular shift of spikes. Ensemble Excitation and Inhibition show correlated and balanced structure that is disrupted by sphering randomization. For details of the randomization techniques, see methods. Fluctuation collapse. B1. Fine-scale collapse curve versus random surrogates. The inset shows the histogram distribution of fluctuation of ensemble activity at the finest scale. B2. Comparison of collapse curves across different scales. Note that at coarsest scale, surrogate collapse curve approaches the data collapse curve. Different random surrogates collapse onto the same curve but distinctively are different from data collapse curve. B3. Comparison of collapse curves for different states (SWS: slow-wave sleep, REM: Rapid Eye Movement, and Awake) across various scales. Dashed circles are visual guides to exemplify that more or less collapse curves of different states have similar behavior in a given scale.

The first method acts like a linear whitening by disrupting the temporal order of aggregate spikes, while the latter method is more akin to a nonlinear whitening as it decorrelates the ensemble through disrupting inter-unit spike timing, thus effectively removing complex higher-order correlations. Ultimately, both randomization techniques disrupt the temporal or spatial structure of the ensemble spike patterns of the neurons, demonstrating the critical role of precise temporal and spatial structuring in maintaining the balanced E-I relationship in neural ensembles. Figure 2A1-3 demonstrates that random permuting of spike times destroys the E/I balance, an expected outcome given the balance at short timescales. This result underscores the importance of temporal structure in maintaining E/I balance and highlights the necessity of preserving temporal correlations for accurate modeling.

.1.4 Limitations of Current Experimental Techniques

It is worth mentioning that the multielectrode array samples spiking activity from an area about the size of a cortical column using a 10x10 electrode grid. Although the density of neurons in a cortical column varies across species [36, 37], the Utah array and other current experimental techniques massively subsample from the neurons in a cortical column, estimated to range between 10,000 to 20,000 per column [38, 39]. As recording technologies advance, they will enable the exploration of excitatory-inhibitory balance with higher spatial resolution across many neighboring columns. Given the underlying functional anisotropy at different spatial scales, such as our prior observation of stereotypical directionality of local field potential traveling waves at the scale of a cortical column [40] and in theoretical models of propagating waves across much larger patches of cortex [41], the investigation of ensemble excitatory-inhibitory translational symmetry will become increasingly feasible. This exploration will guide us in developing a formalism that could link micro and macroscales.

.2 Universality of Balance

.2.1 Collapse Curve Analysis of Ensemble Excitation-Inhibition Fluctuations

We then investigated multiscale properties of fluctuations in cortical ensemble excitation-inhibition through rescaling. We constructed a collapse curve analysis to reveal the self-similar or self-affine properties of these fluctuations across different time scales and states of the wake-sleep cycle. A collapse curve can be obtained by rescaling the axes of different curves of a physical quantity measured at various control parameter values, making them collapse onto a single curve, indicating universality or self-similarity [42]. This method has been shown to effectively demonstrate scaling behaviors in systems exhibiting critical phenomena, such as avalanches and crackling noise [43, 44]. We applied this method to the observed fluctuations of ensemble activity, represented by $f_{obs}=Ens_{E}-Ens_{I}$, where $Ens_{E}$ and $Ens_{I}$ are the ensemble activities of excitatory and inhibitory neurons, respectively. After sorting $f_{obs}$ in ascending order and summing all values, a symmetric parabolic curve emerged, reflecting the overall balance of ensemble excitation and inhibition (Fig.2.B).

To assess self-similarity or self-affine characteristics, we rescaled each curve for different time scales so that it would be bounded between 0 and 1 on the x-axis and then normalized the y-axis by the product of the number of bins. This normalization ensures a fair comparison of timeseries of varied lengths. We systematically identified non-collapsible regions for several comparisons: i) data versus random series, ii) different scales, and iii) wake-sleep states across different scales. The surrogates from various randomization methods showed similar behaviors and failed to collapse onto the original data curve (Fig.2.B1). This highlights that the observed E-I balance isn’t just a statistical byproduct but relates to an emergent collective pattern in ensemble spiking, consistent with shape scaling of avalanches in systems exhibiting crackling noise [44]. At the coarsest scale, the surrogate collapse curve approached the data curve (Fig.2.B2), suggesting that while the E-I balance persists at long timescales, the content might lose its informational structure. Comparing different states, we observed that various states more or less collapsed onto a similar regime (Fig.2.B3), reinforcing the notion of an invariant E-I balance across scales.

On a corollary note, it’s important to highlight that scaling ansatz can help to reveal the universal features of a system and to derive scaling relations among different observables. The ability of different curves of a physical quantity, measured under varying control parameter values, to collapse onto a single curve through rescaling, exemplifies the system’s scaling behavior near a critical point [43, 42]. This has been observed in neuronal avalanche dynamics, where universal critical dynamics are evident in high-resolution data [43]. The scaling collapse we observed intersects with other collective behaviors, particularly neural avalanches. Neural avalanches, characterized by spatiotemporal neuronal activity patterns without a specific size or duration, are often seen as indicators of scale-free dynamics and criticality in the brain. Criticality, a state where the system is poised at the edge of a phase transition between order and disorder, is believed to endow the brain with optimal computational capabilities. However, the nature of neural avalanches in relation to criticality is under debate when considering other factors such as network topology, synaptic plasticity, or external inputs.

Through rigorous analysis, we’ve shown that neural avalanches don’t align with self-organized criticality (SOC) [45], which posits systems naturally tuning to criticality without external parameters. SOC models seem inadequate in explaining the role of inhibition and the variety of complex spatiotemporal patterns observed in neuronal networks. Similarly, maximum entropy models fall short in capturing the complex dynamics of excitation and inhibition [25]. The scaling collapse, combined with our prior observations on Excitatory and Inhibitory populations, points towards a distinct form of criticality in cortical networks, characterized by scale-invariance across a broad range of parameters [42]. This could stem from multiple interacting elements or feedback mechanisms in the networks, potentially leading to various phases or coexisting phases near criticality.

.2.2 Partition Curve Analysis of Neuronal Firing Renewal

The Poisson process, a common model to describe neuronal firing as a continuous-time Markov process with exponentially distributed inter-arrival times [46, 47], falls short in accurately depicting realistic interspike interval distributions. This is primarily due to its inability to account for neuronal refractoriness, a state during which a neuron is unable to fire another spike [48, 49, 50]. This refractory period can result in gaps in the ensemble spike timeseries, complicating the accurate modeling of ensemble behavior. Therefore, taking into account refractoriness is crucial for a more precise representation of neuronal dynamics. To examine the temporal variation of the probability of an event occurring at a given time, we need new methods.

Refractoriness, particularly during in vivo high conductance states, significantly impacts the variability of neuronal responses [51]. To incorporate refractoriness and the resulting response variability in the stochastic description of spike trains, methods like filtering, renewal, or Hawkes processes are employed [52]. Each of these methods have certain drawbacks; for example, filtering can remove short-term fluctuations but also introduces bias.

Renewal theory, extending the Poisson process with arbitrary holding times, is an idealized model for random events under the assumption of independent, identically distributed times between events. It considers the memory of only the last event, useful for modeling ensemble spiking patterns with time since the last event in mind. However, renewal processes might not always be suitable for nonstationary data with trends or patterns violating the independence assumption [53, 52, 54].

Spike trains in typical neurons in vivo exhibit irregular behavior with considerable variability in interspike intervals and spike counts [55]. The source of this irregularity, whether thermal noise, microscopic chaos, or complex neural coding, remains debated [56, 57]. In neurons with strong adaptation, the interspike intervals are not independent, challenging the adequacy of renewal theory. Time-dependent generalizations of renewal theory or Hawkes processes, which model refractoriness and are adaptable to both stationary and non-stationary data, may be more suitable. However, Hawkes processes face limitations with time-varying intensity functions in non-stationary data [58]. These limitations highlight the need for new methods to accurately model the temporal variation of the probability of an event occurring at a given time.

We present a method to analyze the distribution of ensemble activity in a given functional category using a partition curve. The partition curve, a graphical representation of ensemble spike firing rate distribution, plots cumulative percentages of a variable against the population’s cumulative percentages. It indicates the percentage of ensemble activity of a specific size or certain size interval, providing insights into the distribution’s variability equivalence or unequivalence. Analysis of the distribution of ensemble activity using a partition curve can help us to gain a deeper understanding of the relationship between ensemble spiking patterns and renewal processes.

Figure 3.A illustrates a probability–probability plot comparing the distribution of a variable against a hypothetical uniform distribution. The partition curve, convex and lying below the diagonal line, represents the distribution’s variability. The diagonal line represents equivalence of normalized ensemble activity across time. The distance from the diagonal line indicates the degree of unequivalence in ensemble activity distribution. The closer the convex partition curve is to the diagonal line, the lower the ensemble activity unequivalence. The further away the partition curve is from the diagonal line, the greater the unequivalence in the distribution of ensemble activity.

The partition curve show that, notably, spike firing rates are not uniformly distributed, displaying significant variability. The partition curve’s utility extends to nonstationary data, allowing for temporal distribution comparisons. The partition function is a measure of the volume occupied by the system in phase space and can be helpful for identifying trends in the distribution of the variable and understanding how it is affected by changes in other variables. Color-coded curves in Fig.3.A show the partition curves for various scales (darker colors represent finer scales while brighter colors represent coarser scales). As coarse-graining increases, partition curves tend to approach the diagonal, yet they retain a distinct non-equivalence property. Our results show results show that spike firing rates are not evenly distributed and exhibit a large amount of variability in firing rate. This transition reflects the system’s functionally accessible microstates under specific coarse-graining. The black partition curve, representing a normal distribution, stands separate from the data-derived curves, indicating that the observed behavior is not attributable to random processes.

Figure 3: A:Ensemble activity partition curve. In each panel, partition curves for various scales are displayed against the diagonal (magenta) line of equivalence. A perfectly equal ensemble size distribution would be one in which every ensemble size has the same frequency (magenta line). The black curve represents the partition curve for the normal distribution surrogate series. Partition curves for different scales are color-coded, with brighter colors indicating coarser scales and darker colors indicating finer scales. With increasing coarse-graining, the partition curves approach closer to the diagonal and random partition curves but always maintain the characteristic non-equivalence property. B:Multidimensional features of Excitation and Inhibition co-occurrence across different states. In each vertical band, a distinct feature (extracted from the multidimensional co-occurrence matrix of Excitation/Inhibition/Scale) is shown for various states; Awake, REM (rapid-eye movement), LSL (light sleep) and SWS (slow-wave sleep) from left to right. The features belong to three main categories: contrast group (contrast, homogeneity and dissimilarity), descriptive group (correlation), and orderliness group (entropy and energy). Colors brighten as the distance of the matched pairs in the co-occurrence matrix increases. The radius of each circle shows the std. The features define the state-dependent properties of Excitation and Inhibition interaction in a multiscale computational framework. Note that each category of measures is normalized to the maximum across different scales.

.3 Symmetry Breaking

.3.1 Symmetry breaking in wake/sleep states

As we transition towards coarser scales, the joint distribution of variables simplifies. This reflects how macroscopic behaviors, observed at these scales, tend to be more universal than microscopic ones. The concept of a “fixed point of the renormalization group” implies a state of invariance under renormalization transformations, a crucial test for criticality beyond just thermodynamic analogies [59]. Simplification of systems on various length scales can be achieved by integrating out high-momentum degrees of freedom and rescaling (as shown here and previously in [24]). However, this simplicity at a macroscale should not be mistaken for mere reducibility to purely mechanistic microscopic phenomena.

In cortical computation, major functional states emerge from interactions among numerous excitatory and inhibitory neurons. Understanding these dynamics and information processing necessitates studying symmetry breaking, seen as deviations from a balanced excitation-inhibition state at the ensemble level. Both renormalization and symmetry breaking can be explored using multiscale co-occurrence feature analysis in EIS (excitation/inhibition/scale) space. To delve deeper into the joint probability of the fraction of Excitation and Inhibition across multiple temporal scales, we employed a multidimensional feature analysis. This feature analysis operates within the multidimensional space of excitation/inhibition/scale (EIS). Here, we computed a ‘Multiscale Normalized Co-occurrence Matrix’ (MNCM) for various distances and orientations, following the volumetric feature analysis [60, 61], adapted from textural analysis [62]. The MNCM is generated by calculating the frequency of occurrence of pairs of neighboring discrete points/pixels () within the multidimensional EIS space. The spatial relationship between these pixels is defined using an ‘Offset’ parameter, which relates to the spatial co-occurrence of Excitation versus Inhibition in the EIS. Each entry in the MNCM represents the frequency at which a pixel with value $i$ is adjacent to a pixel with value $j$. This spatial dependence matrix, therefore, provides a distribution of similar values (of E and I) at given offsets, essentially quantifying the co-occurrence of specific ensemble E and I values, considering the scale. We then utilized MNCM to assess various properties, including contrast (contrast, homogeneity, dissimilarity), descriptive measures (mean, variance, correlation), and orderliness (maximum probability, entropy, angular second moment or energy).

The analysis of the Multiscale Co-occurrence Matrix (MNCM) employs three distinct groups: Descriptive Statistics, Contrast, and Orderliness. Each group is tailored to specific aspects of MNCM. The Descriptive Statistics Group employs statistical metrics to delineate point intensity distributions at each scale, focusing on correlation measures between pairs. The Contrast Group assesses intensity variations within a scale, utilizing metrics like homogeneity, contrast, and dissimilarity. The Orderliness Group concentrates on scale regularity, measuring features like energy and entropy. These features are computed separately for every scale, then merged into a singular feature vector representing the entire multidimensional space.

Statistical analysis of feature vectors from the Multiscale Co-occurrence Matrix (MNCM) reveals notable state-dependent differences across Descriptive Statistics, Contrast, and Orderliness measures (see Fig. 3.B). For instance, ‘Multiscale Energy’ distinguishes REM from the awake state, and the “Contrast” measure shows maximum distance-dependence in light (LSL) and deep slow-wave sleep (SWS). ‘Homogeneity’ is more prominent in wakefulness, indicating uniformity of excitatory and inhibitory activity, while deep sleep shows less correlation with distance.

In the context of MNCM, high contrast indicates a large amount of variability or differences in the distribution of excitatory and inhibitory activity. If ‘Contrast’ is distance-dependent across domains in light sleep and deep sleep, this could suggest spatial symmetry breaking. Different times of sleep, like UP or DOWN states in SWS, may exhibit different balances of excitation and inhibition, perhaps due to differentially tuned localized neuronal dynamics or inputs. Note that distance-dependence here refers to the pair values in the EIS.

‘Energy’ and ‘Entropy’ show opposing features of the wake-sleep cycle. ‘Energy’, a measure of uniformity or orderliness, is higher in awake or REM, states which have more asynchronous activity. The distinction of REM from the awake state by the ‘Multiscale Energy’ might reflect a state-dependent shift in the balance of excitatory and inhibitory activity, a form of temporal symmetry breaking, that is unique to information processing in dream versus wakefulness. On the other hand, ‘Entropy’ is higher in SWS and LSL, with REM showing the least levels. This could be due to the various dissociation of cortical events from thalamic events during non-awake states. ‘Homogeneity’ and ‘Correlation’ showed less distinctive variations across different states.

‘Homogeneity’, as a measure of the closeness of the distribution of elements in the MNCM to the MNCM diagonal, increases with the uniformity of excitatory and inhibitory activity. The prominence of ‘Homogeneity’ across multiple scales in wakefulness and its lesser correlation with distance in light and deep sleep may reflect a more detailed balance of excitation and inhibition during wakefulness and more localized dynamics (LSL) yet more phase undulation during sleep. This could be due to periods of intense neural activity and silence (UP and DOWN states) that characterize non-rapid eye movement (NREM) sleep.

The observed state-dependent feature differences in excitation-inhibition balance suggest variations in symmetry breaking at the ensemble level. These changes might be influenced by factors like neuronal firing rates, synaptic strengths, or neuromodulator levels during different wake-sleep states. Such distinctive multiscale properties could stem from intrinsic computational characteristics of various stochastic network states like asynchronous irregular (AI), synchronous regular (SR), and asynchronous regular (AR) regimes [63, 64]. Further research, integrating modeling and theoretical studies, is essential to elucidate the relationship between stochastic network states and the multiscale features of the excitatory-inhibitory balance discussed here.

The dynamic modulation of EI balance, and consequently the symmetry of the system, across different sleep states could serve as a mechanism for achieving diverse functional outcomes or maintaining homeostasis. These state-dependent changes in multiscale co-occurrence features may be indicative of these underlying processes. It’s worth noting that the differential Ising model fitting of Excitatory (E) and Inhibitory (I) neurons [25], could potentially manifest as the observed state-dependent differences in symmetry breaking. This reflects different sensitivity to changes in control parameters. Future investigations should include a more comprehensive analysis of spatial and temporal correlations, utilizing a larger dataset of neuron recordings, to test and refine this methodology. This approach will enhance our understanding of the system and contribute to the refinement of this method.

.3.2 Symmetry Breaking and Balance Reset in Pathological States: Insights from Seizures

Our recent study indicates that the balance of excitation and inhibition breaks down during electrographic seizures [24]. This imbalance is characterized by a notable shift in excitatory and inhibitory neuronal firing patterns, as well as significant fluctuations in the oscillatory contents of the local field potential (LFP)(Fig. 4.C1,C2). We observed an increase in high frequencies during a seizure in the average wavelet coherence, followed by a shift to lower frequencies, aligning with established frequency features of seizures in ECOG (electrocorticogram) and LFP recordings [65, 66, 67, 68].

LFP comprises multiple oscillatory patterns simultaneously; simply focusing on the wavelet modulus’ maxima doesn’t truly encapsulate the interaction of oscillations and their alterations during seizure. Wavelet coherence helps in establishing the relationship between LFP’s frequency structure and seizure events [69]. Utilizing the “multiridge” detection approach , we identified ridges of strong components across various frequency ranges, indicating abrupt frequency shifts during seizures. Fig. 4.A illustrates these multiscale LFP coherence ridges, showing a sharp transition concurrent with major fluctuations in ensemble excitatory and inhibitory activity, marking the onset of a seizure.

Figure 4: A: Local field potential (LFP) multiscale during seizure. A1. Multiscale coherence ridges. Crazy-climber algorithm, a relaxation method based on MCMC, was used to extract multiple ridges from the energetic distribution in LFP wavelet coherence. Sharp transition in activity across scales (frequency ranges) shows the emergence of seizure. A2. Multiscale heatmap of synaptic current. Synaptic current (estimated from the normalized ensemble spiking) also shows a multiscale signature of pathological symmetry breaking (during seizure) at the same time that multiridges of wavelet coherence show a sharp transition from low to high frequencies. B: Symmetry breaking in pathological state (seizure). B1. Rasterplot of excitatory (blue) and inhibitory (red) unit activity during a 9 minute epoch, B2. zoom in to the middle 80 sec. B3. Reset of balanced activity after seizure. Normalized cumulative activity of the signed ensembles (in green) shows a deviation from the estimated balanced trajectory (red) at the onset of seizure. Dark and light gray lines show Monte Carlo simulation of trend stationary and difference stationary processes. The normalized cumulative activity shows a return to the balance trajectory toward the end of the epoch.

Local field potentials (LFPs) are slow extracellular electric potential fluctuations in the brain, reflecting the combined activity of neuronal populations [70]. The LFPs are influenced by the synchronous input and output activity of many neurons, including synaptic potentials and other slower transmembrane current flows [71]. However, the amplitude of the LFP is not a simple linear sum of synaptic activity. Due to frequency-dependent filtering by the extracellular space, the LFP is a filtered version of synaptic activity, with higher frequencies attenuated more than lower ones [72, 73, 74]. This has implications for interpreting LFP measurements, as the LFP does not directly reflect raw synaptic activity.

To address this issue, we examined the synaptic current estimated from the ensemble inhibition and excitation . This estimated synaptic current represents a measure of the balance of excitation and inhibition in the neuronal population, with positive values indicating a net excitatory influence and negative values indicating a net inhibitory influence. A heatmap of this synaptic current, estimated from normalized ensemble spiking, reveals a multiscale signature of pathological symmetry breaking during seizure (Fig.4.B). Two distinct bands of maximum difference, corresponding to frequency perturbations observed in LFP oscillations (Fig.4.A), delineate the temporal zone of the seizure. The initial band highlights a clear breakdown of balance, followed by another significant fluctuation, suggesting a rebalancing attempt.

The observed breakdown of excitation-inhibition balance during seizures is followed by a significant phenomenon. A raster plot of neuronal activity over a 540-second epoch reveals the multiscale nature of this imbalance leading to a seizure (Fig.4.C1,C2). This “multiscale” term refers to variations visible across different time (fine to coarse-grained scales) and spatial scales (spiking and LFP). Post-seizure, there’s an apparent reset in the multiscale balance, as indicated by the return of normalized cumulative activity to its balanced trajectory. The extent of this imbalance is marked by the actual behavior’s (red line) deviation from a projected balanced path (green line). Comparisons with simulated time series, representing stochastic (difference-stationary process) and deterministic (trend-stationary) processes, suggest unique adaptive features in the brain facilitating this return to balance. Observed in a limited number of seizure recordings, further research with more data is necessary to confirm these properties and explore their presence in other seizure types.

Discussion

The macroscopic behaviors of a system, observed at larger scales, can often be described by simpler and more universal principles than the microscopic mechanisms that govern the behavior of individual components of the system [75]. Neural field models and neural mass models represent different methodologies for modeling neuronal behavior. Neural field models consider the spatial distribution of neuron groups in a coarse-grained area [76, 77], and explain how neural activity variables change over space and time [78]. Neural mass models, in contrast, focus solely on tracking neural activity over time [79, 80] and do not consider spatial aspects. Both models employ principles of population coding and time-graining [81, 82]. However, mean-field models, typically used in these approaches, are limited. They fail to account for interactions across different scales and assume Gaussian fluctuations, focusing mainly on linear responses to stimuli. Such models overlook crucial nonlinear interactions vital at macroscopic scales. For a more comprehensive understanding of brain dynamics and its various types of criticality, advanced methods that extend beyond traditional mean-field models are necessary [83].

The renormalization group method addresses the limitations of mean-field models in understanding complex brain dynamics. This method integrates out high-momentum degrees of freedom, often associated with fast, microscopic details of the system, and rescales the remaining ones to focus on larger-scale or slower aspects [84, 85, 86, 87]. The invariance of probability distributions under repeated coarse-graining, a fixed point of the renormalization group, is a key test for criticality [59]. This invariance suggests that the system might be in a critical state, indicating the broader applicability of the renormalization group in understanding complex systems like cortex. An example of this is the phenomenological renormalization group (PRG) analysis used on single-neuron recordings from a mouse’s hippocampus, which revealed scaling characteristics in supercritical regimes for finite systems [88].

In our study, we adopted a coarse-graining approach akin to renormalization to analyze neural population activity at multiple temporal scales by aggregating spikes from all excitatory neurons into one time series and normalizing by the number of neurons (N_e), and similarly for inhibitory neurons. This method allows us to identify invariant properties across scales, similar to the goals of renormalization in statistical physics. We then conducted several tests to evaluate the system’s invariance properties. Our findings revealed a consistent decrease in the variance of ensemble Excitatory-Inhibitory (E-I) as the system transitions from finer to coarser scales, highlighting scale-invariance in the correlated fluctuations of excitation and inhibition in the cerebral cortex. Throughout the sleep-wake cycle, we observed a balance between excitatory and inhibitory activity, which was preserved across different scales and states. This E-I balance, which is crucial to cortical computation, remained intact even as the distribution narrowed with coarser temporal scales.

In the context of our study, ‘symmetry breaking’ refers to the transient disruptions in the balance between excitatory and inhibitory activities at the ensemble level. The ‘symmetry’ in this case is the near-equilibrium state of ensemble E/I balance, which is temporarily broken during periods of high neural activity. These transient imbalances are quickly restored, maintaining overall network stability while allowing for flexible information processing.

The novelty lies in framing this balance within the context of symmetry breaking, or symmetry’s edge, underscoring the importance of transient disruptions in this non-static balance for information transmission and computational processes. Unlike static systems where symmetry breaking leads to a new equilibrium, the cortical network dynamically restores the balance after transient unbalances, allowing for flexible and adaptive information processing.

This framing provides several new insights:

1. Dynamic Adaptability: It emphasizes the brain’s ability to quickly restore balance after disruptions, which is crucial for maintaining robust information processing capabilities.

2. Multiscale Balance: It highlights the importance of maintaining E-I balance across multiple scales, ensuring stability and functionality at both micro and macro levels.

3. Pathological States: Understanding symmetry breaking helps elucidate the mechanisms underlying pathological states, such as seizures, where this balance is disrupted.

By examining E-I balance through the lens of symmetry breaking, we can gain valuable insights into the dynamics of cortical processing, linking principles from physics and neuroscience to better understand how the brain maintains its functional stability and adaptability.

It is important to distinguish between the ensemble E/I balance studied here and the input E/I balance discussed in previous research. Input E/I balance, as investigated by [22, 23], focuses on the balance of excitatory and inhibitory inputs received by individual neurons. This type of balance is crucial for maintaining stable firing rates at the single-neuron level. In contrast, our study examines the ensemble E/I balance, which pertains to the collective activity of large populations of excitatory and inhibitory neurons. By aggregating and normalizing their spiking activities, we gain insights into the macroscopic balance of excitation and inhibition across different temporal scales.

Previous studies have demonstrated that the timing of spikes is finely adjusted by dynamic interactions between excitatory and inhibitory conductances [33, 89]. This balance helps regulate the timing of action potentials [34, 90], which is vital for various neural processes. Simulations and experimental studies suggest that the spatial distribution of inhibitory synaptic activity on pyramidal cells primarily controls the system’s activity, keeping it in an excitable state, ready for the next wave of thalamic excitation [35, 91]. In essence, the balance of excitatory and inhibitory influences prepares the cortex to respond swiftly and appropriately to new incoming information, indicating that this balance is a key feature of cortical computation. It’s worth noting that there is a broad parameter space where balance can be maintained while computational processing of inputs takes place [32, 92].

We utilized collapse curve analysis to examine the multiscale attributes of cortical ensemble excitation-inhibition balance fluctuations throughout various wake-sleep cycle stages, offering insights into the system’s self-similarity or universality. Various randomization methods resulted in surrogates that displayed similar behavior but did not align with the actual data curve, underscoring that the observed balance is not a random statistical characteristic but is an emergent property of collective ensemble spiking patterns. This concept is consistent with the notion of emergent phenomena in complex systems, as discussed in the study of crackling noise, where patterns emerge from collective interactions rather than random processes [44].

At the coarsest scale of analysis, the surrogate collapse curve began to align with the data collapse curve, implying that while the balance is maintained, the content of the information might lose significant structure at these extended timescales. This observation is in line with findings in neuronal avalanche dynamics, where scaling behaviors suggest that information content can diminish at larger scales, reflecting a shift from detailed interactions to more coarse-grained descriptions [43].

Furthermore, different states seemed to converge onto a similar regime, underscoring a consistent balance across scales. This exploration of scale-invariant behavior within cortical networks provides a new viewpoint on criticality, suggesting that the system exhibits this behavior across a vast range of parameters. This finding resonates with the concept of universality in critical phenomena, where systems exhibit similar scaling laws across different conditions and parameters, as shown in the study of avalanche shapes and scaling in physical systems [42]. This suggests that the observed criticality in cortical networks may be driven by competing interactions or feedback mechanisms, rather than being confined to a single critical point.

This aligns with our previous suggestions that the balance could be due to multiple interacting factors or feedback mechanisms in the cortical networks [45, 25]. Our methodology and findings draw parallels with studies on crackling noise and avalanche dynamics, further emphasizing the universality of the observed scaling behaviors in neural systems. By comparing our results with established models in other complex systems, we can better understand the underlying dynamics and the nature of criticality in cortical networks.

We have developed a method known as the “partition curve” to analyze the distribution of ensemble activity within groups of neurons. This curve provides a means to evaluate the distribution of ensemble activity based on size or intervals. Importantly, the curve can illuminate the degree of equality, or the absence thereof, in the distribution of ensemble activity. This tool could assist in the investigation of spiking patterns within neuron ensembles. A significant finding from our use of the partition curve was the considerable variability in spike firing rates, which suggests a distribution that is not uniform.

Partition curve is well-suited to track temporal changes in non-stationary data. Further analysis comparing our findings with a normal distribution supported the conclusion that the patterns observed were not simply the result of random processes. This implies that the patterns and distributions identified in the neuronal ensemble activity are potentially meaningful and could provide valuable insights into how these neurons function as a group. At coarser scales, the partition curves appeared to converge towards the diagonal and random partition curves, whilst preserving a unique non-equivalence attribute. This shift signifies the spectrum of microstates - or short-term, specific states of activity - that the system can access within an ensemble at distinct coarse-graining levels.

In statistical mechanics, a degree of freedom is a single scalar number describing the microstate of a system. The specification of all microstates of a system is a point in the system’s phase space. Entropy, in both thermodynamics and information theory, is a measure of the number of possible microstates or configurations that a system can have. In this sense, there is a connection between degrees of freedom, states of the system, and information (in terms of entropy). The more degrees of freedom a system has, the more possible microstates it can have, and the higher its entropy will be. This means that system at finer scales, with more degrees of freedom, can contain more information.

In our study, we introduced multiscale descriptors to analyze degrees of freedom, rescaling, and the simplification of joint distributions under coarse-graining. These descriptors, derived from the Multiscale Normalized Co-occurrence Matrix (MNCM), capture various aspects of the excitation/inhibition/scale multidimensional landscape. Specific attributes, like correlation, contrast, homogeneity, energy, and entropy, offer unique insights into the system’s behavior. These features reduce the complexity of the multiscale data to a set of numbers that describe key aspects of its structure. The MNCM analysis, applied across multiple temporal scales, highlights collective behavior and quantifies the system’s degrees of freedom during rescaling.

In particular, the “Multiscale Energy” metric, which signifies uniformity, distinguishes REM sleep from wakefulness. This suggests that there are unique shifts in the balance of neural activity between these states. The “Contrast” metric, which quantifies variability, demonstrates a strong distance-dependency in both light and deep sleep, hinting at the occurrence of spatial symmetry breaking at the ensemble level. This symmetry breaking is further explored with the “Homogeneity” metric, indicative of uniform activity. It is more pronounced during wakefulness and exhibits less distance-correlation during deep sleep, suggesting the presence of complex neural balances and localized dynamics in the sleep-wake cycle.

The ‘“Entropy” metric, which measures disorder, is highest during deep and light sleep and lowest during REM sleep. This reflects the dissociation of cortical and thalamic events in non-awake states. Symmetry breaking might manifest as a transition from balanced, coordinated activity to a state where certain regimes are favored over others at the ensemble level. The Multiscale Network Community Metric (MNCM) can effectively capture such transitions by quantifying the joint probabilities of different patterns of excitation and inhibition, thereby discerning state-dependent differences in neural activity. This approach could provide valuable insights into the collective behavior of neurons.

.4 Concluding Remarks

Philip Anderson, in his seminal paper More is different, posited that the internal structure of matter doesn’t necessarily mirror its total state’s symmetry and that large systems might not exhibit the symmetry of their governing laws [93]. He further suggested that temporal ordering and periodicity, alongside broken symmetry, are pivotal in biological systems. In our study, we propose that the nervous system occupies a distinct position in computational state space. While the global excitatory-inhibitory balance is scale-invariant and consistent across states (wakefulness, REM, NREM), momentary dominance shifts between excitatory or inhibitory neurons are more pronounced at finer scales and vary with the state. Yet, these symmetry-breaking events, characterized by temporary dominance of either excitatory or inhibitory neurons, are counterbalanced, preserving a scale-invariant overall balance within the ensemble of excitatory-inhibitory neurons. At coarser scales, the possibility of momentary symmetry breaking decreases, indicating a more rigid system. Symmetry breaking, which is fundamental to information [94, 95], plays a crucial role in the cortex. Specifically, momentary symmetry breaking is employed for information processing during wakefulness, memory consolidation during slow-wave sleep (SWS), and the creation of new associations during REM sleep. However, in contrast to systems like crystals or magnets, where symmetry breaking leads to a minimum energy state, the neocortex maintains its information processing capabilities by rebounding from these events, thus sacrificing a minimum energy state for functional flexibility.

The dynamic properties of ensemble excitatory-inhibitory balance in the nervous system show parallels with the concept of symmetry breaking in dissipative structures, as discussed by Prigogine [96, 97]. Previously, the link between the nature of dissipative non-equilibrium and symmetries in dynamical systems prompted considerations of symmetry breaking in large-scale neural activity [98].

This heuristic model of visual hallucination, an extension of [99], despite having certain biologically implausible features such as a homogeneous planar sheet of neurons with radial isotropic connections, is capable of producing a variety of spatial patterns due to the presence of both excitatory and inhibitory neurons and recurrent feedback. This occurs when the homogeneous resting state becomes unstable. A common thread between this model, our quantitative framework, and dissipative structures is the presence of symmetry breaking and unstable symmetry, which serve as their dynamical hallmarks.

In fact, the instability of symmetry in dissipative structures may be the reason why biological systems need to constantly harvest energy to return to the unstable symmetric state [100]. We speculate that perhaps the significant energy usage in cortex is to re-establish the unstable symmetry that is constantly breaking during information processing. Dissipative structures share many attributes with biological systems, such as multistability, sustained oscillations, and spatiotemporal patterns manifesting as propagating waves [101]. Our prior examinations of cortical field potentials show the presence of anisotropic traveling waves and dominant oscillatory events across the wake/sleep cycle [40]. These are signatures of a complex adaptive system with many computational length scales, lending itself to a form of invariant ensemble excitatory-inhibitory balance where momentary symmetry breaking at fine timescales provides the capacity to process information.

Seizures manifest as abnormal symmetry breaking, not confined to momentary excitatory or inhibitory dominance but disrupting ensemble balance across scales. Concurrently, there’s a notable shift in spiking energy in the estimated synaptic current coinciding with transient oscillatory changes in LFP, characterized by an abrupt rise in high frequencies and reduction in slow frequencies. This pathological symmetry breaking halts information processing, leading to consciousness loss. Afterward, a return to balanced state occurs, affirming scale-invariant balance in cortical dynamics.

Let us emphasize that the complexity of neural activity patterns, emerging from fundamental principles like scale-invariant balance, highlights the contrast between the simplicity of laws and their diverse outcomes. By quantifying the catalog of invariances across states and (spatial and temporal) scales, we can search for governing principles and diverse neural activity patterns that emerge from these simple laws.

METHODS

.5 Experiments

For this study, we used recordings from multielectrode array implants in layers II/III of the human temporal cortex. The electrodes in these arrays capture both local field potentials (LFPs) and spike trains of individual neurons. LFPs represent the aggregate electrical activity of neuronal populations, while spike trains are discrete events corresponding to individual neuron firings. The electrodes in these arrays are arranged in a 10x10 grid with $400\mu m$ spacing, covering an area of 4x4 millimeters (see Fig.1.A). The four corner electrodes are used for grounding, leaving a total of 96 active recording electrodes. Data is initially sampled at 30kHz, then filtered and thresholded for spike detection. We used a combination of spike-waveform characteristics and the sign of the short-latency cross-correlogram of spike time pairs to classify units as either excitatory (E) or inhibitory (I). Video, electrocorticogram (ECOG), and scalp EEG were used to classify the 12-hour recordings into state labels: Awake, Light Sleep (stages II and III), Deep Sleep (Slow-Wave Sleep), and Rapid Eye Movement (REM). For further details on recording methods and morpho-functional characterization of excitatory and inhibitory neurons, see the cited references [102, 103]. Patients gave their informed consent to the implantation procedure, which had been approved by the Massachusetts General Hospital Institutional Review Board (IRB) in accordance with the ethical standards of the Declaration of Helsinki.

.6 Renormalization: Spiking coarse-graining and field estimates

.6.1 Ensemble Excitation and Inhibition.

If we define each neuron as an independent particle and its stochastic spiking as random fluctuations in space, the observed randomness can be absorbed by redefining a cluster of neurons (instead of a single neuron) to account for self-interactions, pairwise interactions, and higher-order interactions. This ensemble is then temporally coarse-grained by summing all spikes in a given time-bin, with scales ranging from 1 millisecond to 10 seconds (equally spaced in a logarithmic fashion) to specify the scale of observation. The logarithmic scale is chosen to reduce computational burden, with a denser spread of scales at finer time resolutions.

While these recordings yield a high number of neurons from each array, we face both massive subsampling (in our case, hundreds of neurons out of tens of thousands in the columnar region of recording) and spatial nonuniformity. Although we had a 4/1 match of the number of excitatory and inhibitory neurons with anatomical data and a 4/1 higher spiking rate of inhibitory to excitatory neurons [104, 102], given the subsampling and nonuniformity, the ensemble series are normalized by the number of neurons. This procedure was done independently for excitatory and inhibitory cell categories to yield ensemble fractions of excitation and inhibition.

In this study, we use an approach akin to ‘renormalization’ to describe the process of analyzing neural population activity at multiple scales by aggregating spiking data and normalizing by the number of neurons ($N$). While this approach differs from traditional renormalization in statistical physics, which involves rescaling interactions at different spatial scales, our method serves a similar purpose of identifying invariant properties across scales. By aggregating spikes from all excitatory neurons into one time series and normalizing by the number of excitatory neurons, and similarly for inhibitory neurons, we capture the multiscale dynamics of excitation and inhibition within the cortical network, providing insights into how these dynamics are preserved or altered across different temporal scales.

It is important to note that this “ensemble” E/I balance differs from the “input” E/I balance discussed in some prior theoretical studies [22, 23]. While input E/I balance focuses on the balance between excitatory and inhibitory inputs received by individual neurons, our study examines the collective behavior of populations of neurons, providing a macroscopic view of excitation and inhibition dynamics.

.6.2 Synaptic current.

To explain the multiscale dynamics of spiking excitation and inhibition with multiscale representation of local field potentials (LFPs) using their wavelet transforms, we used unit activity to reconstruct estimated synaptic current. It is suggested that synaptic current generates LFPs [105, 106], which reflect a field measure of the tens of thousands of neurons present in the recording area from which spikes are sampled. The modeling involved convolving each spike with an exponential kernel [107].

$$C(t)=\int_{-\infty}^{\infty}D({t}^{\prime})exp[-(t-{t}^{\prime}/\tau_{s})]d{t}^{\prime}$$ (1)

Here:

$C(t)$ is the estimated synaptic current, $D(t)$ represents the spike trains, $\tau_{s}$ is the time constant for the synaptic current.

Different time constants ($\tau_{s}$) are used for the convolution depending on whether the spike is excitatory (3 milliseconds) or inhibitory (10 milliseconds). This is based on prior studies that found differences in the time constants for excitatory and inhibitory synaptic currents [34]. Synaptic current at any time is a weighted sum (or integral, in the continuous case) of the spike train, where the weights decrease exponentially as the time difference between the spike and the current time increases. After convolving the spike trains with their respective kernels, the estimated synaptic current is obtained by subtracting the ensemble of inhibition from the ensemble of excitation. This result is then corrected by the number of cells in each category (excitatory or inhibitory).

.7 Randomization of ensemble activity

.7.1 Random permutation.

We used this test to verify that randomizing the aggregate spike series by itself cannot mimic the observed excitatory-inhibitory balance. After calculating the inter-spike interval (ISI) for the pooled ensemble series of in a given functional category (Excitatory or Inhibitory), a random permutation of the ISI was performed for that functional category ensemble series. This was followed by a cumulative summation of the ISI to create a new temporal order of ensemble spikes. This process ensures that the randomized ensemble series has the same number of spikes and set of ISIs as the original ensemble series, but with a different arrangement of spikes within the series. The observations verify that the excitatory and inhibitory balance is not reproducible by randomization of ensemble spikes and the observed fluctuations $f_{obs}=Ens_{E}-Ens_{I}$ are not results of random events such as $f_{randperm}=SurrogateEns_{E}-SurrogateEns_{I}$.

.7.2 Fixed-ISI circular shift of spike ensemble.

Before creating the ensemble series, the inter-spike interval (ISI) of each unit’s spike series was calculated (by “unit” here, we refer to a single neuron). The spikes of each unit were then shifted based on a random value between 1 and the maximum ISI of that unit’s spike series. The randomized units were then aggregated to create the randomized ensemble series of a given functional category (Excitation or Inhibition). This process ensures that the resulting ensemble series is made up of units with intact internal spike timing but with disrupted timing between units in that functional category. The results show that by destroying the relation between ensemble inhibition and ensemble excitation while preserving their internal structure, the observed fluctuations and most importantly, the tightly bound relation of Excitation and Inhibition is lost.

.8 Multiscale collapse of ensemble activity

The collapse curve analysis used in this study is inspired by techniques employed in the study of avalanches and crackling noise. These techniques involve rescaling data from different conditions or scales to reveal universal scaling behaviors [44]. Specifically, the collapse curve is obtained by rescaling the axes of different curves of a physical quantity measured at various control parameter values, making them collapse onto a single universal curve, indicating self-similarity or self-affine properties [42]. This approach has been demonstrated to reveal critical scaling behaviors in diverse systems, highlighting the universality of such phenomena [43].

For our analysis, we observed fluctuations in ensemble activity, represented by $f_{obs}=Ens_{E}-Ens_{I}$, where $Ens_{E}$ and $Ens_{I}$ are the ensemble activities of excitatory and inhibitory neurons, respectively. After sorting $f_{obs}$ in ascending order and performing a cumulative summation, we obtained a symmetric parabolic curve, reflecting the overall balance of ensemble excitation and inhibition. We then rescaled each curve for different time scales so that it was bounded between 0 and 1 on the x-axis and normalized the y-axis by the product of the number of bins. This normalization ensures a fair comparison of timeseries of varied lengths. We systematically identified non-collapsible regions for comparisons between: i) data versus random series, ii) different scales, and iii) wake-sleep states across different scales.

The general form of collapse curve for ensemble excitation and inhibition can be written as:

$$C(t)=\frac{F(t)-F_{min}}{F_{max}-F_{min}}$$ (2)

, where $C(t)$ is the collapse curve, $F(t)$ is the cumulative sum of sorted fluctuation of ensemble activity, $F_{min}$ and $F_{max}$ are the minimum and maximum values of $F(t)$, respectively, and $t$ is the time scale.

.9 Ensemble activity partition

Let $A$ denote the fraction of ensemble activity of a given functional category (either excitation or inhibition), let $B$ denote the discrete order of $A$ given by values $b_{1},b_{2}...,b_{n}$ (in non-decreasing order) and their probabilities $g(b_{j}):=\Pr(B=b_{j})$. The partition curve is defined as a continuous piece-wise linear function $(P)$ connecting the points $(C_{i},P_{i})$ such that:

$$C_{i} :=\sum_{j=1}^{i}g(b_{j})$$ (3)

$$D_{i} :=\sum_{j=1}^{i}g(b_{j})\,b_{j}$$

$$P_{i} :={\frac{D_{i}}{D_{n}}}$$

,with $B$ bounded between 0 and 1, the partition curve represents a size-based distribution (as a list of values) that indicates the percentage of ensemble activity of a certain size or in a certain size interval. The partition curve is convex and lies below the diagonal (45 degree) line which represents equivalence of normalized ensemble activity across time. The closer the convex partition curve is to the diagonal line, the lower the ensemble activity unequivalence. The further away the partition curve is from the 45 degree diagonal line, the greater the unequivalence in the distribution of ensemble activity. Color-coded curves in Fig.3.A show the partition curves for various scales (darker colors represent finer scales while brighter colors represent coarser scales). The black partition curve is calculated from a normal distribution for comparison.

.10 Multiscale descriptors of co-occurrence of Excitation and Inhibition

To analyze the properties of the joint probability of the fraction of Excitation and Inhibition across many temporal scales, we implemented a multidimensional textural feature analysis. From the multidimensional space of excitation/inhibition/scale (EIS), we calculated a “Multiscale Normalized Co-occurrence Matrix” (MNCM) for multiple distances and orientations following the volumetric adaptation [60, 61] of textural analysis [62]. The MNCM was created by calculating the frequency of occurrence of pairs of neighboring discrete points (pixels) in the multidimensional EIS space with specific normalized values. The spatial relationship between the pixels can be defined using an “offset” parameter, which relates to the spatial co-occurrence of E vs I in the EIS. Each entry in the MNCM indicates the number of times a pixel with value $i$ is adjacent to a pixel with value $j$. This spatial dependence matrix provides the distribution of similar values (of E and I) at given offsets, indicating the number of times that a given value of ensemble E co-occurred with a given value of ensemble I, considering the EIS space. We then used MNCM to measure features of distinct properties, such as contrast category (contrast, homogeneity, dissimilarity), descriptive measures (mean, variance, correlation), and orderliness features (entropy, angular second moment or energy).

Considering $P(i,j)$ as the joint distribution of $i$ (ensemble excitation dimension) and $j$ (ensemble inhibition dimension) in the original multidimensional matrix, $p(i,j)$ is the $(i,j)th$ entry in a normalized MNCM (i.e., spatial dependence matrix). The normalization of $P(i,j)$ is performed by dividing each entry by the sum of all its values, ensuring the sum of all $p(i,j)$ is 1. $N_{g}$ is the number of distinct normalized levels in the multidimensional space of excitation/inhibition/scale obtained after the operation of MNCM.

The features from the category of “descriptive statistics” were defined as:

$$Correlation=\cfrac{\sum_{i}\sum_{j}(i,j)p(i,j)-u_{x}u_{y}}{\sigma_{x}\sigma_{y}},\\$$ (4)

where $p(i,j)$ represents the normalized joint probability density of observing excitation level $i$ with inhibition level $j$. The terms $u_{x}$ and $u_{y}$ are the means of the marginal distributions $p_{x}(i)$ (probability of excitation level $i$) and $p_{y}(j)$ (probability of inhibition level $j)$, respectively, and $\sigma_{x}$ and $\sigma_{y}$ are their standard deviations. The correlation feature thus quantifies the degree of linear dependency between the excitation and inhibition levels across different scales.

Measures from the “contrast” group were:

$$Contrast=\sum_{i=1}^{N_{g}}\sum_{j=1}^{N_{g}}(i-j)^{2}p(i,j)\\$$ (5)

$$Dissimilarity=\sum_{i=1}^{N_{g}}\sum_{j=1}^{N_{g}}(i-j)p(i,j)\\$$ (6)

(similar to contrast, but instead, with linear weights.)

$$Homogeneity=\sum_{i=1}^{N_{g}-1}\sum_{i=1}^{N_{g}-1}\cfrac{1}{1+(i-j)^{2}}p(i,j)$$ (7)

and the features of “orderliness” were:

$$Energy=\sum_{i=1}^{N_{g}}\sum_{j=1}^{N_{g}}p(i,j)^{2}\equiv sqrt(\>Angular\>Second\>Moment)\\$$ (8)

$$Entropy=-\sum_{i=1}^{N_{g}}\sum_{j=1}^{N_{g}}p(i,j)log(p(i,j))$$ (9)

The “Descriptive Statistics Group”, “Contrast Group”, and “Orderliness Group” capture different aspects of the MNCM: ‘Descriptive Statistics Group is characterized by statistical metrics that describe the overall distribution of point intensities in a given scale. One key feature of this group is “Correlation”, which measures the joint probability occurrence of specified pairs. Contrast Group includes several measures of contrast or variation in intensity within an scale. “Homogeneity” measures the closeness of the distribution of elements in the MNCM to the MNCM diagonal. “Contrast” measures the local variations in the MNCM. “Dissimilarity” measures the variation of pairs in the MNCM. Orderliness Group includes measures of order or regularity within a given scale. “Energy” (also known as second angular moment) provides the sum of squared elements in the MNCM, and “Entropy” measures the complexity or disorder in a given scale. The application to volumetric (multiscale) data involves computing these features for each scale independently and then combine the resulting feature vectors into a single feature vector representing the entire multidimensional space.

.11 Wavelet Coherence and Ridge Detection

To relate the frequency structure of the local field potential (LFP) with the estimated synaptic current, we used wavelets to calculate coherence [69]. The wavelet cross-spectrum is defined as the product of the wavelet transform of signal $(x)$ and the complex conjugate of the wavelet transform of signal $(y)$. The squared absolute value of the smoothed cross-spectrum normalized by the smoothed individual power spectra defines wavelet coherence.

$$W_{n}^{XY}(s)=W_{n}^{X}(s)W_{n}^{Y^{*}}(s)\\$$ (10)

$$C_{n}^{2}(s)=\frac{\left|\left\langle W_{n}^{XY}(s)\right\rangle.s^{-1}\right|^{2}}{\left\langle\left|W_{n}^{XX}(s)\right|.s^{-1}\right\rangle\left\langle\left|W_{n}^{YY}(s)\right|.s^{-1}\right\rangle}$$ (11)

To extract ridges from the wavelet, we used the crazy-climber algorithm. This stochastic relaxation method is based on Markov chain Monte Carlo (MCMC) and can extract multiple ridges from analyzed signals [108, 109]. In its algorithmic procedure, the energetic distribution in the time-frequency representation of the LFP is used to extract multiple ridges across frequency ranges of LFP.

References