Scale-Change Symmetry in the Rules Governing Neural Systems.

Similar universal phenomena can emerge in different complex systems when those systems share a common symmetry in their governing laws. In physical systems operating near a critical phase transition, the governing physical laws obey a fractal symmetry; they are the same whether considered at fine or coarse scales. This scale-change symmetry is responsible for universal critical phenomena found across diverse systems. Experiments suggest that the cerebral cortex can also operate near a critical phase transition. Thus we hypothesize that the laws governing cortical dynamics may obey scale-change symmetry. Here we develop a practical approach to test this hypothesis. We confirm, using two different computational models, that neural dynamical laws exhibit scale-change symmetry near a dynamical phase transition. Moreover, we show that as a mouse awakens from anesthesia, scale-change symmetry emerges. Scale-change symmetry of the rules governing cortical dynamics may explain observations of similar critical phenomena across diverse neural systems.

Introduction

The ongoing collective population activity of neurons in cerebral cortex exhibits complex spatiotemporal fluctuations. This ongoing activity is responsible for the majority of the brain's energy consumption (Buzsáki et al., 2007), is closely related to past experiences (Han et al., 2008, Kenet et al., 2003, Luczak et al., 2009), contributes to memory consolidation (Gupta et al., 2010, Ji and Wilson, 2007), and modulates ongoing cortical processing (Arieli et al., 1996, Fox et al., 2007, Fox et al., 2006, Petersen et al., 2003). Thus understanding how ongoing cortical activity is organized is an important goal of systems neuroscience. Comparing ongoing cortical activity across diverse species, measured with different experimental methods over the past decade, a common phenomenon has been found with surprising consistency. The spatiotemporal sizes of ongoing fluctuations follow a specific statistical law; they are distributed according to a power law probability density function with the same exponent near −1.5. This phenomenon has been observed in functional magnetic resonance imaging (fMRI) in humans (Haimovici et al., 2013, Tagliazucchi et al., 2012); magnetoencephalography (MEG) in humans (Shriki et al., 2013); voltage imaging in mice (Fagerholm et al., 2016, Scott et al., 2014); multielectrode electrophysiology in monkeys (Petermann et al., 2009, Yu et al., 2017), cats (Hahn et al., 2017), rats (Gautam et al., 2015, Gireesh and Plenz, 2008), and even turtles (Clawson et al., 2017, Shew et al., 2015); and in vitro brain-in-a-dish systems (Beggs and Plenz, 2003, Shew et al., 2011). What can explain such shared phenomena observed across such diverse cortical systems, measured with differing spatial and temporal resolution?

One hypothesis is that the cerebral cortex can operate near the critical point of a phase transition (Beggs and Timme, 2012, Hesse and Gross, 2014, Plenz et al., 2014, Shew and Plenz, 2013). This hypothesis builds upon well-established physics of critical phase transitions; at criticality, multiple material properties are expected to be power law distributed according to universal scaling laws that are insensitive to many details of the physical system (Stanley, 1987, Stanley, 1971, Wilson, 1979, Wilson, 1975). From theory and computational modeling work, it is clear that neural systems can be tuned into a variety of different dynamical regimes or phases (e.g., asynchronous, oscillatory, bursting), often with distinct boundaries separating different regimes (Brunel, 2000, Gautam et al., 2015, Haldeman and Beggs, 2005, Poil et al., 2012, Wang et al., 2011). Similar to physical systems, near certain regime boundaries these neural models exhibit power-law-distributed collective dynamics with the same power law exponents, despite many detailed differences among models (Gautam et al., 2015, Haldeman and Beggs, 2005, Wang et al., 2011). What does the physics of critical phenomena tell us about where this universality comes from?

Fundamental insight into the origins of universal critical phenomena in physical systems (and a Nobel Prize) came from the realization that the basic laws governing the system obey a peculiar symmetry. The laws are the same across different scales; they have scale-change symmetry. Critical phenomena, including power-law-distributed observables, stem directly from this bizarre fractal symmetry of the governing physical laws (Wilson, 1979, Wilson, 1975). Motivated by this basic fact about critical phenomena in physical systems, here we hypothesize that the governing laws for cortical dynamics also conform to scale-change symmetry. Despite more than a decade of intense research on criticality in neural systems, this basic question has not been addressed. How does one go about testing this hypothesis?

