Neural criticality from effective latent variables

Abstract
eLife assessment
Introduction
Results
Discussion
Materials and methods
Data availability
References
Article and author information
Metrics

Abstract

Observations of power laws in neural activity data have raised the intriguing notion that brains may operate in a critical state. One example of this critical state is ‘avalanche criticality’, which has been observed in various systems, including cultured neurons, zebrafish, rodent cortex, and human EEG. More recently, power laws were also observed in neural populations in the mouse under an activity coarse-graining procedure, and they were explained as a consequence of the neural activity being coupled to multiple latent dynamical variables. An intriguing possibility is that avalanche criticality emerges due to a similar mechanism. Here, we determine the conditions under which latent dynamical variables give rise to avalanche criticality. We find that populations coupled to multiple latent variables produce critical behavior across a broader parameter range than those coupled to a single, quasi-static latent variable, but in both cases, avalanche criticality is observed without fine-tuning of model parameters. We identify two regimes of avalanches, both critical but differing in the amount of information carried about the latent variable. Our results suggest that avalanche criticality arises in neural systems in which activity is effectively modeled as a population driven by a few dynamical variables and these variables can be inferred from the population activity.

eLife assessment

This paper provides a simple example of a neural-like system that displays criticality, but not for any deep reason; it's just because a population of neurons are driven (independently!) by a slowly varying latent variable, something that is common in the brain. Moreover, criticality does not imply optimal information transmission (one of its proposed functions). The work is likely to have an important impact on the study of criticality in neural systems and is convincingly supported by the experiments presented.

https://doi.org/10.7554/eLife.89337.3.sa0

Introduction

The neural criticality hypothesis – the idea that neural systems operate close to a phase transition, perhaps for optimal information processing – is both ambitious and banal. Measurements from biological systems are limited in the range of spatial and temporal scales that can be sampled, not only because of the limitations of recording techniques but also due to the fundamentally non-stationary behavior of most, if not all, biological systems. These limitations make proving that an observation indicates critical behavior difficult. At the same time, the idea that brain networks are critical echoes the anthropic principle: tuned another way, a network becomes quiescent or epileptic and in either state, seems unlikely to support perception, thought, or flexible behavior, yet these observations do not explain how such fine-tuning could be achieved. Further muddying the water, researchers have reported multiple kinds of criticality in neural networks, including through analysis of avalanches (Beggs and Plenz, 2003; Plenz et al., 2021; O’Byrne and Jerbi, 2022; Girardi-Schappo, 2021) and of coarse-grained activity (Meshulam et al., 2019), as well as of correlations (Dahmen et al., 2019). How these flavors of critical behavior relate to each other or any functional network mechanism is unknown.

The phenomenon that we will refer to as ‘avalanche criticality’ appears remarkably widespread. It was first observed in cultured neurons (Beggs and Plenz, 2003) and later studied in zebrafish (Ponce-Alvarez et al., 2018), turtles (Shew et al., 2015), rodents (Ma et al., 2019), monkeys (Petermann et al., 2009), and even humans (Poil et al., 2008). The standard analysis, described later, requires extracting power-law exponents from fits to the distributions of avalanche size and of duration and assessing the relationship between exponents. There is debate over whether these observations reflect true power laws, but within the resolution achievable from experiments, neural avalanches exhibit power laws with exponent relationships predicted from theory developed in physical systems (Perkovic et al., 1995).

Avalanche criticality is not the only form of criticality observed in neural systems. Zipf’s law, in which the frequency of a network state is inversely proportional to its rank, appears in systems as diverse as fly motion estimation and the salamander retina (Mora and Bialek, 2011; Schwab et al., 2014; Aitchison et al., 2016). More recently, Meshulam et al., 2019 reported various statistics of population activity in the mouse hippocampus, including the eigenvalue spectrum of the covariance matrix and the activity variance. These were found to scale as populations were ‘coarse-grained’ through a procedure in which neural activities were iteratively combined based on similarity. Similar observations have been reported in spontaneous activity recorded across a wide range of brain areas in the mouse (Morales et al., 2023). Simple neural network models of such data explain neither Zipf’s law nor coarse-grained criticality (Meshulam et al., 2019).

Even though these three forms of criticality are observed through different analyses, they may originate from similar mechanisms. Numerous studies have reported relatively low-dimensional structure in the activity of large populations of neurons (Mazor and Laurent, 2005; Ahrens et al., 2012; Mante et al., 2013; Pandarinath et al., 2018; Stringer et al., 2019; Nieh et al., 2021), which can be modeled by a population of neurons that are broadly and heterogeneously coupled to multiple latent (i.e. unobserved) dynamical variables. Using such a model, we previously reproduced scaling under coarse-graining analysis within experimental uncertainty (Morrell et al., 2021). Zipf’s law has been explained by a similar mechanism (Schwab et al., 2014; Aitchison et al., 2016; Humplik and Tkačik, 2017). A single quasi-static latent variable has been shown to produce avalanche power laws, but not the relationships expected between the critical exponents (Priesemann and Shriki, 2018), while a model including a global modulation of activity can generate avalanche criticality (Mariani et al., 2021), but has not demonstrated coarse-grained criticality (Morrell et al., 2021). It is not known under what conditions the more general latent dynamical variable model generates avalanche criticality.

Here, we examine avalanche criticality in the latent dynamical variable model of neural population activity. We find that avalanche criticality is observed over a wide range of parameters, some of which may be optimal for information representation. These results demonstrate how criticality in neural recordings can arise from latent dynamics in neural activity, without need for fine-tuning of network parameters.

Results

Critical exponents values and crackling noise

We begin by defining the metrics used to quantify avalanche statistics and briefly summarize experimental observations, which have been reviewed in detail elsewhere (Plenz et al., 2021; O’Byrne and Jerbi, 2022; Girardi-Schappo, 2021). Activity is recorded across a set of neurons and binned in time. Avalanches are then defined as contiguous time bins, in which at least one neuron in the population is active. The duration of an avalanche is the number of contiguous time bins and the size is the summed activity during the avalanche. The distributions of avalanche size and duration are fit to power laws ( $P (S) \sim S^{- τ}$ for size $S$ , and $P (D) \sim D^{- α}$ for duration $D$ ) using standard methods (Clauset et al., 2009).

Power laws can be indicative of criticality, but they can also result from non-critical mechanisms (Touboul and Destexhe, 2017; Priesemann and Shriki, 2018). A more stringent test of criticality is the ‘crackling’ relationship (Perkovic et al., 1995; Touboul and Destexhe, 2017), which involves fitting a third power-law relationship, $\bar{S} (D) \sim D^{γ_{fit}}$ , and comparing $γ_{fit}$ to the predicted exponent $γ_{pred}$ , derived from the size and duration exponents, $τ$ and $α$ :

γ_{f i t} \overset{?}{=} γ_{p r e d} \equiv \frac{α - 1}{τ - 1} .

Previous work demonstrating approximate power laws in size and duration distributions through the mechanism of a slowly changing latent variable did not generate crackling (Touboul and Destexhe, 2017; Priesemann and Shriki, 2018).

Measuring power-laws in empirical data is challenging: it generally requires setting a lower cut-off in the size and duration, and the power-law behavior only has limited range due to the finite size and duration of the recording itself. Nonetheless, there is some consensus (Shew et al., 2015; Fontenele et al., 2019; Ma et al., 2019) that even if $τ$ and $α$ vary over a wide range (1.5 to about 3) across recordings, the values of $γ_{fit}$ and $γ_{pred}$ stay in a relatively narrow range, from about 1.1 to 1.3.

Avalanche scaling in a latent dynamical variable model

We study a model of a population of neurons that are not coupled to each other directly but are driven by a small number of latent dynamical variables – that is, slowly changing inputs that are not themselves measured (Figure 1A). We are agnostic as to the origin of these inputs: they may be externally driven from other brain areas, or they may arise from large fluctuations in local recurrent dynamics. The model was chosen for its simplicity, and because we have previously shown that this model with at least about five latent variables can produce power laws under the coarse-graining analysis (Morrell et al., 2021). In this paper, we examine avalanche criticality in the same model.

Figure 1

