Neural Criticality Tutorial
Loading tutorial…
From Physics to Implementation — Understanding Critical Neural Dynamics
Imagine you're in a dark room and someone whispers. If your brain is too quiet (subcritical), you won't hear it. If it's too noisy (supercritical), the whisper gets drowned out. But at just the right balance—the critical point—your brain can detect the faintest signals while avoiding chaos.
This tutorial explores criticality in neural networks. We'll learn:
You know that water freezes at 0°C. Above this temperature, water molecules move freely (liquid). Below it, they lock into a crystal structure (solid). The temperature of 0°C is a critical point—a boundary between two phases.
Imagine building a sand pile one grain at a time:
Pile is too flat. Grains just sit there. Nothing interesting happens.
Pile at perfect slope. Adding one grain can trigger avalanches of any size.
Pile too steep. Constant massive landslides. Complete chaos.
Every phase transition has an order parameter—a quantity that is zero in one phase and non-zero in the other:
| System | Order Parameter | Subcritical | Supercritical |
|---|---|---|---|
| Magnet | Magnetization | 0 (disordered) | >0 (aligned) |
| Water/Ice | Density diff. | 0 (liquid) | >0 (solid) |
| Neural Net | Firing Rate | 0 (silent) | >0 (active) |
The brain has two main types of neurons:
~80% of neurons. Release glutamate. Increase activity of targets.
~20% of neurons. Release GABA. Decrease activity of targets.
Each connection between neurons has a weight determining its strength:
This single parameter w₀ captures the net balance between excitation and inhibition:
Inhibition wins
Activity dies out. Poor signal transmission.
Critical balance
Optimal dynamics. Maximum information capacity.
Excitation wins
Runaway activity. Seizure-like behavior.
Beggs & Plenz (2003) discovered neuronal avalanches in cortical slices—cascades of neural activity with power-law distributed sizes. This was the first direct evidence that the brain might operate near criticality.
The original Wilson-Cowan model (1972) is deterministic—given the same initial conditions, it always produces the same output. But real neural networks have noisefrom finite numbers of neurons.
Each neuron can be active (firing, state = 1) or quiescent(silent, state = 0). Transitions occur stochastically with rates determined by two key parameters:
Depends on synaptic input s (no activation if s ≤ 0)
Constant rate (refractory period)
This is equivalent to si = Σj wij × aj + h from the original Wilson-Cowan formulation, but here we make each population explicit. In the mean-field approximation, each neuron receives the same average input from each population:
Where:
The mean-field assumption is that every neuron sees the same average input based on population activity fractions, rather than the specific pattern of which neurons are active.
This is exact for all-to-all connectivity with uniform weights. Forhomogeneous random networks (where each neuron has the same statistical connectivity pattern), the mean-field approximation becomes increasingly accurate as N → ∞ because each neuron samples from a large enough pool that fluctuations in the input average out.
In other words: as N grows, each neuron's input converges to the population mean, and the mean-field equations become exact in the thermodynamic limit.
The critical point w₀c is where the system transitions from absorbing (silent) to active. We can derive this analytically by linearizing the mean-field equations around the quiescent state.
The formula w₀c = α/β is analytically derived from the mean-field equations in the thermodynamic limit (N → ∞). It represents the exact bifurcation point where the quiescent fixed point loses stability.
However, for finite networks, fluctuations shift the effective critical point. To verify criticality at any specific N, one must empiricallysweep w₀ and check avalanche statistics. As N increases, the empirically-observed critical point converges toward 0.1.
The thermodynamic limit is when N → ∞. In this limit, the mean-field approximation becomes exact, fluctuations become negligible relative to average values, and the critical exponents (τS, τT, γ) attain their exact theoretical values (3/2, 2, 2).
For finite networks, we are not at the thermodynamic limit. Fluctuations shift the effective critical point—much like how a real circle can never perfectly achieve π for its circumference-to-diameter ratio, but approaches it in the limit of infinite precision. The finite-N deviation is called afinite-size effect.
To find the effective critical point for your network size:
Empirical observations of convergence toward the thermodynamic limit:
As N grows, the macroscopic properties (critical point, exponents) converge to their theoretical values—this convergence is the approach to the thermodynamic limit.
We discretize time into bins of width δ (typically ~1 ms). An avalanche is:
Each dot is a neuron firing. Red = excitatory, blue = inhibitory.
Total activations per time bin. Avalanches are bursts of activity separated by silent periods.
Total number of neurons that fired (0→1 transitions) during the avalanche.
Duration of each avalanche as measured by the number of bins.
| Regime | w₀ | Avalanche Behavior |
|---|---|---|
| Subcritical | < 0.1 | Small avalanches, activity dies quickly |
| Critical | ≈ 0.1 | Power-law distributed sizes |
| Supercritical | > 0.1 | Giant avalanches, system-spanning |
Most distributions you encounter have a characteristic scale. Human heights cluster around 170cm—you never see someone 10 meters tall. But some phenomena arescale-free: earthquakes, city populations, word frequencies. These followpower laws.
In words: the probability of observing a value x is proportional to x raised to the power −τ. The exponent τ (Greek letter "tau") characterizes the distribution.
Power laws are hard to see on linear plots because most data clusters near small values. Let's see this directly with our avalanche data:
Data clusters at small values—hard to see structure
Same data—now the linear relationship is clear!
Let's derive why this works by taking the logarithm of both sides:
This is the equation of a straight line with slope −τ! On a log-log plot, a power law appears as a linear relationship. The exponent τ is simply the negative of the slope.
Size exponent τS ≈ 1.5
Duration exponent τT ≈ 2.0
Scaling exponent γ ≈ 2.0
On these log-log plots, the probability of increasingly large avalanches fallslinearly with the size (or duration) of that avalanche. This linear relationship in log-log space is the defining signature of a power law.
Concretely: if τS = 1.5, then an avalanche 10× larger is 101.5 ≈ 32× less likely. This holds across all scales—from tiny ripples to massive cascades.
| Exponent | Distribution | Theory | Measured |
|---|---|---|---|
| τS | P(S) ∝ S−τS | 3/2 = 1.50 | 1.43 |
| τT | P(T) ∝ T−τT | 2.0 | 1.94 |
| γ | S ∝ Tγ | 2.0 | 2.19 |
Different physical systems can have identical critical exponents, even if their microscopic details are completely different. Systems with the same exponents belong to the same universality class.
| Class | Examples | τS | τT |
|---|---|---|---|
| Ising (2D) | Magnets, alloys | 1.06 | 1.5 |
| Dir. Percolation | Epidemics, fires | 1.11 | 1.45 |
| Mean-Field | Neurons, quakes | 1.50 | 2.00 |
The stochastic Wilson-Cowan model belongs to the mean-field branching processuniversality class. Here's where the exponents come from:
Think of each neuron activation as a "parent" that can produce "children" (downstream activations):
At criticality (branching ratio σ = 1), the probability theory of branching processes gives exact results:
These aren't adjustable parameters—they emerge from the mathematics of branching processes. Let's see exactly how:
"Mean-field" means each neuron feels the average effect of all others, not the specific pattern of of incoming spikes. This is valid when:
For sparse networks with specific topology, different universality classes (like directed percolation) may apply.
The three exponents are not independent! They are connected by a fundamental relationship called the scaling relation (or "crackling noise relation"):
This relation comes from dimensional analysis: if S ∝ Tγ and both S and T have power-law distributions, then the exponents must be related this way. It's aself-consistency check for criticality.
Rearranging, we can derive τS from τT and γ:
The scaling relation provides the most reliable estimate of τS because:
Comparison of methods:
| MLE: | 1.329 |
| Log-binning: | 1.534 |
| Scaling-implied: | 1.428 |
| Theory: | 1.500 |
If your system is truly critical, the exponents must satisfy this relation. If they don't:
We've collected avalanche data from our simulation—thousands of sizes and durations. To check if the system is at criticality, we need to extract the exponentsfrom this data and compare them to the theoretical values (1.5, 2.0, 2.0).
The challenge: we have raw data, but we need the exponent τ that describes the probability distribution P(x) ∝ x−τ. There are several methods to estimate τ.
The simplest approach: since log P(x) = −τ · log(x) + const, we can:
Pros: Intuitive, easy to visualize, works well for large datasets
Cons: Binning choices affect the result, biased for small samples, doesn't use all information in the data
A statistically rigorous approach that uses all data points, following Clauset et al. (2009):
Where n is the number of observations and xmin is the minimum value above which the power law holds.
Use the scaling relation with separately measured τT and γ:
| xmin | MLE | Log-Bin | Implied |
|---|---|---|---|
| 10 | 1.283 | 1.482 | 1.373 |
| 50 | 1.308 | 1.469 | 1.388 |
| 100 | 1.329 | 1.534 | 1.396 |
| 200 | 1.333 | 1.556 | 1.436 |
| 500 | 1.369 | 1.582 | 1.594 |
| 1000 | 1.399 | 1.602 | 1.800 |
| Theory | 1.500 | ||
A naive approach is to update all neurons at each timestep. But this has problems:
The Gillespie algorithm (1977) provides an exact simulation of continuous-time Markov processes:
There are only 4 types of events (using mean-field approximation):
| Event | Rate | Effect |
|---|---|---|
| E activates | (N − k) × f(s) | k → k + 1 |
| E deactivates | k × α | k → k − 1 |
| I activates | (N − l) × f(s) | l → l + 1 |
| I deactivates | l × α | l → l − 1 |
The time bin width δ defines avalanche boundaries. De Candia et al. used δ = 0.64 ms.
Power laws only hold above some minimum value xmin. The Clauset method finds xmin by minimizing the Kolmogorov-Smirnov distance between data and fit.
Rule of thumb: Start with xmin = 10-100 for sizes, and xmin = 5-10 for durations. Check that your exponent is stable as you vary xmin.
Small avalanches are most affected by finite-size effects and bin-size choices. The paper shows that for moderate avalanche sizes, measured exponents depend on the time bin width δ. However, for large avalanches (S > 10³), the scaling becomes independent of δ.
Why bin-size independence at criticality? At the critical point, the dynamics are scale-invariant (self-similar): there is no characteristic length or timescale. If you zoom in or out in time (which is what changing bin size does), the statistical properties should look the same. This only manifests clearly for large avalanches that span many timescales—small avalanches are too short to "see" the scale-free structure.
(Note: "Critical slowing down" is a related but distinct phenomenon—it refers specifically to the relaxation time diverging at criticality, meaning perturbations take longer and longer to decay. Scale invariance and critical slowing down both stem from the diverging correlation length, but they describe different aspects of critical behavior.)
This means: if your network is truly at criticality, increasing xminshould cause the measured exponents to converge toward theoretical values, not diverge. Large avalanches better approximate the true power-law behavior because they span many time bins and many neurons, averaging out finite-size artifacts.
Practical test: Plot your measured τS as a function of xmin. If you're at criticality, it should stabilize near 1.5 for large xmin. If it diverges or keeps changing, you may not be at the true critical point, or you may need more data.
For finite networks, the "effective" critical point shifts:
| N | w₀ | τS | τT | γ |
|---|---|---|---|---|
| 10,000 | 0.12 | 1.40 | 2.04 | 1.92 |
| 100,000 | 0.10 | 1.37 | 2.00 | 2.09 |
| 1,000,000 | 0.10 | 1.36 | 1.97 | 2.11 |
| ∞ | 0.10 | 1.50 | 2.00 | 2.00 |
An important clarification: w₀ = WE − WI is the difference in synaptic strengths, not a ratio of population sizes.
The model works with any E/I ratio (50/50, 80/20, etc.). The critical point w₀c = α/β = 0.1 is determined by the deactivation and activation rates, not by the number of E vs I neurons.