Here, we develop an approach inspired by renormalization group theory, but with a focus on practical applicability to real neural data. Renormalization group theory is the mathematical approach used to understand scale-invariance at criticality in equilibrium physical systems (Stanley, 1999, Stanley and Wong, 1972, Wilson, 1979, Wilson, 1975), and in certain non-equilibrium dynamical systems (Loreto et al., 1995, Tauber, 2014, Vespignani et al., 1997, Vespignani et al., 1995). However, renormalization group ideas have not been developed in the context of neural systems. In brief, the idea begins with a coarse-graining procedure that transforms all the system variables and governing laws at one spatiotemporal scale to a new set of variables and laws at a coarser scale. If the system obeys the scale-change symmetry, then an appropriately chosen coarse-graining procedure will leave the governing laws unchanged and leave the system variables with identical statistics.

Here, we first applied our approach to two computational models, one simple and the other more biologically realistic. For both models, we found that dynamical rules were indeed maximally scale-invariant when operating near a phase transition. Next, we applied our approach to experimental data. We found that, in the awake state, the apparent rules governing dynamics of mouse cerebral cortex were more scale-invariant than in the anesthetized state. This suggests that the unconscious cortex deviates further from criticality than does the conscious cortex.

Results

Simple Model

We developed our approach using a network of binary, probabilistic, excitable nodes similar to that used in previous studies of non-equilibrium critical phenomena in neural systems (Haldeman and Beggs, 2005, Kinouchi and Copelli, 2006, Larremore et al., 2011, Shew et al., 2009), but on a two-dimensional L X L square lattice (L = 400) with nearest-neighbor connections (including self). We interpret each node in the network as the aggregate state of a group of neurons, analogous to the signals measured using experimental techniques with somewhat coarse spatial resolution (e.g., local field potential, Petermann et al., 2009; voltage imaging, Scott et al., 2014; or fMRI, Haimovici et al., 2013, Tagliazucchi et al., 2012). The state X(t+1) of a node at time t+1 is either active (1) or inactive (0) depending on the number of neighbors n(t) that were active at time t.and C defines the strength of coupling between all connected pairs of nodes and p defines the probability of activation without any active neighbors (p = 0.001 unless stated otherwise). We interpret p as an external source of input to the neural system. The states of all the network nodes are updated synchronously. One way to characterize the collective network dynamics is by measuring the fraction of active nodes at a time, or network activity, given by . By tuning the coupling strength, we tuned the model through a critical phase transition at C = C∗, separating a low-activity (small S) from a high-activity (large S) regime (Figure 1). Most of our results consider the range of C from 0.15 to 0.35.

Figure 1

Phase Transition in a Simple Neural Model

(A) Each panel shows the two-dimensional lattice of nodes at a single time step. Each pixel represents one node (yellow, active; blue, inactive). A subset of the full lattice is shown for clarity.

(B) As coupling strength C increases a sharp increase in time-averaged network activity occurs at a critical coupling strength near  = 0.23. is averaged over 104 time steps excluding a transient period of 102 time steps.

Our first goal was to examine how the dynamical rules governing the system change when considered at different scales of observation. For this, we consider six activation probabilities with n = 0, 1, 2, 3, 4, or 5 (n = 5 means all 4 neighbors plus self). Together, these six probabilities completely specify the rules of the system dynamics (actually they over-specify the rules, which are completely specified by C and p, but over-specification can be helpful with noisy data). Importantly, we can estimate these probabilities directly from data, whether simulated or experimentally measured. For instance, is estimated as the fraction of instances with n(t) = 1 that lead to activation at time t+1. Based on 104 time steps of data from the model (excluding the transient of duration 102 time steps), we can estimate all six probabilities with precision ≤ 2% for 0.15 < C < 0.35 (except , which was less accurate for C > 0.3; see Figure S1).

The scale-change symmetry at the heart of critical phenomena in physical systems suggests that the could be invariant across observational scales when C equals its critical value C*, but not if C deviates from C*. The next step toward testing this possibility was to identify an appropriate coarse-graining renormalization transformation to map one observational scale onto a coarser scale.

Coarse-Graining Scheme

