L4 · Chapter 7

Why a power law cannot prove the brain is critical

Sandbox Physics · 2026-10-31

The last chapter found that a branching process at its critical point gives avalanche sizes falling with exponent three halves. Experiments often argue the other way: the recorded avalanches look like a power law, so the brain must be critical. This chapter checks that argument. First it fits power laws properly; then, in simulations where the truth is known, it asks which tests can be fooled and by what; finally it looks at a preprint that offers another reason why fitted models so often sit near a critical point.

Fit the power law properly first

A straight line on log–log axes is not a test. Clauset, Shalizi and Newman give a three-step recipe. For each candidate lower cutoff $s_{\min}$ , estimate the exponent by maximum likelihood, and keep the cutoff that minimises the Kolmogorov–Smirnov distance $D$ between the data and the fitted law. Then draw many synthetic data sets from that law and count the fraction whose distance exceeds the measured one; that fraction is the goodness-of-fit $p$ . A value $p\ge0.1$ means the power law is not rejected, not that it is proved.

For discrete data the exponent has an approximate formula, given below, which they report accurate to about 1% once $s_{\min}$ is above about 6. We tested the Lab’s fit on synthetic data of known exponent, 40 sets of 3000 values each: a true exponent of 1.5 gives $1.4969\pm0.0020$ , 2 gives $1.988\pm0.006$ and 2.5 gives $2.434\pm0.015$ . The last is 4.5 standard errors low, because the median selected cutoff is only 3, below the range where the approximation holds.

The Lab’s $p$ takes one shortcut: the synthetic sets keep the cutoff chosen for the measured data, whereas the original method chooses a new cutoff for every synthetic set. Without that optimisation the synthetic distances tend to be larger, so the Lab’s $p$ is more lenient than the published one.

\hat\alpha\simeq1+n\Bigl[\sum_{i=1}^{n}\ln\frac{s_i}{s_{\min}-\frac12}\Bigr]^{-1},\qquad \sigma_{\hat\alpha}=\frac{\hat\alpha-1}{\sqrt n},\qquad D=\max_{s\ge s_{\min}}\bigl|S(s)-P(s)\bigr|

Left: the approximate maximum-likelihood exponent for discrete data, using the

n

events with

s_i\ge s_{\min}

. Middle: its statistical error. Right: the Kolmogorov–Smirnov distance between the cumulative distribution of the data

S

and the fitted law

P

Source: A. Clauset · J. Touboul and A. Destexhe

Power laws without critical cascades

Touboul and Destexhe found that thresholding a random signal, or shuffling a recording into a surrogate, can produce scaling much like neuronal avalanches; a straight line on a log–log plot does not survive testing. Priesemann and Shriki go further: a Poisson process with a constant rate gives no power law, but once the input rate switches among a few levels, avalanche sizes come close to a power law with exponent near 3/2 over a limited range. Morrell, Sederberg and Nemenman let independently firing neurons share a few slowly varying inputs and reproduced the scaling under coarse-graining that had been read as a sign of criticality.

The Lab has two such mimics, with no coupling between units at all. The one with a shared slow input, at the default array and a bin of 4 steps, gives $\hat\alpha=2.13$ and $p=0.87$ over 2.7 decades, and multistep regression gives 0.997. The one with a drifting input gives a multistep estimate of 1.000.

Multistep regression is fooled for a simple reason. It only measures how long activity stays correlated in time, and a slowly varying input keeps activity correlated for a long time without any unit triggering another. Priesemann and Shriki conclude that passive analysis cannot rule out an external drive that happens to mimic critical activity; small, controlled perturbations are needed.

Source: J. Touboul and A. Destexhe · V. Priesemann and O. Shriki · M. C. Morrell · J. Wilting and V. Priesemann

A false-positive table

To see what fools each test we took four models with a known truth: a coupled network near the critical point, a subcritical coupled network, and two uncoupled mimics, one with a drifting input and one with a shared slow input. Each model was analysed with 3 electrode arrays (covering 60.8%, 10% and 2% of the units), 3 bin widths (1, 4 and 8 steps) and 10 random seeds: 90 analyses per model.

A plausible power law ( $p\ge0.1$ over at least 1.5 decades): 0.178 for the near-critical network, 0.378 for the shared slow input, 0.022 for the drifting input and 0 for the subcritical network. The shared-slow-input mimic passes more than twice as often as the network that really is near critical. That network fails in 82% of analyses: all 60 with the 10% and 2% arrays fail, and with the largest array 3, 8 and 5 of 10 pass at bins of 1, 4 and 8 steps. So failing a power-law test does not show that a network is not critical either.

The naive verdict, a plausible power law with exponent within 0.2 of 1.5, passes only the near-critical network, at 0.178; the mimics’ plausible fits all have exponents between 1.75 and 2.74. A multistep estimate of at least 0.95 passes both mimics every time, the near-critical network 0.433 of the time and the subcritical one never.