Download asset Open asset

Latent dynamical variable model produces avalanche criticality.

Simulated network is $N = 1024$ neurons. Other parameters in Table 1. (A) Model structure. Latent dynamical variables $h_{μ} (t)$ are broadly coupled to neurons $s_{i} (t)$ in the recorded population. (B) Raster plot of a sample of activity binned at 3 ms resolution across 128 neurons with five latent variables, each with correlation timescale $τ_{F} = 15 s$ . (C) Projection of activity into a simulated field of view for illustration. (**D-F**) Avalanche analysis in a network (parameters $N_{F} = 5$ , $τ_{F} = 10^{4}$ , $η = 4$ and $ϵ = 12$ ), showing size distribution (D), duration distribution (E), and size with duration scaling (F). Lower cutoffs used in fitting are shown with vertical lines and their values are indicated in the figures. There are $N_{o b s} = 42725$ avalanches of size $S \geq S_{m i n}$ in this simulated dataset. Estimated values of the critical exponents are shown in the titles of the panels.

Specifically, we model the neurons as binary units ( $s_{i}$ ) that are randomly ( $J_{i μ} \sim N (0, 1)$ ) coupled to dynamical variables $h_{μ} (t)$ . The probability of any pattern ${s_{i}}$ , given the current state of the latent variables, is

P (s_{i} | {h_{μ} (t)}) = \frac{1}{Z ({h_{μ} (t)})} \exp (- η \sum_{μ = 1}^{N_{F}} s_{i} J_{i μ} h_{μ} (t) - ϵ s_{i}),

where the parameter $η$ controls the scaling of the variables and $ϵ$ controls the overall activity level. We modeled each latent variable as an Ornstein-Uhlenbeck process with the time scale $τ_{F}$ (see Materials and methods). Thus our model has four parameters: $η$ (input scaling), $ϵ$ (activity threshold), $τ_{F}$ (dynamical timescale), and $N_{F}$ (number of latent variables).

Distributions of avalanche size and avalanche duration within this model followed approximate power laws (Figure 1C; see Materials and methods). In the example shown ( $N_{F} = 5$ , $τ_{F} = 10^{4}$ , $η = 4$ and $ϵ = 12$ ), we found exponents $τ = 1.89 \pm 0.02$ (size) and $α = 2.11 \pm 0.02$ (duration). Further, the average size of avalanches with fixed duration scaled as $S \sim D^{γ}$ , with the fitted $γ_{fit} = 1.24 \pm 0.02$ , in agreement with the predicted value $γ_{pred} = 1.24 \pm 0.02$ . Thus, our model could generate avalanche scaling, at least for some parameter choices. In the following sections, we examine how avalanche scaling depends on model parameters ( $N_{F}$ , $τ_{F}$ , $η$ and $ϵ$ ; see Table 2). We first focus on two sets of simulations: one set with $N_{F} = 1$ latent variable, which does not generate scaling under coarse-graining (Morrell et al., 2021), and one set with $N_{F} = 5$ latent variables, which can generate such scaling for some values of parameters $τ_{F}$ , $η$ , and $ϵ$ (Morrell et al., 2021; Table 1).

Table 1

Simulation parameters for Figure 1.

Parameter	Description	Value
$ϵ$	bias towards silence	$ϵ = 12$
$η$	variance multiplier	$η = 4.0$
$N_{F}$	number of latent fields	$N_{F} = 5$
$τ_{F}$	latent field time constant	$τ = 10^{4}$
$N$	number of cells	$N = 1024$

Avalanche scaling depends on the number of latent variables

We analyzed avalanches from one- and five-variable simulations, each with fixed latent dynamical timescale ( $τ_{F} = 5 \times 10^{3}$ time steps; see Table 2 for parameters). In the following sections, time is measured in simulation time steps, see Materials and methods for converting time steps to seconds. We used established methods for measuring empirical power laws (Clauset et al., 2009). The minimum cutoffs for size ( $S_{min}$ ) and duration ( $D_{min}$ ) are indicated by vertical lines in Figure 2. For the population coupled to a single latent variable, the avalanche size distribution was not well fit by a power law (Figure 2A). With a sufficiently high minimum cut-off ( $D_{min}$ ), the duration distribution was approximately power-law (Figure 2B).

Table 2

Simulation parameters for Figure 2.

Parameter	Description	Value
$ϵ$	bias towards silence	$ϵ = 8$ (for $N_{F} = 1$ ) or $ϵ = 12$ (for $N_{F} = 5$ )
$η$	variance multiplier	$η = 4.0$
$N_{F}$	number of latent fields	$N_{F} = 1 or 5$
$τ_{F}$	latent field time constant	$τ_{F} = 10^{3}, . . {.10}^{5}$
$N$	number of cells	$N = 1024$

Figure 2 with 4 supplements see all

Download asset Open asset

Multiple latent variables generate avalanche scaling at shorter timescales than a single latent variable.

Simulated network is $N = 1024$ neurons. Other parameters used for simulations for this figure are found in Table 2. (**A-C**) Scaling analysis for one variable models where the dynamic timescale is equal to 5×10³ time steps. (A) Distribution of avalanche sizes. MLE value of exponent for best-fit power law is $τ = 1.98$ (0.02 SE), with lower cutoff indicated by the vertical line. (B) Distribution of avalanche duration. MLE value of $α$ is 1.81 (0.02 SE). (C) Average size plotted against avalanche duration (blue points), with power-law fit (black line) and predicted relationship (yellow line) from MLE values for exponents in A and B. Gray bar on the horizontal axis indicates range, over which a power law with $γ = 1.72$ fits the data (see Materials and methods). (**D-F**) Analysis of avalanches from a simulation of a population coupled to five independent latent variables where the dynamic timescale is equal to 5×10³ time steps. (G) Exponents $τ$ for avalanche size distributions across timescales for one-variable (blue) and five-variable (red) simulations. Each circle is a simulation with independently drawn coupling parameters. Simulations had to show scaling over at least two decades to be included in panels (**G–J**). (H) Exponents $α$ for avalanche duration distributions for simulations in G. (I) Fitted values of $γ$ for simulations in G. (J) Difference between fitted and predicted $γ$ values. Five-variable simulations produce crackling over a wider range of timescales than single-variable simulations.

We next assessed whether the simulation produced crackling. If so, the value $γ_{fit}$ obtained by fitting $\bar{S} (D) \sim D^{γ_{fit}}$ would be similar to $γ_{pred} = \frac{α - 1}{τ - 1}$ . In many cases, such as the one-variable example shown in Figure 2C, the full range of avalanche durations were not fit by a single power law. Therefore, we determined the largest range, over which a power law was a good fit to the simulated observations. In this case, slightly over two decades of apparent scaling were observed starting from avalanches with minimum duration slightly less than 100 time steps (Figure 2C), with a best-fit value of $γ_{f i t} \in [1.69, 1.74]$ . As we did not find a power-law in the size distribution, calculating $γ_{pred}$ is meaningless. If we do it anyway, we obtain $γ_{p r e d} = 0.83 \pm 0.03$ (yellow line in Figure 2C), which clearly deviates from the fitted value of $γ$ . Thus, for the single latent dynamical variable model ( $τ_{F} = 5000$ ), power-law fits are poor, and there is no crackling.

The five-variable model produces a different picture. We now find avalanches, for which size and duration distributions are much better fit by power-law models starting from very low minimum cutoffs (Figure 2D–E, Figure 2—figure supplement 2). Average size scaled with duration, again over more than two decades, with $γ_{f i t} = 1.27 \pm 0.03$ , which was in close agreement with $γ_{p r e d} = 1.25 \pm 0.02$ (Figure 2F). Holding other parameters constant, we thus found that scaling relationships and crackling arise in the multi-variable model but not the single-variable model.

Avalanche scaling depends on the time scale of latent variables