Motivated by the limited spatiotemporal resolution of typical experimental data and real-space block renormalization schemes originally used to gain insight into the critical phenomena in physical systems (Kadanoff, 1993, Stanley and Wong, 1972, Wilson, 1979), we devised a spatiotemporal transformation scheme described as follows. Each r × r ×τ spatiotemporal block of nodes (r × r in space and of duration τ) at scale b transforms into a single node at coarser scale b + 1. Thus, an L × L × T lattice at observational scale b transforms into a coarse L/r × L/r × T/τ lattice at observational length scale b + 1, as illustrated for r = 4 and τ = 1 in Figures 2A and 2B. The state X of a coarse node depends on the average state of the corresponding τr2 nodes at the finer scale according to

Figure 2

Scale-Invariance of Dynamical Rules Peaks at Criticality

(A) Cartoon illustration of coarse-graining scheme. Each block of nodes at fine scale b is transformed probabilistically to one node at the coarse scale b + 1.

(B) Examples of activity snapshots before and after coarse graining.

(C) Upon coarse graining, the dynamical rules change the least ( is minimal) at criticality. Inset shows the coarse-graining transformation function with (k; x0) = (76; 0.22). Block size was r = 8.

(D) Shown are optimal coarse-graining functions for three example C values and six block sizes (legend in E specifies different values of r and τ).

(E) Using the optimal coarse-graining function for each C resulted in the strongest scale-invariance of dynamical rules, i.e., lowest around . This result held for multiple choices of block size and duration (see legend).

(F) The valley in as a function of coupling strength C became broader as p was increased. For A–E p = 0.001.

By tuning the transformation parameters k and x0, we can explore a family of logistic function coarse-graining procedures including majority rules and other previously studied functions. Similar coarse-graining schemes have been used in previous work on deep learning (Hinton and Salakhutdinov, 2006).

To quantify the change in the dynamical rules due to coarse graining, we define a parameter ζwhere is the activation probability estimate at scale b and is the activation probability estimate at the coarser scale b + 1. Note that root-mean-square differences or any other monotonic function of the differences of and would not change our following conclusions. We ran the model for a range of coupling strengths 0.15 < C < 0.35, simulating 104 time steps for each C. In line with expectations from critical phenomena in physical systems, we found that was minimized for C near C∗ = 0.23 (Figure 2C). This result was obtained using transformation parameters k = 76, x0 = 0.22, and block size r×r×τ = 8 × 8 × 1 (Figure 2C, inset). This means that the dynamical rules governing the network activity are least changed at C∗.

However, some caution is appropriate; if a different transformation function f was used, perhaps the minimum would occur at a different C. Therefore, to draw a more definitive conclusion, we next systematically searched the two-dimensional space of all possible functions, seeking the function that minimizes for each C independently. We found that, when using such optimal f functions, the dynamical rules remained most scale-invariant (i.e., lowest ) for C near C∗ (Figures 2D and 2E). Thus our first conclusion is that scale-invariance of the governing laws at criticality, one of the most fundamental concepts of critical phenomena in physical systems, can also manifest in neural systems.

We also verified our results for coarse-graining blocks with different spatial sizes r = 4, 8, 16 (Figures 2D and 2E), different durations τ = 1, 2, 4, 8 (Figures 2D and 2E), and different levels of noise p = 10−2, 10−3, 10−4 (Figure 2F). Optimal x0 increased with C and r; optimal k increased with r and was lowest near C = C* (Figure 2D). Changes in r, τ, and p did not qualitatively change our conclusion that coarse graining causes the least change in the dynamical rules near criticality. However, increasing p, decreasing , and increasing τ all had the effect of broadening the minimum in and shifting the minimum toward slightly smaller values of C. For increasing p, this shift in the minimal C may be due to the phase transition becoming “smeared out,” less sharp with increased noise (Williams-García et al., 2014), but further investigation is needed to test this possibility. We chose the range 0.15 < C < 0.35 because at lower values (C < 0.15) the rates of activity became so small that our estimates of the became poor due to subsampling, which made an unreliable measure (see Figure S2). A skeptical reader may note that is not zero at criticality for our model; the scale-invariance we observe is imperfect. Why might this be? In the calculation of , we found that the largest contribution to came from the term in the sum, which represents spontaneous activation of a node. This is because the coarse-graining procedure tends to cause periods with very low S at scale b to become S = 0 at scale b+1, thus creating excessive apparent spontaneous activation at scale b+1. When is calculated by excluding it still minimizes for C near C∗ but with a much lower value of , indicating higher degree of scale-invariance for the dynamical rules governing interactions, i.e., those with n > 0 (see Figure S3).