The last test is an intervention: pause the external input, kick one recorded unit, and ask whether more than 2 recorded spikes follow on average. The near-critical network passes all 90 times, neither mimic ever does, and the subcritical network passes 0.667 of the time. A kick shows whether there is coupling; how close to critical the network is needs the number of spikes one kick triggers in the whole network.

Source: A. Clauset · V. Priesemann and O. Shriki · J. Wilting and V. Priesemann

Why fitted models so often sit near a critical point

Another line of evidence comes from maximum-entropy models: fit an Ising-type model using only each neuron’s mean activity and the pairwise correlations, and the parameters often land near a phase transition. A 2026 preprint by Carcamo and Lynn offers an explanation. In the Curie–Weiss model, which keeps only the mean activity $m$ and a correlation strength $\chi$ , a small patch of parameters near the critical point maps onto a large region of statistics, and that region grows with the number of neurons. Inferring parameters from data therefore pushes models toward the critical line without any fine-tuning.

We checked their Eq. (3) on the model only; we did not reanalyse any recording. Take $m=-0.8$ and $\chi=2$ , already above the mean-field upper limit of 0.712. Solving the inverse problem exactly gives $\lambda$ of 1.478, 1.374 and 1.371 for $N=20$ , 100 and 1000, all above the critical value 1, and $h$ of $-0.110$ , $-0.033$ and $-0.0048$ , approaching 0 as $N$ grows, that is, approaching the critical line. The double-well approximation of Eq. (3) gets $\lambda$ with an error falling from 1.8% to 0.29%.

Page 11 of the preprint says the inferred interactions exceed the critical strength but writes $\lambda\ge0$ in parentheses. By Eq. (3), $\lambda=\tanh^{-1}(m_0)/m_0$ is at least 1 for every $0<m_0<1$ , and the critical strength of this model is 1, so it should read $\lambda\ge1$ . The paper is still a preprint and has not been peer reviewed.

P(\mathbf x)\propto\exp\Bigl[h\sum_{i}x_i+\frac{\lambda}{2N}\Bigl(\sum_{i}x_i\Bigr)^{2}\Bigr],\qquad m_0=\sqrt{m^2+\chi/N},\quad h\simeq\frac{1}{N m_0}\tanh^{-1}\frac{m}{m_0},\quad \lambda\simeq\frac{1}{m_0}\tanh^{-1}m_0

Left: the Curie–Weiss model with

x_i=\pm1

, external field

h

and coupling strength

\lambda

. Right: the double-well approximation of the preprint’s Eq. (3), which infers the parameters from the mean activity

m

and the correlation strength

\chi

(

N

times the variance of the population mean activity).

Source: D. P. Carcamo and C. W. Lynn · M. C. Morrell

Worked example: test the Lab’s default record end to end

Problem. In the Lab’s default record the selected cutoff is $s_{\min}=2$ , $n=1990$ avalanches lie at or above it, $\sum\ln[s_i/(s_{\min}-\tfrac12)]=4306.3$ , the largest avalanche is 7461 and the goodness of fit is $p=0.21$ . Find the exponent and its error and decide whether the power law is plausible. Then repeat for the shared-slow-input mimic: $n=646$ , $s_{\min}=11$ , a sum of 572.08, largest 5401, $p=0.87$ . Finally draw a conclusion from the kicks.

Model. Use the approximate maximum likelihood and the error formula above; the number of decades is the base-ten logarithm of the largest value over the cutoff.

Numbers. Default record: $\hat\alpha=1+1990/4306.3=1.462$ , $\sigma_{\hat\alpha}=0.462/\sqrt{1990}=0.0104$ , spanning $\log_{10}(7461/2)=3.57$ decades with $p=0.21$ : a plausible power law, 0.04 from 1.5. Mimic: $\hat\alpha=1+646/572.08=2.129$ , $\sigma_{\hat\alpha}=0.0444$ , spanning $\log_{10}(5401/11)=2.69$ decades with $p=0.87$ : also a plausible power law.

Explanation. Both records pass the power-law test and differ only in the exponent. The kicks settle it. In the default network each kick yields 97 recorded spikes on average and 155 in the whole network, a cascade ratio of $1-1/155=0.9936$ . In the mimic each kick yields only the kicked spike. The first is a coupled cascade near the critical point; the second has no cascades at all.

Check. The two errors, 0.0104 and 0.0444, match the Lab readouts. Had the mimic’s exponent happened to be near 1.5, neither the power-law test nor the exponent could tell the two apart, and only the kick would remain.

Source: A. Clauset · V. Priesemann and O. Shriki

Try it in the Lab

A power law fools the fit, not the kick

Predict first: choose Shared slow input, no coupling. Will the power-law fit be rejected? What will multistep regression give?
Choose the preset and read the exponent, the p-value, the multistep estimate and the kick readout. Then go back to Cascades near the critical point and compare the four readouts.
Explain: the mimic’s power law is just as plausible and its multistep estimate is even closer to one, but a kick yields only the kicked spike. In the coupled network one kick triggers nearly a hundred recorded spikes.

Slow input breaks multistep regression