Based on simulations in the previous section, we surmised that the five-variable simulation generated scaling more readily due to creating an ‘effective’ latent variable that had slower dynamics than any individual latent variable. We reasoned that at any moment in time, the latent variable state $h_{μ} (t)$ is a vector in the latent space. Because coupling to the latent variables is random throughout the population, only the length ( $\sim \sqrt{N_{F}}$ ) and not the direction of this vector matters, and the timescale of changes in this length would be much slower than $τ_{F}$ , the timescale of each of the components $h_{μ} (t)$ . We therefore speculated that increasing the timescale of dynamics of the latent variables should eventually lead to scaling and crackling in the single-variable model as well as the five-variable one. To examine the dependence of avalanche scaling on this timescale, we simulated one-variable and five-variable networks at fixed $η$ and $ϵ$ coupled to latent variables with the correlation time of their Ornstein-Uhlenbeck dynamics of $τ_{F} \in [10^{3}, 10^{5}]$ time steps, spanning from a factor of 10 faster to a factor of 10 slower than the original $τ_{F}$ in Figure 1. Simulations were replicated five times at each combination of parameters by drawing new latent variable coupling values ( $J_{i μ}$ ), as well as new latent variable dynamics and instances of neural firing. For simulations that passed the criteria to be fitted by power laws, we plot the fitted values of $τ$ , $α$ , $γ_{fit}$ , and $γ_{fit} - γ_{pred}$ (Figure 2G–J). Missing points are those for which distributions did not pass the power law fit criteria.

In the single-variable model, best-fit exponents changed abruptly for latent variable timescale around $τ_{F} = 10^{4}$ (Figure 2G and H), while in the five-variable model, exponents tended to increase gradually (Figure 2G and H, red). The discontinuity in the single-variable case reflected a change in the lower cutoff values in the power-law fits: size and duration distributions generated with faster latent dynamics could be fit reasonably well to a power law by using a high value of the lower cutoff (Figure 2—figure supplement 3). For time scales greater than ∼10⁴, the minimum cutoffs dropped, and the single-variable model generated power-law distributed avalanches and crackling (Figure 2J), similar to the five-variable model. In summary, in the latent dynamical variable model, introducing multiple variables generated scaling at faster timescales. However, by slowing the timescale of the latent dynamics, the model generated signatures of critical avalanche scaling for both multi- and single-variable simulations.

Avalanche criticality, input scaling, and firing threshold

In the previous section, we found that a very slow single latent dynamical variable generated avalanche criticality in the simulation population. Thus, from now on, we simplify the model in order to characterize avalanche statistics across values of input scaling $η$ and firing threshold $ϵ$ . Specifically, we modeled a population of $N = 128$ neurons coupled to a single quasi-static latent variable. We simulated 10³ segments of 10⁴ steps each and drew a new value of the latent variable ( $h \sim N (0, 1)$ ) for each segment. Ten replicates of the simulation were generated at each of the combinations of $η$ and $ϵ$ (see Materials and methods).

Almost independent of $η$ and $ϵ$ , we found quality power law fits and crackling. Figure 3 shows the average (across $n = 10$ network realizations) of the exponents extracted from size ( $τ$ , Figure 3A) and duration ( $α$ , Figure 3C) distributions. At small firing threshold ( $ϵ = 2$ ), we do not observe scaling because the system is always active, and all avalanches merge into one. At large firing threshold $ϵ$ and low input scaling $η$ , we do not observe scaling because activity is so sparse that all avalanches are small. At intermediate values of the parameters, the simulations generated plausible scaling relationships in size and duration. The difference between $γ_{fit}$ and $γ_{pred}$ was typically less than 0.1 (Figure 4D–F), which was consistent with previously reported differences between fit and predicted exponents (Ma et al., 2019). Thus, there appears to be no need for fine-tuning to generate apparent scaling in this model, at least in the limit of (near) infinite observation time. Wherever $η$ and $ϵ$ generate avalanches, there are approximate power-law distributions and crackling.

Figure 3

Download asset Open asset

Exponents across network simulations for networks of $N = 128$ neurons.

Each parameter combination $η, ϵ$ was simulated for ten replicates, each time drawing randomly the couplings $J_{i}$ , the latent variable values, and the neural activities. Other parameters in Table 3. (A) Average across replicates for the size exponent $τ$ . (B) Scatter plot of $α$ vs. $τ$ for each network replicate for parameter combinations indicated in A. Linear relationships between $τ$ and $α$ , corresponding to the minimum and maximum values of $γ_{f i t}$ from panel E, are shown to guide the eye. (C) Same as A, for duration exponent $α$ . (D) Predicted exponent, $γ_{pred}$ , derived from A and C. (E) Value of $γ_{f i t}$ from fit to $\bar{S} \sim D^{γ}$ . (F) Difference between $γ_{pred}$ and $γ_{fit}$ .

Figure 4 with 1 supplement see all

Download asset Open asset

Avalanches in the latent dynamical variable model with a single quasistatic variable.

Parameters in Table 3. (A) Number of avalanches in simulations from Figure 3 as a function of the calculated probability of avalanches at fixed $η$ across values of $ϵ$ and latent variable $h$ . Line indicates equality. (B) Probability of avalanches with $η = 2$ across values of $ϵ$ and $h$ . The latent variable $h$ is normally distributed with mean 0 and variance 1. Where the distribution of $h$ overlaps with regions of high probability (black), avalanches occur. (C) Probability of avalanches at $ϵ = 8$ across values of $η$ and $h$ . Increasing $η$ narrows the range of $h$ that generates avalanches. (D) Probability of avalanches at $h = 0$ for a populations of 128 neurons (black line) and for a varying $ϵ$ . Size distributions corresponding to simulations marked by the green and orange crosses are in E, F. (E) Example of size distribution with $ϵ < ϵ_{0}$ (orange marker in D). Size cutoff is close to 100. (F) Example of size distribution with $ϵ > ϵ_{0}$ (green marker in D). Size cutoff is < 10.

To determine where avalanches occur, we derive the avalanche rate across values of the latent variable $h$ . In the quasi-static model, the probability of an avalanche initiation is the probability of a transition from the quiet to an active state. Because all neurons are conditionally independent, this is $P_{ava} = P_{silence} (1 - P_{silence})$ . Then the expected number of avalanches ${\hat{N}}_{ava}$ is obtained by integrating $P_{ava}$ over $h$ at each value of $η$ and $ϵ$ :

{\hat{N}}_{a v a} = \int P_{a v a} (ϵ, η, h; J_{i}, N) p (h) d h = \int \prod_{i} (\frac{1}{1 + e^{- η J_{i} h - ϵ}}) (1 - \prod_{i} (\frac{1}{1 + e^{- η J_{i} h - ϵ}})) p (h) d h,

where $p (h)$ is the standard normal distribution. This probability tracks the observed number of avalanches across simulations, Figure 4A.

To gain an intuition for the conditions under which avalanches occur, we show two slices of the avalanche probability, at fixed $η$ (Figure 4B) and at fixed $ϵ$ (Figure 4C). Black regions indicate where avalanches are likely to occur. If, for a given value of $ϵ$ and $η$ , there is no overlap between high avalanche probability regions and the distribution of $h$ , then there will be no avalanches. For large $ϵ$ , avalanches occur because neurons with large coupling to the latent variable ( $η | J_{i} | >> 1$ , recall $J_{i} \sim N (0, 1)$ ) are occasionally activated by a value of the latent variable $h$ that is sufficient to exceed $ϵ$ (Figure 4B). Thus, the scaling parameter $η$ controls the value of $h$ for which avalanches occur most frequently (Figure 4C). As $ϵ$ decreases, avalanches occur for smaller and smaller $h$ until avalanches primarily occur when $h = 0$ .

To calculate the probability of avalanches, we must integrate over all values of $h$ , but we can gain a qualitative understanding of which avalanche regime the system is in by examining the probability of avalanches at $h = 0$ . At $h = 0$ , the avalanche probability (see Materials and methods) is

P_{a v a} (ϵ, η, h = 0; J_{i}, N) = {(\frac{1}{1 + e^{- ϵ}})}^{N} (1 - {(\frac{1}{1 + e^{- ϵ}})}^{N}),