Realistic Model

The abstract binary model we have discussed so far is simple to interpret and directly comparable to previous work on critical phenomena in neural systems using similar models. Moreover, it is relatively similar to models of spreading dynamics in previous renormalization group studies (e.g., sand piles, Vespignani et al., 1995, and forest fires, Loreto et al., 1995). In this sense, it is not terribly surprising that our simple model conformed to scale-change symmetry of governing rules, which is expected from theory of critical phenomena. However, our simple binary model is not very biologically realistic. For example, the dynamics of real cerebral cortex never exhibits a sustained high-firing regime as we see in the simple model with C > C∗. In reality, if the excitability of real cortex is enhanced (or other parameters analogous to C are enhanced) the dynamics tends to exhibit large, repetitive bursts of activity that are suppressed by depressive adaptive mechanisms before reaching a sustained high-firing regime (Gautam et al., 2015, Shew et al., 2009). Indeed, activity-dependent adaptive effects can act to make the critical regime more robust (Levina et al., 2007, Shew et al., 2015), which could make the possibility of scale-invariant dynamical rules even more plausible in real brains. Next, we set out to test our ideas in a more realistic model of a neural network that generates more realistic network activity. Similar to other recent model studies (di Santo et al., 2018, Gautam et al., 2015), the phase transition we examine in our more realistic model separates an asynchronous phase from an oscillatory ordered phase, which better matches experimental observations.

Building on our simple model, we kept binary, probabilistic, integrate-and-fire neurons on a two-dimensional square lattice (L × L = 160 × 160), but now we introduced inhibitory neurons (20% of all neurons), spike-frequency adaptation, refractoriness, and different distance-dependent connectivity for excitatory and inhibitory neurons. At each time t, the binary state X(t) = 1 with probability where

Input from other neurons is and the activity-dependent adaptation is modeled by where τ = 80. If this sum is zero, we set h to 1. At each time step a neuron can also be activated due to external sources with probability p = 0.001 (treated independently of p). After a neuron fires, its state is set to 0 for a refractory period of 1 time step. The default synaptic weight matrix W is constructed with long-range inhibition relative to shorter-range excitation. Weights of inhibitory synapses were and those of excitatory synapses were , where d is the distance from the presynaptic neuron to the postsynaptic neuron. We consider C = 3 and C = 2 (in units of lattice spacing) and fixed excitatory input strength E = 1. By tuning the inhibitory input strength I, we tuned the model through a phase transition near I = 0.65 separating a synchronous oscillatory firing state (high mean pairwise correlation) from asynchronous low-firing (low mean pairwise correlation) regime (Figures 3A–3C).

Figure 3

Scale-Invariance of Dynamical Rules Peaks at Phase Transition in a More Biologically Plausible Model

(A) Each panel shows the two-dimensional lattice of neurons at a single time step. Each pixel represents one neuron (yellow, active; blue, inactive). The spatiotemporal dynamics was limited to small scales for strong inhibition (I = 2.0, bottom row), exhibited massive propagating waves and oscillations for weak inhibition (I = 0.01, top row), and had more complexity near the transition between these extremes (I = 0.65, middle row).

(B) Time series of network activity reveals the prominent oscillatory activity of the weak inhibition regime (red).

(C) As inhibition is increased, the boundary of the oscillatory regime near I = 0.65 (dashed line) is revealed by the drop in mean pairwise correlations.

(D) Scale-invariance of dynamical rules peaked ( is minimal) near the onset of the oscillatory regime. This held for blocks with different spatial sizes and durations (see legend).

Like the simple model, we found that is minimized near the phase transition (Figure 3D). We emphasize that this more realistic model is outside the bounds of well-understood critical phenomena in physics. Although it shares some features with simple models of spreading dynamics (e.g., excitable nodes with refractory periods), the presence of inhibition and nonlocal connectivity makes it quite different from these previously studied models. Thus there is no guarantee that the phase transition we consider here has anything to do with criticality. In this sense, it is substantially more surprising that we observe minimal scale change of dynamical rules near this phase transition. For the results in Figure 3D, we only considered with n = 1, 2, 3, 4 for calculation, because tended to be most prone to error as discussed above (Figure S3). We did not include because the refractory period precluded the occurrence of at scale b.