Predict first: in Drifting input, no coupling, the units have no connections at all. What will multistep regression give?
Choose the preset and read the multistep estimate. Then drag the analysis bin from 8 steps down to 1 and watch the exponent and the multistep estimate.
Explain: the multistep estimate is about 1.000 because the slow drift of the input keeps activity correlated for a long time. The exponent changes a lot with the bin width, which is itself a sign that there are no cascades.

Exercises with explained answers

A record has 500 avalanches above the cutoff and a fitted exponent of 1.5. What is the statistical error? Roughly how many avalanches would bring it down to 0.005?
Show the answer
$0.5/\sqrt{500}=0.022$ ; from $0.5/\sqrt n=0.005$ we get $n=10^4$ . This is the statistical error only; choosing the cutoff and the discrete approximation add systematic errors on top.
Why does keeping the measured cutoff for the synthetic sets make the goodness-of-fit p-value lenient?
Show the answer
The cutoff for the measured data was chosen to make its distance as small as possible. Synthetic sets that do not get the same choice tend to have larger distances, so a larger fraction of them exceed the measured distance and p comes out larger. The original method chooses a new cutoff for every synthetic set too, so both sides receive the same optimisation.
In the Lab’s subcritical network one kick triggers 4.84 spikes in the whole network on average. What cascade ratio does that give? Compare with the roughly 0.79 that multistep regression gives on all units.
Show the answer
$1-1/4.84=0.793$ , in line with 0.79. In the subcritical network a single cascade in a quiet network and the branching ratio during ongoing activity differ little; near the critical point collisions matter more and the two differ more (0.9936 against 0.9855).
Use Eq. (3) of the preprint to compute $\lambda$ and $h$ for $m=-0.8$ , $\chi=2$ and $N=1000$ .
Show the answer
$m_0=\sqrt{0.64+0.002}=0.8012$ , $\lambda=\tanh^{-1}(0.8012)/0.8012=1.3755$ , $h=\tanh^{-1}(-0.8/0.8012)/(1000\times0.8012)=-0.00447$ . The exact solution is 1.3715 and $-0.00476$ , so $\lambda$ differs by 0.29%.

Sources and what we recomputed

J. Wilting and V. Priesemann, Inferring collective dynamical states from widely unobserved systems, arXiv:1608.07035v3; Nature Communications 9, 2325 (2018)
Version 3, main text pp. 1 to 3; the supplementary methods were not read. The branching process with drive, the conventional estimator as the slope of the next count against the current count, the result that subsampling multiplies every lagged slope by the same factor, and the multistep regression estimator that fits that exponential across lags. Their examples: a process with m = 0.9 sampled at 10 % gives a conventional estimate of 0.312.
V. Priesemann and O. Shriki, Can a time varying external drive give rise to apparent criticality in neural systems?, PLoS Computational Biology 14, e1006081 (2018)
Introduction, Results and Discussion read; Methods not read. How avalanches are defined from binned events. A Poisson process with a constant rate gives no power law; one whose rate switches among a few levels gives sizes close to a power law with exponent near 3/2 over a limited range. The measures that tell the two apart, and their conclusion that passive analysis cannot rule out a drive that mimics critical activity, so a small controlled perturbation is needed.
A. Clauset, C. R. Shalizi and M. E. J. Newman, Power-law distributions in empirical data, arXiv:0706.1062v2; SIAM Review 51, 661–703 (2009)
Version 2, read in part: the abstract, Box 1, Section 3.1 on discrete maximum likelihood, Section 3.3 on choosing the lower cutoff by the Kolmogorov–Smirnov distance and the p-value parts of Section 4. Discrete approximation for the exponent with the half-integer shift; p above 0.1 means plausible, not proven; a large p does not show that the data follow a power law.
J. Touboul and A. Destexhe, Can power-law scaling and neuronal avalanches arise from stochastic dynamics?, arXiv:0910.0805v3; PLoS ONE 5, e8982 (2010)
Version 3, read in part: the abstract, introduction and part of the methods. A straight line on a log–log plot is not a robust test; thresholded stochastic signals and shuffled surrogates can show similar scaling.
M. C. Morrell, A. J. Sederberg and I. Nemenman, Latent dynamical variables produce signatures of spatiotemporal criticality in large biological systems, arXiv:2008.04435v1; Physical Review Letters 126, 118302 (2021)
Version 1, read the abstract and the model section. Neurons that fire independently given a few slowly varying shared inputs reproduce the scaling under coarse-graining that had been read as a sign of criticality.
D. P. Carcamo and C. W. Lynn, Emergence of criticality in models of real neurons, arXiv:2609.09438v1 (8 September 2026), preprint, not peer reviewed
Version 1, read in full. Models fitted to recordings of many neurons tend to sit near critical points because a small region of parameters maps onto a large region of possible statistics; the effect grows with the number of neurons. Their mean-activity model, Eq. (3), gives a double-well approximation for the fitted field and coupling. They report that it converges to the exact parameters for about 100 neurons or more. We solved the same inverse problem exactly and compared. On p. 11 the text says the inferred couplings exceed the critical strength and writes this as a coupling of at least 0; the critical strength of this model is 1, so we read it as at least 1.