which is maximized at $ϵ_{0} = - \log (2^{1 / N} - 1)$ , independent of $J_{i}$ and $η$ . After some algebra, we find that $ϵ_{0} \sim \log N$ for large $N$ . The dependence on $N$ reflects that a larger threshold is required for larger networks: large networks ( $N \to \infty$ ) are unlikely to achieve complete network silence, therefore preventing avalanches from occurring. Similarly, small networks ( $N \sim 1$ ) are unlikely to fire consecutively and thus are unlikely to avalanche.

We plot $P_{ava} (ϵ, η; J_{i}, N, h = 0)$ as a function of $ϵ$ in Figure 4D. The peak at $ϵ_{0}$ divides the space into two regions. For $ϵ < ϵ_{0}$ , a power-law is only observed in the large-size avalanches, which are rare (Figure 4E, green). By contrast, when $ϵ > ϵ_{0}$ , minimum size cutoffs are low (Figure 4F, orange). Both regions, $ϵ < ϵ_{0}$ and $ϵ > ϵ_{0}$ , exhibit crackling noise scaling. If observation times are not sufficiently long (estimated in Figure 4—figure supplement 1), then scaling will not be observed in the $ϵ < ϵ_{0}$ region, whose scaling relations arise from rare events. Insufficient observation times may explain experiments and simulations where avalanche scaling was not found.

Inferring the latent variable

Our analysis of $P_{ava} (ϵ, η, h)$ at $h = 0$ suggested that there are two types of avalanche regimes: one with high activity and high minimum cutoffs in the power law fit (Type 1), and the other with lower activity and size cutoffs (Type 2). Further, when $P_{ava}$ drops to zero, avalanches disappear because the activity is too high or too low. We now examine how information about the value of the latent variables represented in the network activity relates to the activity type. To delineate these types, we calculated numerically $ϵ^{*} (η)$ , the value of $ϵ$ , for which the probability of avalanches is maximized, and the contours of $P_{ava}$ (Figure 5A). Curves for $ϵ^{*} (η)$ and $ϵ_{0}$ and $P_{ava} = 10^{- 3}$ are shown in Figure 5A and B.

Figure 5

Download asset Open asset

Information in the neural activity about the latent variable is higher in the low- $ϵ$ avalanche region, compared to high- $ϵ$ avalanche or high-rate avalanche-free activity.

(A) Probability of avalanche per time step across values of $η$ and $ϵ$ . Solid magenta curve follows $ϵ^{*} (η)$ , the value of $ϵ$ maximizing the probability of avalanches at fixed $η$ . Dashed magenta line indicates $ϵ_{0}$ , calculated analytically, which matches $ϵ^{*}$ at $η = 0$ . (B) Information about latent variable, calculated from maximum likelihood estimate of $h$ using population activity. MLE approximation is invalid in the dark-blue region bounded by gray curve. Magenta line marks the maximum values of $P_{a v a}$ , reproduced from A. Dashed black curve indicates $P_{a v a} = 0.001$ . The highest information region falls between $ϵ^{*} (η)$ and the contour for $P_{ava} = 0.001$ . (**C - E**) Slices of B, showing $I_{M L E} (ϵ)$ for $η = {2, 5, 9}$ . Magenta and dashed black lines again indicate $ϵ^{*}$ and $P_{ava} = 0.001$ , respectively, as in B. Black dashed line marks the approximate boundary between the high-activity/no avalanche and the high-cutoff avalanche, and magenta line marks boundary between high-cutoff and low-cutoff avalanche regions.

We expect that the more cells fire, the more information they would convey, until the firing rate saturates, and inferring the value of the latent variable becomes impossible. Figure 5B supports the prediction: generally, information is higher in regions with more activity (lower $ϵ$ , higher $η$ ), but only up to a limit: as $ϵ \to 0$ , information decreases. This decrease begins approximately where the probability of avalanches drops to nearly zero (dashed black lines, Figure 5B–E) because all of the activity merges into a few very large avalanches. In other words, the Type-1 avalanche region coincides with the highest information about the latent variable.

The critical brain hypothesis suggests that the brain operates in a critical state, and its functional role may be in optimizing information processing (Beggs, 2008; Chialvo, 2010). Under this hypothesis, we would expect the information conveyed by the network to be maximized in the regions we observe avalanche criticality. However, we see that critical regions do not always have optimal information transmission. In Figure 5, the region that displays crackling noise is that where avalanches exist ( $P_{ava} > 0.001$ ), which corresponds to any $η$ value and $ϵ ≳ 3$ . This avalanche region encompasses both networks with high information transmission and networks with low information transmission. In summary, observing avalanche criticality in a system does not imply a high-information processing network state. However, the scaling can be seen at smaller cutoffs, and hence with shorter recordings, in the high-information state. This parallels the discussion by Schwab et al., 2014, who noticed that the Zipf’s law always emerges in neural populations driven by quasi-stationary latent fields, but it emerges at smaller system sizes when the information about the latent variable is high.

Discussion

Here, we studied systems with distributed, random coupling to latent dynamical variables and we found that avalanche criticality is nearly always observed, with no fine-tuning required. Avalanche criticality was surprisingly robust to changes in input gain and firing rate threshold. Loss of avalanche criticality could occur if the latent process was not well-sampled, either because the simulation was not long enough or the dynamics of the latent variables were too fast. Finally, while information about the latent variables in the network activity was higher where avalanches were generated compared to when they were not, there was a range of information values across the critical avalanche regime. Thus, avalanche criticality alone was not a predictor of optimal information transmission.

Explaining experimental exponents

A wide range of critical exponents has been found in ex vivo and in vivo recordings from various systems. For instance, the seminal work on avalanche statistics in cultured neuronal networks by Beggs and Plenz, 2003 found size and duration exponents of 1.5 and 2.0 respectively, along with $γ = 2$ , when time was discretized with a time bin equal to the average inter-event interval in the system. A subset of the in vivo recordings analyzed from anesthetized cat (Hahn et al., 2010) and macaque monkeys (Petermann et al., 2009) also exhibited a size distribution exponent close to 1.5. By contrast, a survey of many in vivo and ex vivo recordings found power-law size distributions with exponents ranging from 1 to 3 depending on the system (Fontenele et al., 2019). Separately, Ma et al., 2019 reported recordings in freely moving rats with size exponents ranging from 1.5 to 2.7. In these recordings, when the crackling relationship held, the reported value of $γ$ was near 1.2 (Fontenele et al., 2019; Ma et al., 2019).

A model for generating avalanche criticality is a critical branching process (Beggs and Plenz, 2003), which predicts size and duration exponents of 1.5 and 2 and scaling exponent $γ$ of 2. However, there are alternatives: Lombardi et al., 2023 showed that avalanche criticality may also be produced by an adaptive Ising model in the sub-critical regime, and in this case, the scaling exponent $γ$ was not 2 but close to 1.6. Our model, across the parameters we tested that produced exponents consistent with the scaling relationship, generated $τ$ values that ranged from 1.9 to about 2.5. Across those simulations, we found values $γ$ within a narrow band from 1.1 to 1.3 (see Figure 2I and J and Figure 3H). While the exponent values our model produces are inconsistent with a critical branching process ( $γ = 2$ ), they match the ranges of exponents estimated from experiments and reported by Fontenele et al., 2019. In this context, it might be useful to explore if our model and that of Lombardi et al., 2023 might be related, with the adaptive feedback signal of the latter viewed as an effective, latent variable of the former.

A genuine challenge in comparing exponents estimated from different experiments with different recording modalities (spiking activity, calcium imaging, LFP, EEG, or MEG) arises from differences in spatial and temporal scale specific to a particular recording, as well as the myriad decisions made in avalanche analysis, such as defining thresholds or binning in time. Thus, one possible reason for differences in exponents across studies may derive from how the system is sub-sampled in space or coarse-grained in time, both of which systematically change exponents $τ$ and $α$ (Beggs and Plenz, 2003; Shew et al., 2015) and could account for differences in $γ$ (Capek et al., 2023). The model we presented here could be used as a test bed for examining how specific analysis choices affect exponents estimated from recordings.