So far, we have applied our approach to spike data from our two models. However, many experimentally measured signals are continuous and more closely related to membrane potential, particularly at the larger spatial scales (e.g., local field potential (LFP), fMRI, electroencephalography, voltage imaging, wide-field calcium imaging). How might we apply our approach to such continuous signals? To address this question, here, we first used our more realistic model. We generated a continuous signal from the model, using the “membrane potential” from the model neurons. We defined the membrane potential of the i neuron to be χ(t)h(t)−1, as defined above. Next, we binarized the continuous membrane potential signal. Following a previously established approach (Scott et al., 2014, Tagliazucchi et al., 2012), we defined time points of activation when the voltage imaging time course crossed above a threshold of 0.5 SD beyond the mean from below (results consistent for different thresholds 0.25 SD and 1 SD) (Figure 4). We found that, like the spike data presented above, this binarized continuous signal also supported our main claim: minimal near the phase transition.

Figure 4

Applying Our Approach to Continuous Synaptic Input

(A) Mean pairwise correlations of binarized membrane potential for the realistic model.

(B) Change in dynamical rules governing the binarized membrane potential as a function of inhibition strength I for r = 8 (left), r = 16 (right), and different binarization thresholds (color). For all the cases shown τ = 1 and network size, L × L = 160 × 160.

Mouse Cerebral Cortex

One advantage of our approach is that it can readily be applied to experimental data, provided the data have sufficient spatial resolution and coverage. Here, we demonstrate this for measurements of cerebral cortex in a mouse awakening from anesthesia. Previous work (Scott et al., 2014) with these data suggested that as the mouse awakens the cortical dynamics transitioned from a supercritical regime, similar to the oscillatory regime in our realistic model, toward criticality. The brain activity was measured over the dorsal surface of nearly an entire hemisphere of mouse cortex with high spatial resolution using genetically encoded voltage-sensitive fluorescence imaging. Each dataset was acquired during 1 min with 50 Hz sample rate and 33 × 33μm2 per pixel spatial resolution. Multiple such 1-min recordings were performed over a period of 200 min as the mouse awoke. Each pixel represents the aggregate activity of many neurons within cortical layers 2 and 3. This voltage imaging signal is a continuous signal. To apply our approach, we followed the same binarization approach described above for the membrane potential analysis of the realistic model (Figure 4). We converted each pixel voltage time series into binary form (Figure 5A). For each pixel and each 1-min recording, time points of activation were defined at times when the voltage imaging time course crossed above a threshold of 0.5 SD beyond the mean from below (results consistent for different thresholds 0.25 SD and 1 SD).

Figure 5

Increase in Scale-Invariance of Cortical Dynamical Rules as Mouse Awakens

(A) Genetically encoded voltage-sensitive fluorescence imaging was done to measure the spatiotemporal dynamics across one hemisphere of mouse cortex as it awoke from anesthesia. Each panel shows a snapshot of binarized activity (yellow, active; blue, inactive). The signal of each pixel arises from many neurons within a 33 × 33 μm2 area.

(B) Time series of binary network activity datasets. Under anesthesia (red), the dynamics exhibited relatively large-scale bursts, whereas the awake dynamics (blue) tended to be more diverse.

(C) Mean pairwise correlation decreases as the mouse awakens.

(D) Scale-invariance of dynamical rules increases (ζ decreases) as the mouse awakens. Results were qualitatively consistent for three different binarization thresholds (yellow, red, and blue) and two different coarse-graining block sizes (r = 8 and 16).

Under the effect of anesthesia, we found that the mouse exhibited synchronous burst firing, similar to the oscillatory regime of our realistic model. As the mouse woke up, more asynchronous firing was observed (Figure 5B), similar to the activity near the onset of the oscillatory regime of our realistic model. This change manifested as a decrease in pairwise correlations as the mouse awoke from anesthesia (Figure 5C). Next, we assumed nearest-neighbor interactions among pixels and proceeded to estimate the activation probabilities . Then, we applied our course-graining procedure and estimated . As with our analysis of the realistic model, only four activation probabilities, n = 1, 2, 3 4 were used to estimate . We found that decreased as the mouse awoke from anesthesia. This finding demonstrates that the rules governing cerebral cortex dynamics approach scale-change symmetry during the transition from anesthetized to awake.

