This simulator illustrates the core setup of the paper: how dynamical regimes at the micro level (inside each agent) relate to collective behavior at the macro level (across interacting agents).
Agents are placed in a 2D toroidal space and communicate only through binary light signals. Each agent contains an internal probabilistic reservoir; its light state is read from a designated neuron. At each macro step, agents receive egocentric exteroceptive observations of illuminated neighbors within vision radius $v_r$, and these inputs perturb their internal reservoir dynamics.
The key question is whether near-critical micro dynamics are sufficient to produce macro-level critical-like avalanche statistics, or whether collective behavior is primarily constrained by interaction topology and activity propagation.
Paper: Emergent Macro-Criticality from Micro-Critical Agents (arXiv:2605.01818). For quantitative claims and code, see the paper repository.
Each agent contains an independent reservoir of $N_{\mathrm{neurons}}$ binary neurons. Micro-dynamics are defined by the discrete-time branching-process update introduced by Beggs and Plenz [1] for neuronal avalanches. The neuron state is $s_{t,i}\in\{0,1\}$, and connectivity is a transmission-probability matrix $W_{i,j}\in[0,1]$. For each active pre-synaptic neuron, transmission events are sampled as $X_{t,i,j}\sim\mathrm{Bernoulli}(W_{i,j})$.
A neuron activates at $t+1$ if it receives at least one successful transmission:
$$s_{t+1,j}=\left[\sum_{i=1}^{N_{\mathrm{neurons}}} s_{t,i}X_{t,i,j}\ge1\right]$$
We remove self-connections ($W_{i,i}=0$) and use homogeneous transmission:
$$W_{i,j}=\begin{cases}p^{\mathrm{micro}}, & i\ne j\\0, & i=j\end{cases}$$
Under this assumption, the branching ratio $\sigma=(N_{\mathrm{neurons}}-1)p^{\mathrm{micro}}$ is the expected number of activations generated by one active neuron at the next micro-step. Intuitively: if $\sigma=1$, each active neuron replaces itself with about one new active neuron on average (activity is balanced); if $\sigma<1$, activity tends to shrink and die out; if $\sigma>1$, activity tends to expand and can become dense/saturated in finite systems. The equivalent critical point in terms of the transmission parameter is $p_c^{\mathrm{micro}}=\frac{1}{N_{\mathrm{neurons}}-1}$. Note: in this animation, the slider varies $\sigma$ directly; this is equivalent to varying $p^{\mathrm{micro}}$ through $\sigma=(N_{\mathrm{neurons}}-1)p^{\mathrm{micro}}$. The raster below illustrates these qualitative regimes.
What is this? A single-reservoir raster (rows = neurons, columns = time) under the micro rule. Use $\sigma$ and $N_{\mathrm{neurons}}$ to explore regimes, and press Play/Pause to run. When activity dies out, the demo auto-reseeds; the brief flashing border marks a reseeding event.
At the macro level, the system is a population of $N_{\mathrm{agents}}$ light-signaling agents embedded in a 2D toroidal space. Agents interact only through local visual perception of illuminated neighbors, which induces a spatial interaction structure.
Each agent emits a binary light state $l_{t,i}\in\{0,1\}$, taken to be the state of a single designated readout neuron in its reservoir (one binary unit, not a pooled or vector-valued readout). At each macro step, neighbors within vision radius $v_r$ are projected into $E$ egocentric angular sectors, producing a binary sensory vector. Those $E$ inputs are clamped to designated input neurons, distinct from the readout.
What is this? A representation of one agent and its internal reservoir: the inner graph is the reservoir, and the surrounding wedges are exteroceptive vision sectors at radius $v_r$. The designated readout (light) neuron is drawn at the center of the reservoir ring; while the agent’s light is ON, that center node is drawn in bright yellow. Agents only receive observations from other agents when those agents are ON and inside the focal agent’s vision field ($v_r$ and angular sectors). The circular pointer represents another agent’s light source: gray means inactive, yellow means ON. Click or hold in a sector to continuously stream that external light input and trigger reservoir activity. Adjust $\sigma$, $v_r$, and $E$ to see how sensing geometry and micro dynamics change the response.
Avalanche events are extracted from the population activity signal:
$$A_l(t)=\sum_{i=1}^{N_{\mathrm{agents}}} l_i(t)$$
An avalanche starts at $t_s$ when $A_l(t_s-1)=0$ and $A_l(t_s)>0$, and ends at $t_e$ when $A_l(t_e)>0$ and $A_l(t_e+1)=0$. Its duration and size are:
$$T=t_e-t_s+1,\qquad S=\sum_{t=t_s}^{t_e}A_l(t)$$
Avalanches still active when the recording stops are censored and excluded from analysis, because the finite observation window does not show whether the cascade would later return to silence.
What is this? A single-cascade extraction demo. Press Fire one cascade to seed one active neuron at $t=0$ (left); $T$ and $S$ use the same population-activity definitions as above (here: neuron count per timestep). If the cascade is still active after 260 steps it is discarded (censored by the finite window).
We start from the list of avalanche observables and analyze size $S$ and duration $T$ separately. Writing either one as $x\in\{S,T\}$, we follow the Clauset procedure [2] and fit a power law only to the tail of the empirical distribution. The reason is that the smallest events are typically the least well described by scaling: they are strongly affected by discreteness, finite-size effects, and other non-power-law structure. Therefore a lower cutoff $x_{\min}$ is introduced, and only observations with $x\ge x_{\min}$ are used in the fit.
Because both avalanche size and duration are integer-valued, the fitted model is a discrete probability mass: $P(x)$ is the probability assigned to each integer $x$ in the tail,
$$P(x)=\frac{x^{-\alpha}}{\zeta(\alpha,x_{\min})},\qquad x=x_{\min},x_{\min}+1,\ldots$$
where $\alpha$ is the scaling exponent. The denominator
$$\zeta(\alpha,x_{\min})=\sum_{n=x_{\min}}^\infty n^{-\alpha}$$
is the Hurwitz zeta function. It appears because the raw decay $x^{-\alpha}$ is not yet a probability distribution; it must be normalized over the allowed tail values $x_{\min},x_{\min}+1,\ldots$. Dividing by $\zeta(\alpha,x_{\min})$ makes the probabilities sum to $1$, so this renormalized expression is a proper discrete model.
The cutoff $x_{\min}$ is chosen from the data. For each candidate value of $x_{\min}$, we keep only the tail observations $x_i\ge x_{\min}$, estimate $\alpha$ by maximum likelihood on that restricted sample, and compare the empirical and fitted cumulative probabilities. For integer $x\ge x_{\min}$, the running total of the model masses is
$$F(x)=\sum_{t=x_{\min}}^{x}P(t)$$
so $F(x)$ is the probability (under the fitted tail model) that the next draw is at most $x$. In the same way, $F_{\mathrm{emp}}(x)$ is the fraction of tail observations in the sample that are $\le x$, and $F_{\mathrm{model}}(x)$ is the value of the expression above using the fitted $P(x)$.
Their discrepancy is measured with the Kolmogorov–Smirnov statistic
$$\mathrm{KS}(x_{\min})=\max_{x\ge x_{\min}}\bigl|F_{\mathrm{emp}}(x)-F_{\mathrm{model}}(x)\bigr|.$$
This is the largest vertical separation between the step functions $F_{\mathrm{emp}}$ and $F_{\mathrm{model}}$ over the tail. The selected cutoff is the value of $x_{\min}$ that minimizes this KS distance. In other words, we choose the threshold above which the data are best matched by a discrete power-law model. Once that optimal cutoff is found, the corresponding maximum-likelihood estimate of $\alpha$ is taken as the fitted exponent. The same procedure is carried out independently for avalanche size $S$ and avalanche duration $T$.
IMPORTANT NOTE: For speed and responsiveness, the simulation in this web version fixes $x_{\min}=2$ and estimates $\alpha_S$ and $\alpha_T$ with the same discrete maximum likelihood on the raw tail data, without re-running the $\mathrm{KS}$ cutoff search on every frame. The visualizations below show both methodologies.
Demo: Run many trials; each yields a pair $(S,T)$. The two panels use the same binned log–log data but different rules for the tail cutoff $x_{\min}$.
Fixed $x_{\min}$
Best $x_{\min}$ for $S$ and for $T$ via $\mathrm{KS}$ sweep (paper methodology)
After power-law fitting, each $(p^{\mathrm{micro}},v_r)$ configuration yields exponents $\alpha_S$ and $\alpha_T$ for avalanche size and duration, together with the Kolmogorov–Smirnov distances $\mathrm{KS}_S$ and $\mathrm{KS}_T$ at the selected cutoffs. To summarize how close the statistics are to reference critical scaling, we compare $\alpha_S$ to $\frac{3}{2}$ and $\alpha_T$ to $2$. In mean-field critical branching processes, $\frac{3}{2}$ is the canonical tail exponent for avalanche sizes and $2$ the corresponding benchmark for avalanche durations; those values align with the exponents highlighted in experimental studies of neuronal avalanches [1] and recur in broader discussions of criticality in living neural systems [3] and in controlled spiking neuromorphic networks [4]. The scalar loss used in the phase-diagram analysis is
$$\mathcal{L}(p^{\mathrm{micro}},v_r)=|\alpha_S-1.5|+|\alpha_T-2|+\mathrm{KS}_S+\mathrm{KS}_T$$
where the absolute-value terms measure mismatch to the reference exponents (written as $1.5$ and $2$ above for direct comparison with numerical estimates) and the $\mathrm{KS}$ terms retain the tail goodness-of-fit penalties from the Clauset procedure [2]. Lower $\mathcal{L}$ therefore indicates exponents nearer the branching-process benchmarks together with smaller KS discrepancies.
In the live viewer, the effective branching ratio trace and the criticality-score $\mathcal{L}$ heatmap are produced from the stored outputs of the paper simulations and follow the paper’s methodology. The on-screen power-law readouts (avalanche log–log panel and the $\alpha_S$, $\alpha_T$ reported beside it) instead use the simplified fitting procedure summarized above—fixed $x_{\min}$ and no full $\mathrm{KS}$ search every frame—so the page stays fast and responsive.
For the same reason, each interactive run is capped at 500 macro steps per adjacency realization, shorter than the full recording used in the paper pipeline. Implementation details, parameter choices, and how the archived curves and heatmap were generated are documented in the paper’s repository: https://github.com/nhbess/emergent-macro-criticality.
To explore the system dynamics, $\mathcal{L}$ heatmap, and power-law readouts together in one interface, Go to live simulation.