A second possible explanation for differences in exponents is that different experiments study similar, but distinct biological phenomena. For instance, networks that were cultured in vitro may differ from those that were not, whether they are in vivo or ex vivo (i.e. brain slices), and sensory-processing networks may have different dynamics from networks with different processing demands. It is possible that certain networks develop connections between neurons such that they truly do produce dynamics that approximate a critical branching process, while other networks have different structure and resulting dynamics and thus can be better understood as coupled neurons receiving feedback (Lombardi et al., 2023) or as a system coupled to latent dynamical variables. This is especially true in sensory systems, where coupling to (latent) external stimuli in a way that the neural activity can be used to infer the stimuli is the reason for the networks’ existence (Schwab et al., 2014).

Relationship to past modeling work

Our model is not the first to produce approximate power-law size and duration distributions for avalanches from a latent variable process (Touboul and Destexhe, 2017; Priesemann and Shriki, 2018). In particular, Priesemann and Shriki, 2018 derived the conditions, under which an inhomogeneous Poisson process could produce such approximate scaling. The basic idea is to generate a weighted sum of exponentially distributed event sizes, each of which are generated from a homogeneous Poisson process. How each process is weighted in this sum determines the approximate power-law exponent, allowing one to tune the system to obtain the critical values of 1.5 and 2. Interestingly, this model did not generate non-trivial scaling of size with duration ( $S \sim D^{γ}$ ). Instead, they found $γ = 1$ , not the predicted $γ = 2$ . Our results differ significantly, in that $γ$ was typically between 1.1 and 1.3 and it was nearly always close to the prediction from $α$ and $τ$ . We speculate that this is due to nonlinearity in the mapping from latent variable to spiking activity, as doubling the latent field $h$ does not double the population activity, but doubling the rate of a homogeneous Poisson process does double the expected spike count. As biological networks are likely to have such nonlinearities in their responses to common inputs, this scenario may be more applicable to certain kinds of recordings.

Summary

Latent variables – whether they are emergent from network dynamics (Clark et al., 2023; Sederberg and Nemenman, 2020) or derived from shared inputs – are ubiquitous in large-scale neural population recordings. This fact is reflected most directly in the relatively low-dimensional structure in large-scale population recordings (Stringer et al., 2019; Pandarinath et al., 2018; Nieh et al., 2021). We previously used a model based on this observation to examine signatures of neural criticality under a coarse-graining analysis and found that coarse-grained criticality is generated by systems driven by many latent variables (Morrell et al., 2021). Here, we showed that the same model also generates avalanche criticality, and that when information about the latent variables can be inferred from the network, avalanche criticality is also observed. Crucially, finding signatures of avalanche criticality required long observation times, such that the latent variable was well-sampled. Previous studies showed that Zipf’s law appears generically in systems coupled to a latent variable that changes slowly relative to the sampling time, and that the Zipf’s behavior is easier to observe in the higher information regime (Schwab et al., 2014; Aitchison et al., 2016). However, this also suggests that observation of either scaling at modest data set sizes indeed points to some fine-tuning — namely to the increase of the information in the individual neurons (and, since neurons in these models are conditionally independent, also in the entire network) about the value of the latent variables. In other words, one would expect a sensory part of the brain, if adapted to the statistics of the external stimuli, to exhibit all of these critical signatures at relatively modest data set sizes. In monocular deprivation experiments, when the activity in the visual cortex is transiently not adapted to its inputs, scaling disappears, at least for recordings of a typical duration, and is restored as the system adapts to the new stimulus (Ma et al., 2019). We speculate that the observed recovery of criticality by Ma et al., 2019 could be driven by neurons adapting to the reduced stimuli state, for instance, by adjusting $η$ (input scaling) and $ϵ$ (firing rate threshold).

Taken together, these results suggest that critical behavior in neural systems – whether based on the Zipf’s law, avalanches, or coarse-graining analysis – is expected whenever neural recordings exhibit some latent structure in population dynamics and this latent structure can be inferred from observations of the population activity.

Materials and methods

Simulation of dynamic latent variable model

Request a detailed protocol

We study a model from Morrell et al., 2021, originally constructed as a model of large populations of neurons in mouse hippocampus. In the original version of the model, neurons are non-interacting, receiving inputs reflective of place-field selectivity as well as input current arising from a random projection from a small number of latent dynamical variables, representing inputs shared across the population of neurons that are not directly measured or controlled. In the current paper, we incorporate only the latent variables (no place variables), and we assume that every cell is coupled to every latent variable with some randomly drawn coupling strength.

The probability of observing a certain population state ${s_{i}}$ given latent variables $h_{μ} (t)$ at time $t$ is

P ({s_{i}} | {h_{μ}}) = \frac{1}{Z ({h_{μ}})} e^{- H ({s_{i}}, {h_{μ}})},

where $Z$ is the normalization, and $H$ is the ‘energy’:

H = \sum_{i, μ = 1}^{N, N_{f}} η h_{μ} (t) J_{i μ} s_{i} + ϵ s_{i} .

The latent variables $h_{μ} (t)$ are Ornstein-Uhlenbeck processes with zero mean, unit variance, and time constant $τ_{m}$ . Couplings $J_{i μ}$ are drawn from the standard normal distribution.

Parameters for each figure are laid out in Tables 1–3. For the infinite time constant simulation, we draw a value $h \sim N (0, 1)$ and simulate for 10000 time steps, then repeat for 1000 draws of $h$ .

Table 3

Simulation parameters for Figures 3 and 4.

Parameter	Description	Value
$ϵ$	bias towards silence	$ϵ \in {2, 4, . . .14}$
$η$	variance multiplier	$η \in {1, 2, . . .10}$
$N_{F}$	number of latent fields	$N_{F} = 1$
$τ_{F}$	latent field time constant	quasistatic
$N$	number of cells	$N = 128$

Time step units

Request a detailed protocol

Most results were presented using arbitrary time units: all times (i.e. $τ_{F}$ and avalanche duration $D$ ) are measured in units of an unspecified time step. Specifying a time bin converts the probability of firing into actual firing rates, in spikes per second, and this choice determines which part of the $η$ - $ϵ$ phase space is most relevant to a given experiment.

The time step is the temporal resolution at which activity is discretized, which varies from several to hundreds of milliseconds across different experimental studies (Beggs and Plenz, 2003; Fontenele et al., 2019; Ma et al., 2019). In physical units and assuming a bin size of 3 ms to 10 ms, our choice of $η$ and $ϵ$ in Figure 2 would yield physiologically realistic firing rate ranges (Hengen et al., 2016), with high-firing neurons reaching averages rates of 20-50 spikes/second and median firing-rate neurons around 1-2 spikes/second. The timescales of latent variables examined range from about 3 s to 3000 s, assuming 3-ms bins. Inputs with such timescales may arise from external sources, such as sensory stimuli, or from internal sources, such as changes in physiological state.

Simulations were carried out for the same number of time steps (2×10⁶), which would be approximately 1 to 2 ‘hours’, a reasonable duration for in vivo neural recordings. Note that at large values of $τ_{F}$ , the latent variable space is not well sampled during this time period.

Analysis of avalanche statistics

Setting the threshold for observing avalanches

Request a detailed protocol

In our model, we count avalanches as periods of continuous activity (in any subset of neurons) that is book-ended by time bins with no activity in the entire simulated neural network. For real neural populations of modest size, this method fails because there are no periods of quiescence. The typical solution is to set a threshold, and to only count avalanches when the population activity exceeds that threshold, with the hope that results are relatively robust to that choice. In our model, this operation is equivalent to changing $ϵ$ , which shifts the probability of firing up or down by a constant amount across all cells independent of inputs. Our results in Figure 3 show that $α$ and $τ$ decrease as the threshold for detection is increased (equivalent to large $| ϵ |$ ), but that the scaling relationship is maintained. The model predicts that $γ_{p r e d} - γ_{f i t}$ would initially increase slightly with the detection threshold before decreasing back to near zero.

Following the algorithm laid out in Clauset et al., 2009, we fit power laws to the size and duration distributions from simulations generating avalanches. We use least-squares fitting to estimate $γ_{f i t}$ , the scaling exponent for size with duration, assessing the consistency of the fit across decades.

Reading power laws from data

Request a detailed protocol