We note that, for higher values of τ, there was increased noise in the trend relating to time since drug delivery. As τ increases, we are left with fewer time points with which we estimate the rules and calculate . This decrease in samples may be responsible for poorer estimates of the rules and a noisier versus time trend. Future studies with longer-duration experimental recordings could better test this possibility.

Relating Rules to Dynamics

The approach we have developed here offers a new way to learn about whether a system is operating near criticality by examining how the rules governing dynamics change with scale. In contrast, the traditional way to assess whether a system is near criticality is by examining the dynamics, rather than the rules. For instance, the experimental data we analyze here (Figure 5) has previously been shown to exhibit dynamics that are consistent with an approach to criticality as the mouse wakes up (Fagerholm et al., 2016, Scott et al., 2014). This claim was largely based on examining cascades of propagating neural activity, often called neuronal avalanches. At criticality, theory predicts that different sizes of avalanches occur with probability that is related to size according to a power law with exponent near −1.5. In contrast, in the ordered phase (like the oscillatory phase in our realistic model) avalanche sizes are expected to deviate from a power law distribution, with very large cascades becoming more prominent. When Scott et al. analyzed our experimental data, they found that the anesthetized state was consistent with an ordered phase, whereas the awake state was more consistent with criticality. This finding is consistent with our conclusion that scale-change symmetry increases as the mouse wakes up. However, we emphasize that a power law distribution of avalanche sizes is certainly not equivalent to our finding of scale-change symmetry of rules. Indeed, power law distributions of observed dynamics can arise due to mechanisms that are totally unrelated to criticality (Beggs and Timme, 2012, Mitzenmacher, 2003, Reed and Hughes, 2002, Sornette, 1998, Stumpf and Porter, 2012). If criticality is responsible for a power law distribution of an observable, then we would expect to also find scale-change symmetry of rules. To clarify how scale-change symmetry of rules is related to power law avalanche size distributions, we studied our models further. We examined avalanche size distributions for the simple model as C was tuned through C*, and for the more realistic model as inhibition (I) was tuned from strong (I = 2) to weak (I = 0.01). We used the previously developed measure κ to quantify deviation from a power law (Scott et al., 2014, Shew et al., 2009). In brief, κ = 1 for a perfect match to the reference power law with exponent −1.5, κ> 1 when large avalanches become prominent (as expected in the ordered phase), and κ< 1 when large avalanches become rare (as expected in the asynchronous phase). Further description of κ is given in the Methods section.

As expected for the simple model, we found that the distribution of avalanche sizes was close to a power law (κ near 1) when C was tuned near C* (Figures 6A and 6B). Thus, in this case, was also smallest near κ = 1 (Figure 6C). Comparing κ and , we found that κ is more sensitive to changes in C for C > C*, whereas is more sensitive for C < C*. For the more complex model, the oscillatory regime (low inhibition) showed a clear bump in the tail of the avalanche distributions similar to previous studies of dynamics in the ordered (supercritical) phase (Figure 6D). As inhibition was tuned from weak to strong, the avalanche size distributions did approach approximate, but rather imperfect power laws. The changes in avalanche distributions we observed for increasing inhibition is similar to what has been shown previously for our experimental data as the mouse wakes up (Scott et al., 2014). Despite the imperfect power laws found in the realistic model, we computed κ and found a clear tendency for minimal near κ = 1 (Figures 6E and 6F). These results demonstrate a close relationship between avalanche statistics and the scale-change symmetry of rules governing a neural system.

Figure 6

Scale-Invariance of Rules Versus Avalanche Size Distributions

(A) Shown are avalanche size distributions obtained from the simple model with different values of coupling, C. The probability for large avalanches is prominent for strong coupling and dramatically lower for weak coupling. Distributions are shifted vertically for visual comparison. Black dashed line indicates a power law with exponent −1.5.

(B) The parameter κ measures deviation between a measured avalanche size distribution and −1.5 power law. Near C = C *, we found minimal deviation from power law (κ = 1).

(C) We found minimal change in rules near κ = 1.