We want, from each simulation, a quantification of the quality of scaling (how many decades, minimally) and an estimate of the scaling exponents ( $τ$ for the size distribution, $α$ for the duration distribution). We first compile all avalanches observed in the simulation and for each avalanche, calculate its size (total activity across the population during the avalanche) and its duration (number of time bins). Following the steps outlined by Clauset et al., 2009, we use the maximum-likelihood estimator to determine the scaling exponent. This is the solution to the transcendental equation

\frac{ζ^{'} (\hat{α}, x_{m i n})}{ζ (\hat{α}, x_{m i n})} = - \frac{1}{n} \sum_{i = 1}^{n} \ln x_{i}

where $ζ (α, x_{min})$ is the Hurwitz zeta function and $x_{i}$ are observations; that is, either the size or the duration of each avalanche i. For values of $x_{min} < 6$ , a numerical look-up table based on the built-in Hurwitz zeta function in the symbolic math toolbox was used (MATLAB2019b). For $x_{min} > 6$ we use an approximation (Clauset et al., 2009),

\hat{α} = 1 + n {(\sum_{i} \ln \frac{x_{i}}{x_{m i n} - \frac{1}{2}})}^{- 1} .

To determine $x_{min}$ , we computed the maximum absolute difference between the empirical cumulative density ( $S (x)$ ) function and model’s cumulative density function $P (x)$ (the Kolmogorov-Smirnov (KS) statistic; $D = {max}_{x \geq x_{min}} | S (x) - P (x) |$ ). The KS statistic was computed between for power-law models with scaling parameter $\hat{α}$ and cutoffs $x_{min}$ . The value of $x_{min}$ that minimizes the KS statistic was chosen. Occasionally the KS statistic had two local minima (as in Figure 2—figure supplement 1), indicating two different power-laws. In these cases, the minimum size and duration cutoffs were the smallest values that were within 10% of the absolute minimum of the KS statistic. Note that the statistic is computed for each model only on the power-law portion of the CDF (i.e. $x_{i} \geq x_{min}$ ). We do not attempt to determine an upper cut-off value.

To assess the quality of the power-law fit, Clauset et al., 2009 compared the empirical observations to surrogate data generated from a semi-parametric power-law model. The semi-parametric model sets the value of the CDF equal to the empirical CDF values up to $x = x_{min}$ and then according to the power-law model for $x > x_{min}$ . If the KS statistic for the real data (relative to its fitted model) is within the distribution of the KS statistics for surrogate datasets relative to their respective fitted models, the power-law model was considered a reasonable fit.

Strict application of this methodology could give misleading results. Much of this is due to the loss of statistical power when the minimum cutoff is so high that the number of observations drops. For instance, in the simulations shown in Figure 2, the one-variable duration distribution passed the Clauset et al., 2009 criterion, with a minimum KS statistic of 0.03 when the duration cutoff was 18 time steps. However, for the five-variable simulation in Figure 2, a power-law would be narrowly rejected for both size and duration, despite having much smaller KS statistics: for $τ$ , the KS statistic was 0.0087 (simulation range: 0.0008 to 0.0082; number of avalanches observed: 58,787) and for $α$ it was 0.0084 (simulation range: 0.0011 to 0.0075). Below we discuss this problem in more detail.

Determining range over which avalanche size scales with duration

Request a detailed protocol

For fitting $γ$ , our aim was to find the longest sampled range, over which we have apparent power-law scaling of size with duration. Because our sampled duration values have linear spacing, error estimates are skewed if a naive goodness of fit criterion is used. We devised the following algorithm. First, the simulation must have at least one avalanche of size 500. We fit $S = c D^{γ}$ over one decade at a time. We chose as the lower duration cutoff the value of minimum duration, for which the largest number of subsequent (longer-duration) fits produced consistent fit parameters (Figure 2—figure supplements 3 and 4, top row). Next, with the minimum duration set, we gradually increased the maximum duration cut-off, and we determined whether there was a significant bias in the residual over the first decade of the fit. We selected the highest duration cutoff, for which there was no bias. Finally, over this range, we re-fit the power law relationship and extracted confidence intervals.

Our analysis focused on finding the apparent power-law relationship that held over the largest log-scale range. A common feature across simulation parameters ( $τ_{F}$ , $N_{F}$ ) was the existence of two distinct power-law regimes. This is apparent in Figure 2I, which shows that when $N_{F} = 1$ at small $τ_{F}$ , the best-fit $γ$ (that showing the largest range with power-law-consistent scaling) is much larger ( 1.7), and then above $τ_{F} \sim 3000$ , the best-fit $γ$ drops to around 1.3.

Statistical power of power-law tests

Request a detailed protocol

In several cases, we found examples of power-law fits that passed the rejection criteria commonly used to determine avalanche scaling relationships because of limited number of observations. A key example is that of the single latent variable simulation shown in Figure 2B, where we could not reject a power law for the duration distribution. Conversely, strict application of the surrogate criteria would reject a power law for distributions that were quantitatively much closer to a power-law (i.e. lower KS statistic), but for which we had many more observations and thus a much stronger surrogate test (Figure 2). This points to the difficulty of applying a single criterion to determining a power-law fit. In this work, we adhere to the criteria set forth in Clauset et al., 2009, with a modification to control for the unreasonably high statistical power of simulated data. Specifically, the number of avalanches used for fitting and for surrogate analysis was capped at 500,000, drawn randomly from the entire pool of avalanches.

Additionally, we found examples, in which a short simulation was rejected, but increasing the simulation time by a factor of five yielded excellent power-law fits. We speculate that this arises due to insufficient sampling of the latent space. These observations raise an important biological point. Simulations provide the luxury of assuming the network is unchanging for as long as the simulator cares to keep drawing samples. In a biological network, this is not the case. Over the course of hours, the effective latent degrees of freedom could change drastically (e.g. due to circadian effects [Aton et al., 2009], changes in behavioral state [Fu et al., 2014], plasticity [Hooks and Chen, 2020], etc.), and the network itself (synaptic scaling, firing thresholds, etc.) could be plastic (Hengen et al., 2016). All these factors can be modeled in our framework by determining appropriate cutoffs (in duration of recording, in time step sizes, for activity distributions) based on specific experimental timescales.

Calculation of avalanche regimes

Request a detailed protocol

In the quasistatic model, we derive the dependence of the avalanche rate on $η$ , $ϵ$ and number of neurons $N$ , finding that there are two distinct regimes, in which avalanches occur. Each time bin is independent, conditioned on the value of $h$ . For an avalanche to occur, the probability of silence in the population (i.e. all $s_{i} = 0$ ) must not be too close to 0 (or there are no breaks in activity) or too close to 1 (or there is no activity). At fixed $h$ , the probability of silence is

P_{silence} (ϵ, η; J_{i}, N, h) = \prod_{i} \frac{1}{1 + \exp (- η J_{i} h + ϵ)} .

An avalanche occurs when a silent time bin is followed by an active bin, which has probability $P_{ava} (ϵ, η; J_{i}, N, h) = P_{silence} (1 - P_{% silence})$ .

Information calculation

Maximum-likelihood decoding

Request a detailed protocol

For large populations coupled to a single latent variable, we estimated the information between population spiking activity and the latent variable as the information between the maximum-likelihood estimator $h^{*}$ of the latent variable $h$ and the latent variable itself. This approximation fails at extremes of network activity levels (low or high).

Specifically, we approximated the log-likelihood of $h^{*}$ given $h_{true}$ near $h^{*}$ by $\log L (h - h^{*}) \approx \log L_{m a x} - \frac{1}{2} \frac{(h - h^{*})^{2}}{σ_{h^{*}}^{2}}$ . Thus we assume that $h^{*}$ is normally distributed about $h_{true}$ with variance $σ^{2} (h_{true})$ . The variance is then derived from the curvature of the log-likelihood at the maximum. The information between two Gaussian variables, here $P (h^{*} | h) = N (h, σ_{h^{*}}^{2})$ and $p (h) = N (0, 1)$ , is

I (h; {\bar{s}}_{i, T}) \approx \frac{1}{2} {⟨ \log \frac{T}{σ_{h_{t r u e}}^{2}} ⟩}_{h_{t r u e}},

where the average is taken over $h_{true} \sim N (0, 1)$ .

Given a set of $T$ observations of the neurons ${s_{i}}$ , the likelihood is

P ({s_{i}}_{t} | h) = \prod_{i, t}^{N, T} P (s_{i} | h) = \prod_{i, t}^{N, T} \frac{e^{- η s_{i} J_{i} h - ϵ s_{i}}}{1 + e^{- (η J_{i} h + ϵ)}} .

Maximizing the log likelihood gives the following condition:

0 = \frac{\partial (\log P)}{\partial h} |_{h^{*}} = \frac{\partial}{\partial h} (\sum_{i, t} ((- η s_{i} J_{i} h - ϵ s_{i}) - \log (1 + e^{- (η J_{i} h + ϵ)}))) |_{h^{*}}

= \sum_{i} - η {\bar{s}}_{i} J_{i} T + \frac{T J_{i} η}{1 + e^{η J_{i} h^{*} + ϵ}},

where ${\bar{s}}_{i} = \frac{1}{T} \sum_{t} s_{i t}$ is the average over observations $t$ . The uncertainty in $h^{*}$ is $σ_{h}$ , which was calculated from the second derivative of the log likelihood:

\frac{1}{σ_{h^{*}}^{2}} = - \frac{\partial^{2} (\log P)}{\partial h^{2}}

= - \frac{\partial}{\partial h} (\sum_{i} - η {\bar{s}}_{i} J_{i} T + \frac{T J_{i} η}{1 + e^{η J_{i} h + ϵ}}) |_{h^{*}}

= \sum_{i} \frac{T (η J_{i})^{2} e^{η J_{i} h^{*} + ϵ}}{(1 + e^{η J_{i} h^{*} + ϵ})^{2}}

= \sum_{i} \frac{T (η J_{i})^{2}}{4 \cosh^{2} (\frac{η J_{i} h^{*} + ϵ}{2})} .

This expression depends on the observations ${\bar{s}}_{i}$ only through the maximum-likelihood estimate $h^{*}$ . When $h^{*} \to h_{true}$ , then the variance is

\frac{1}{σ_{h^{*}}^{2}} = \sum_{i} \frac{T (η J_{i})^{2}}{4 \cosh^{2} (\frac{η J_{i} h_{t r u e} + ϵ}{2})} \equiv \frac{T}{σ_{h_{t r u e}}^{2}} .

To generate Figure 5, we evaluated Equation 10 using Equation 18.

Code availability

Request a detailed protocol

Simulation code was adapted from our previous work (Morrell et al., 2021). Code to run simulations and perform analyses presented in this paper is uploaded as Source code 1 and also available from https://github.com/ajsederberg/avalanche (copy archived at Sederberg, 2024).

Data availability

The current manuscript is a computational study, so no data were generated for this manuscript. Modelling code is uploaded as Source code 1.

References

1. Ahrens MB
2. Li JM
3. Orger MB
4. Robson DN
5. Schier AF
6. Engert F
7. Portugues R
(2012) Brain-wide neuronal dynamics during motor adaptation in zebrafish
Nature 485:471–477.

https://doi.org/10.1038/nature11057
- PubMed
- Google Scholar
(2016) Zipf’s law arises naturally when there are underlying, unobserved variables
PLOS Computational Biology 12:e1005110.

https://doi.org/10.1371/journal.pcbi.1005110
- PubMed
- Google Scholar
1. Aton SJ
2. Seibt J
3. Dumoulin M
4. Jha SK
5. Steinmetz N
6. Coleman T
7. Naidoo N
8. Frank MG
(2009) Mechanisms of sleep-dependent consolidation of cortical plasticity
Neuron 61:454–466.

https://doi.org/10.1016/j.neuron.2009.01.007
- PubMed
- Google Scholar
1. Beggs JM
2. Plenz D
(2003) Neuronal avalanches in neocortical circuits
The Journal of Neuroscience 23:11167–11177.

https://doi.org/10.1523/JNEUROSCI.23-35-11167.2003
- PubMed
- Google Scholar
1. Beggs JM
(2008) The criticality hypothesis: how local cortical networks might optimize information processing
Philosophical Transactions of the Royal Society A 366:329–343.

https://doi.org/10.1098/rsta.2007.2092
- Google Scholar
1. Capek E
2. Ribeiro TL
3. Kells P
4. Srinivasan K
5. Miller SR
6. Geist E
7. Victor M
8. Vakili A
9. Pajevic S
10. Chialvo DR
11. Plenz D
(2023) Parabolic avalanche scaling in the synchronization of cortical cell assemblies
Nature Communications 14:2555.

https://doi.org/10.1038/s41467-023-37976-x
- PubMed
- Google Scholar
1. Chialvo DR
(2010) Emergent complex neural dynamics
Nature Physics 6:744–750.

https://doi.org/10.1038/nphys1803
- Google Scholar
(2023) Dimension of activity in random neural networks
Physical Review Letters 131:118401.

https://doi.org/10.1103/PhysRevLett.131.118401
- PubMed
- Google Scholar
(2009) Power-law distributions in empirical data
SIAM Review 51:661–703.

https://doi.org/10.1137/070710111
- Google Scholar
1. Dahmen D
2. Grün S
3. Diesmann M
4. Helias M
(2019) Second type of criticality in the brain uncovers rich multiple-neuron dynamics
PNAS 116:13051–13060.

https://doi.org/10.1073/pnas.1818972116
- Google Scholar
(2019) Criticality between Cortical States
Physical Review Letters 122:208101.

https://doi.org/10.1103/PhysRevLett.122.208101
- PubMed
- Google Scholar
1. Fu Y
2. Tucciarone JM
3. Espinosa JS
4. Sheng N
5. Darcy DP
6. Nicoll RA
7. Huang ZJ
8. Stryker MP
(2014) A cortical circuit for gain control by behavioral state
Cell 156:1139–1152.

https://doi.org/10.1016/j.cell.2014.01.050
- PubMed
- Google Scholar
1. Girardi-Schappo M
(2021) Brain criticality beyond avalanches: open problems and how to approach them
Journal of Physics 2:031003.

https://doi.org/10.1088/2632-072X/ac2071
- Google Scholar
1. Hahn G
2. Petermann T
3. Havenith MN
4. Yu S
5. Singer W
6. Plenz D
7. Nikolic D
(2010) Neuronal avalanches in spontaneous activity in vivo
Journal of Neurophysiology 104:3312–3322.

https://doi.org/10.1152/jn.00953.2009
- PubMed
- Google Scholar
(2016) Neuronal firing rate homeostasis is inhibited by sleep and promoted by wake
Cell 165:180–191.

https://doi.org/10.1016/j.cell.2016.01.046
- PubMed
- Google Scholar
1. Hooks BM
2. Chen C
(2020) Circuitry underlying experience-dependent plasticity in the mouse visual system
Neuron 106:21–36.

https://doi.org/10.1016/j.neuron.2020.01.031
- PubMed
- Google Scholar
1. Humplik J
2. Tkačik G
(2017) Probabilistic models for neural populations that naturally capture global coupling and criticality
PLOS Computational Biology 13:e1005763.

https://doi.org/10.1371/journal.pcbi.1005763
- PubMed
- Google Scholar
(2023) Statistical modeling of adaptive neural networks explains co-existence of avalanches and oscillations in resting human brain
Nature Computational Science 3:254–263.

https://doi.org/10.1038/s43588-023-00410-9
- PubMed
- Google Scholar
1. Ma Z
2. Turrigiano GG
3. Wessel R
4. Hengen KB
(2019) Cortical circuit dynamics are homeostatically tuned to criticality in vivo
Neuron 104:655–664.

https://doi.org/10.1016/j.neuron.2019.08.031
- PubMed
- Google Scholar
(2013) Context-dependent computation by recurrent dynamics in prefrontal cortex
Nature 503:78–84.

https://doi.org/10.1038/nature12742
- PubMed
- Google Scholar
1. Mariani B
2. Nicoletti G
3. Bisio M
4. Maschietto M
5. Oboe R
6. Leparulo A
7. Suweis S
8. Vassanelli S
(2021) Neuronal avalanches across the rat somatosensory barrel cortex and the effect of single whisker stimulation
Frontiers in Systems Neuroscience 15:709677.

https://doi.org/10.3389/fnsys.2021.709677
- PubMed
- Google Scholar
1. Mazor O
2. Laurent G
(2005) Transient dynamics versus fixed points in odor representations by locust antennal lobe projection neurons
Neuron 48:661–673.

https://doi.org/10.1016/j.neuron.2005.09.032
- PubMed
- Google Scholar
1. Meshulam L
2. Gauthier JL
3. Brody CD
4. Tank DW
5. Bialek W
(2019) Coarse graining, fixed points, and scaling in a large population of neurons
Physical Review Letters 123:178103.

https://doi.org/10.1103/PhysRevLett.123.178103
- PubMed
- Google Scholar
1. Mora T
2. Bialek W
(2011) Are biological systems poised at criticality?
Journal of Statistical Physics 144:268–302.

https://doi.org/10.1007/s10955-011-0229-4
- Google Scholar
(2023) Quasiuniversal scaling in mouse-brain neuronal activity stems from edge-of-instability critical dynamics
PNAS 120:e2208998120.

https://doi.org/10.1073/pnas.2208998120
- PubMed
- Google Scholar
(2021) Latent dynamical variables produce signatures of spatiotemporal criticality in large biological systems
Physical Review Letters 126:118302.

https://doi.org/10.1103/PhysRevLett.126.118302
- PubMed
- Google Scholar
1. Nieh EH
2. Schottdorf M
3. Freeman NW
4. Low RJ
5. Lewallen S
6. Koay SA
7. Pinto L
8. Gauthier JL
9. Brody CD
10. Tank DW
(2021) Geometry of abstract learned knowledge in the hippocampus
Nature 595:80–84.

https://doi.org/10.1038/s41586-021-03652-7
- PubMed
- Google Scholar
1. O’Byrne J
2. Jerbi K
(2022) How critical is brain criticality?
Trends in Neurosciences 45:820–837.

https://doi.org/10.1016/j.tins.2022.08.007
- PubMed
- Google Scholar
1. Pandarinath C
2. O’Shea DJ
3. Collins J
4. Jozefowicz R
5. Stavisky SD
6. Kao JC
7. Trautmann EM
8. Kaufman MT
9. Ryu SI
10. Hochberg LR
11. Henderson JM
12. Shenoy KV
13. Abbott LF
14. Sussillo D
(2018) Inferring single-trial neural population dynamics using sequential auto-encoders
Nature Methods 15:805–815.

https://doi.org/10.1038/s41592-018-0109-9
- PubMed
- Google Scholar
(1995) Avalanches, Barkhausen noise, and plain old criticality
Physical Review Letters 75:4528–4531.

https://doi.org/10.1103/PhysRevLett.75.4528
- PubMed
- Google Scholar
(2009) Spontaneous cortical activity in awake monkeys composed of neuronal avalanches
PNAS 106:15921–15926.

https://doi.org/10.1073/pnas.0904089106
- Google Scholar
1. Plenz D
2. Ribeiro TL
3. Miller SR
4. Kells PA
5. Vakili A
6. Capek EL
(2021) Self-organized criticality in the brain
Frontiers in Physics 9:639389.

https://doi.org/10.3389/fphy.2021.639389
- Google Scholar
(2008) Avalanche dynamics of human brain oscillations: relation to critical branching processes and temporal correlations
Human Brain Mapping 29:770–777.

https://doi.org/10.1002/hbm.20590
- PubMed
- Google Scholar
1. Ponce-Alvarez A
2. Jouary A
3. Privat M
4. Deco G
5. Sumbre G
(2018) Whole-brain neuronal activity displays crackling noise dynamics
Neuron 100:1446–1459.

https://doi.org/10.1016/j.neuron.2018.10.045
- PubMed
- Google Scholar
1. Priesemann V
2. Shriki O
(2018) Can a time varying external drive give rise to apparent criticality in neural systems?
PLOS Computational Biology 14:e1006081.

https://doi.org/10.1371/journal.pcbi.1006081
- PubMed
- Google Scholar
(2014) Zipf’s law and criticality in multivariate data without fine-tuning
Physical Review Letters 113:068102.

https://doi.org/10.1103/PhysRevLett.113.068102
- PubMed
- Google Scholar
1. Sederberg A
2. Nemenman I
(2020) Randomly connected networks generate emergent selectivity and predict decoding properties of large populations of neurons
PLOS Computational Biology 16:e1007875.

https://doi.org/10.1371/journal.pcbi.1007875
- PubMed
- Google Scholar
Software
1. Sederberg A
(2024) Avalanche, version swh:1:rev:c2e91a5341d1aa5978650cb2c22227c7e52997dc
Software Heritage.

https://archive.softwareheritage.org/swh:1:dir:de294ef532e0e07507488df691f10a229fd29665;origin=https://github.com/ajsederberg/avalanche;visit=swh:1:snp:83b5ae292da441d80eb6e99588c5dd3ac359a2fd;anchor=swh:1:rev:c2e91a5341d1aa5978650cb2c22227c7e52997dc
1. Shew WL
2. Clawson WP
3. Pobst J
4. Karimipanah Y
5. Wright NC
6. Wessel R
(2015) Adaptation to sensory input tunes visual cortex to criticality
Nature Physics 11:659–663.

https://doi.org/10.1038/nphys3370
- Google Scholar
(2019) Spontaneous behaviors drive multidimensional, brainwide activity
Science 364:255.

https://doi.org/10.1126/science.aav7893
- PubMed
- Google Scholar
1. Touboul J
2. Destexhe A
(2017) Power-law statistics and universal scaling in the absence of criticality
Physical Review. E 95:012413.

https://doi.org/10.1103/PhysRevE.95.012413
- PubMed
- Google Scholar

Article and author information

Author details

Mia C Morrell

Department of Physics, New York University, New York, United States

Contribution
Conceptualization, Software, Formal analysis, Validation, Investigation, Visualization, Methodology, Writing – original draft, Writing – review and editing

Competing interests
No competing interests declared
Ilya Nemenman

Department of Physics, Department of Biology, Initiative in Theory and Modeling of Living Systems, Emory University, Atlanta, United States

Contribution
Conceptualization, Software, Formal analysis, Validation, Investigation, Visualization, Methodology, Writing – review and editing

Competing interests
No competing interests declared
Audrey Sederberg

Department of Neuroscience, University of Minnesota Medical School, Minneapolis, United States

Present address
School of Psychology and School of Physics, Georgia Institute of Technology, Atlanta, United States

Contribution
Conceptualization, Software, Formal analysis, Validation, Investigation, Visualization, Methodology, Writing – original draft, Writing – review and editing

For correspondence
audrey.sederberg@gatech.edu

Competing interests
Reviewing editor, eLife

"This ORCID iD identifies the author of this article:" 0000-0003-4458-3773

Funding

Simons Foundation (Simons-Emory Consortium on Motor Control)

Ilya Nemenman

National Science Foundation (BCS/1822677)

Ilya Nemenman

National Institute of Neurological Disorders and Stroke (2R01NS084844)

Ilya Nemenman

National Institute of Mental Health (1RF1MH130413)

Audrey Sederberg

Simons Foundation (Investigator Program)

Ilya Nemenman

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

IN was supported in part by the Simons Foundation Investigator program, the Simons-Emory Consortium on Motor Control, NSF grant BCS/1822677 and NIH grant 2R01NS084844. AS was supported in part by NIH grant 1RF1MH130413-01 and by startup funds from the University of Minnesota Medical School. The authors acknowledge the Minnesota Supercomputing Institute (MSI) at the University of Minnesota for providing resources that contributed to the research results reported within this paper.

Version history

Sent for peer review: May 30, 2023
Preprint posted: June 8, 2023
Reviewed Preprint version 1: September 12, 2023
Reviewed Preprint version 2: February 2, 2024
Version of Record published: March 12, 2024

Cite all versions

You can cite all versions using the DOI https://doi.org/10.7554/eLife.89337. This DOI represents all versions, and will always resolve to the latest one.

Copyright

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.