(D) For the realistic model, avalanche size distributions exhibited high probability for large avalanches when inhibition was small (blue) and approximate power law distributions for stronger inhibition.

(E) Near the onset of the oscillatory phase, we found the smallest deviation from power law (κ near 1).

(F) Change in rules was minimal near κ = 1.

Discussion

Here we have shown that in neural network models, the rules that govern the system behavior obtain a degree of scale-change symmetry that is most pronounced when the system operates near criticality. In the spirit of renormalization group theory in physical systems, our approach offers a potential explanation of why diverse neural systems can exhibit very similar critical phenomena.

Another useful outcome of our work is that the tools we develop provide a way to assess whether a change in system dynamics takes the system closer to or further from criticality. Our approach complements traditional approaches for seeking evidence for criticality based on power laws and scaling relations. Such power law distributions constitute necessary, but insufficient, evidence for the criticality hypothesis (Beggs and Timme, 2012, Klaus et al., 2011, Stumpf and Porter, 2012). Our results suggest that a change toward criticality should be accompanied by an increase in scale-invariance of effective dynamical rules. To be more specific, considering our model work (both models, Figures 2, 3, and 4), it is clear that, for a wide range of parameters (C or I) around criticality, a decrease in ς does not happen unless the parameters have been tuned toward criticality—tuning C toward C* (Figure 2) or tuning I toward the value that corresponds to the onset of oscillations (Figures 3 and 4). However, for a given observed decrease in , we cannot make precise conclusions about how much closer to criticality the system has shifted. This is not possible considering the quantitative differences in the shape of the versus C and versus I curves and how they depend on the details of analysis parameters (e.g., block size). The stronger conclusions that can be made are about (1) the direction of change and (2) about the relative size of changes within one system. For example, our methods could be used to ascertain which kinds of pharmacological manipulations result in a bigger shift away or toward criticality in a single system with consistent measurement tools. This is an important value of our approach—measuring relative changes in proximity to criticality.

Limitations of Study

One limitation of our approach is that it currently requires binary data (Figures 1, 2, and 3) or binarization of continuous data (Figures 4 and 5). Although spike data are well suited to a binary approach, it would be useful for future work to generalize the approach to continuous data, which is more common in experiments. However, we emphasize that our results based on binarization of continuous data from our realistic model (Figure 4) were qualitatively consistent with our results based on spikes. This suggests that the binarization of continuous data provides a useful strategy for studying scale-change symmetry of dynamical rules.

Another limitation of our approach is that we have based our assessment of scale-change symmetry on one step of renormalization. Renormalization group theory, in contrast, emphasizes that many steps of renormalization may be required before true scale-invariance is found. Each step of renormalization is thought to remove irrelevant details that could corrupt scale-invariance. This fact may contribute to why our scale change measure is not very close to zero even when our simple model is at criticality. However, the fact that we do find a prominent minimum in near criticality indicates that we can detect the scale-change symmetry in spite of only performing one step of renormalization. This is important, because, when working with real data from finite systems, it is not plausible to do many steps of coarse graining especially with a large block size.

We also note that there are some trivial types of scale-change symmetry that we are not interested in and have not highlighted in our results above. For instance, for a completely saturated system—all neurons firing at every time step—the dynamics and rules would be identical upon coarse graining. The coarse level would still be completely saturated. Likewise, a completely silent network would also exhibit such trivial scale-invariance. These trivial types of invariance are not the subject of our work, but should be noted for completeness.

Our initial application of our approach to experimental data provided interesting results. We showed that scale-invariance of the rules governing cortical network dynamics increased during the transition from unconsciousness to consciousness. This finding is in line with other recent work that suggests that the anesthetized cortex deviates from critical dynamics and approaches criticality as it wakes up (Bellay et al., 2015, Fagerholm et al., 2016, Scott et al., 2014). A recent coarse-graining study of experimental data obtained in hippocampus also revealed interesting scaling, in line with the scale-change symmetry we report here (Meshulam et al., 2018). Our finding raises interesting questions about the functional consequences of scale-invariant rules. What does it mean that local interactions, say among cortical columns, are governed by the same rules as larger-scale interactions, say among cortical regions? We expect that our approach will be useful for future studies of these questions and how the dynamical rules governing brain dynamics differ across scales and brain states.

Methods

All methods can be found in the accompanying Transparent Methods supplemental file.