Your new experience awaits. Try the new design now and help us make it even better

ORIGINAL RESEARCH article

Front. Comput. Neurosci., 23 January 2026

Volume 19 - 2025 | https://doi.org/10.3389/fncom.2025.1724190

This article is part of the Research TopicExploring the Role of Criticality in Neural SystemsView all 5 articles

Measurement effects on critical scaling in neural systems

  • 1School of Physics, Georgia Institute of Technology, Atlanta, GA, United States
  • 2School of Psychology, Georgia Institute of Technology, Atlanta, GA, United States

The recently developed phenomenological renormalization group (pRG) analysis has uncovered scale-free properties in large-scale neural population recordings across recording modalities, including extracellular electrophysiology and calcium imaging. The convergence of these properties across the datasets hints at universal neural behavior. Yet, it is unknown how differences in temporal resolution and measurement details affect pRG scaling. Here, we use a network model known to produce scaling under pRG analysis as a testbed to assess how recording and analysis choices shape inferred scaling exponents. We show that scaling properties depend on the choices of temporal binning, measurement nonlinearities, and deconvolution, and that the quality of scaling for cluster covariance eigenvalues is particularly sensitive to measurement effects. Moreover, all scaling exponents shift substantially with these transformations, even when the underlying neural dynamics are identical. Together, these results show how experimental choices can change pRG scaling and provide a framework for separating scaling driven by neural dynamics from that introduced by the recording method.

1 Introduction

Criticality is a leading candidate framework for understanding brain function and dysfunction (Chialvo, 2010; Hengen and Shew, 2025), potentially accounting for what dynamical states are optimal for learning and computation (Legenstein and Maass, 2007; Kinouchi and Copelli, 2006; Shew and Plenz, 2013; Safavi et al., 2024; Cambrainha et al., 2025) and what goes awry in disease (Rocha et al., 2022; Xu et al., 2024; McGregor et al., 2024). Over the past few decades, criticality has been extensively explored in neural data, including spiking activity in small populations, LFP, and M/EEG (Beggs and Plenz, 2003; Petermann et al., 2009; Shew et al., 2015; Fontenele et al., 2019; Ponce-Alvarez et al., 2018), and much work has examined the challenges of interpreting such analyses (Touboul and Destexhe, 2017; Aitchison et al., 2016; Priesemann and Shriki, 2018; Levina et al., 2022; Morrell et al., 2024). Neural populations numbering in the thousands are now routinely recorded by modern imaging and extracellular electrophysiology (Stringer et al., 2019; Chung et al., 2019; Steinmetz et al., 2021; Demas et al., 2021; Angelaki et al., 2025), and new datasets open new approaches to the study of criticality. While these high-dimensional datasets are complex, insight from statistical physics suggests that not all details matter: the renormalization group (RG) shows us that microscopic details can be irrelevant to macroscopic behavior. In this spirit, a phenomenological renormalization group (pRG) approach was recently adapted from statistical physics to neural data (Meshulam et al., 2019, 2023). This analysis systematically coarse-grains population activity to reveal scale-invariant relationships across levels of description. The pRG analysis has now been applied across multiple recording modalities, including electrophysiology across brain areas (Morales et al., 2023), cortical states (Castro et al., 2024), and behavioral tasks (Cambrainha et al., 2025), as well as calcium imaging in diverse species (Munn et al., 2024) and fMRI in humans (Calvo et al., 2025). However, even when robust scaling is observed, the meaning of scaling exponents, and their interpretation across diverse systems and recording methods, remain uncertain. In particular, it is unknown which details of experimental design and measurement matter for pRG analysis.

The first work on RG-based approaches to the analysis of large neural populations was based on recordings of 2000+ neurons in the mouse hippocampus using two-photon calcium imaging (Meshulam et al., 2019). Briefly, their analysis proceeds by iteratively clustering neurons into groups based on activity correlations and assessing cluster activity distributions, variance, covariance matrix spectrum, and autocorrelation across scales. When scaling of these quantities is observed, the outcome of this analysis is a set of scaling exponents that reflect the spatial and temporal correlation structure of the system. Meshulam et al. (2019) showed that several statistical and dynamical properties of neural activity collapsed onto the same forms, which were invariant across two orders of magnitude in scale.

The pRG analyses have now been applied across diverse neural systems and recording modalities, and in several cases the resulting scaling exponents appear surprisingly consistent. For example, Munn et al. (2024) examined PRG scaling across five different species, using recordings from very different experimental set-ups, and found nearly identical scaling of variance. Additionally, similar values of scaling exponents were reported in mouse hippocampus using two-photon calcium imaging (Meshulam et al., 2019) and later using high-density extracellular electrophysiology (Morales et al., 2023), despite substantial differences in temporal resolution, processing pipelines, and measurement physics. Each of these factors could influence the estimated exponents. Understanding the extent to which such observational factors shape pRG scaling is essential for determining when comparisons across modalities, species, or experiments reflect underlying neural dynamics rather than differences in data acquisition or processing. Recent work has begun to explore aspects of this issue, including supplemental analyses in multiple studies examining how exponents vary with temporal coarse-graining (Morales et al., 2023; Castro et al., 2024), though these results are not yet consistent across datasets. Together, these considerations underscore the need for a systematic framework to characterize how pRG outcomes depend on the transformations introduced by measurement and processing.

Directly addressing these questions with experimental data is challenging, as it would require simultaneously recording large populations of neurons with both electrophysiology and imaging. Instead, here we use a model of spiking activity whose pRG scaling properties match those observed experimentally (Morrell et al., 2021). Our aim is to examine how observational transformations affect the inferred exponents. We instantiate a experimentally tuned forward model using GCaMP6f-like kinetics (Wei et al., 2020) and apply commonly used out-of-the-box deconvolution tools (Giovannucci et al., 2019). With these controlled conditions, we examine three classes of observational transformations:

1. Temporal resolution: how the choice of spike-train bin size shapes the apparent existence and quality of scaling.

2. Nonlinear measurement dynamics: how the fluorescence transients associated with calcium indicators alter scaling exponents.

3. Deconvolution and processing: whether “raw” pRG exponents are recoverable from standard imaging pipelines.

In the remainder of the paper, we first provide an overview of the pRG analysis and our modeling testbed. We then examine how the scaling of variance (α), free energy (β), and covariance structure (μ) are modulated by observational transformations. We observe that scaling quality is maintained for variance and free energy, but degrades for the covariance structure, and we find changes in all exponents with these transformations. These results establish that observational transformations can differentially impact the exponents recovered by pRG analyses. Taken together, this work provides a framework that can be used to distinguish scaling signatures that reflect neural dynamics from those arising from experimental and analytical constraints, and for clarifying how methodological choices impact efforts to identify criticality in neural systems.

2 Scaling and temporal binning

2.1 Overview of the pRG analysis

Meshulam et al. (2019) introduced the pRG procedure to quantify the self similarity of the system across scales. In the “real space” version of the pRG analysis, neurons are iteratively clustered based on maximum pairwise correlations, and properties of these clusters are tracked across scales (see Section 5.2). Several properties scale with the cluster size: the variance of cluster activity (σk2~Kα), the “free energy” of the cluster (F=-ln Psilence~K-β), the eigenvalues of the cluster covariance (λ(r) ~ (r/K)−μ), and the time constant of cluster activity autocorrelation (τc~Kz~). Here, we focus on the three “stationary” scaling relations (α, β, μ). The exponents characterizing these relationships depend on the correlation structure within the population, and they deviate from the trivial scaling relationship expected in any independent system (i.e., α = 1, β = 1, μ = 0). The fourth exponent, z~, is dynamical in nature and depends on the intrinsic time constants of the activity, and we discuss this in more detail in Section 4.

Cross-correlations between neurons computed from spike-counts are well-known to depend on the spike count bin size (Δ) (Kass and Ventura, 2006; Cohen and Kohn, 2011). Throughout the pRG analysis, correlations play a large role: they are used to determine coarse-grained clusters, and one of the scaling relationships observed derives from the eigenvalue spectrum of the cluster covariance matrix. Thus, at multiple points in the coarse-graining analysis, time bin size has a role. In this section, we first introduce a network model with known scaling behavior (Morrell et al., 2021) and then examine the sensitivity of the coarse-graining analysis to the choice of temporal binning.

2.2 A model testbed for scaling behavior

The model has been described in detail elsewhere (Morrell et al., 2021, 2024). Briefly, individual neurons are modeled as binary units for which the probability of firing in any time bin is determined by the activation of NF latent fields, weighted by the coupling Jim~N(0,NF-1) between the neuron i and each of the fields m. Each latent field is modeled as a mean-reverting random walk with autocorrelation time τF. The overall strength of coupling to latent fields is controlled by the parameter η (“coupling”). A global bias ϵ adjusts the average population firing rate. We simulated activity in N = 210 neurons with η = 3, ϵ = 8 with NF = 10 fields, equivalent to the parameters used in Morrell et al. (2021) that produced good scaling (Table 1; see Supplementary Figure 1).

Table 1
www.frontiersin.org

Table 1. Dynamic latent variable model parameters.

Because the model yields the probability of activity in each time step (see Section 5.1, Equation 1), setting the step size dt fixes the range of firing rates in the network. We set the simulation step size dt = 0.01 s, and so the distribution of neural firing rates ranged from 0.1 to 10 Hz. For simplicity, all latent fields have the same time constant (τF = 1 s). We defined the spike count bin to be Δ, and this ranged from 0.01 s (= dt) to 10 s. A sample of binned spike counts from 200 neurons over a 30 s period are shown for bin sizes Δ = 0.01, 0.1, 1 s in Figures 1AC. The total simulation length was T = 20, 000 s (Section 5.1; see Supplementary Figure 2 for finite-length analysis).

Figure 1
Graphs A, B, and C show neuronal spike activity over 30 seconds with time-bins of 0.01s, 0.1s, and 1s, respectively. Panel D displays probability density of correlation for random, 0.01s, and 10s time bins. Panels E-G depicts scaling for activity at 0.01s and 1.0s timbinning. E depicts variance as a function of network size K. Panel F illustrates free energy vs. K for different time intervals. Panel G shows eigenvalues against r/K.

Figure 1. Coarse-graining analysis results depend on temporal binning. (A–C) Illustration of temporal binning, for Δ = 0.01, 0.1, 1 s, on a 30-s sample of activity in 200 neurons simulated using the dynamic latent variable model (details, Section 5.1). (D) Distribution of correlations for simulated data with time bins ranging from Δ = 0.01 s (orange) to Δ = 10 s (dark blue). Gray: distribution of correlations for Gaussian random data with number of samples matched to the full simulation binned at Δ = 10 s. (E–G) pRG analysis (details, Section 5.2) results for Δ = 0.01 s and Δ = 1 s. (E) Cluster variance vs. cluster size, on log-log scale. (F) Free energy (defined as the negative natural logarithm of the probability of silence in a cluster) vs. cluster size, on log-log scale. (G) Eigenvalues vs. normalized rank, on log-log scale.

We examined the distribution of correlations as Δ was increased (Figure 1D). While coupling in the model is random, the correlations are always much larger than expected from randomized data matched to simulation length (gray curve, Figure 1D). When the Δ ≪ τF, correlations are determined by the coupling Jim to the latent field hm(t), but correlations are small because spiking is probabilistic and spike rates are low. As Δ increases, noise in the spike generation is averaged out, and correlations increase. Thus, we expect that as Δ increases, exponents will move away from the random (independent) model values: α and μ will increase, while β will decrease.

2.3 Loss of scaling quality with larger time bins

When coarse graining was performed at the simulation resolution (Figure 1A), we found clear scaling in the cluster variance (Figure 1E, orange), free energy (Figure 1, orange), and eigenvalue spectrum (Figure 1D, orange). Using the same simulation, we re-binned the spiking activity with Δ = 1 s and repeated the coarse graining analysis. Qualitatively, scaling was preserved, but scaling relationship specifics changed (blue lines, Figures 1EG). As expected, α and μ increased, while β decreased.

To examine this in greater detail, we next examined the scaling relationships and quality of scaling as a function of bin size Δ (Figure 2). Quality of scaling was assessed by R2 of the linear fit in log-log space computed over log-spaced points (see Section 5.2) and, for eigenvalue scaling, by examining the spread in slope (μ) across cluster sizes K (Figure 2D). At the simulation resolution (Δ = dt = 0.01 s), the quality of scaling of the variance with the cluster size (σ2 ~ Kα ) was very high (R2 > 0.999), with scaling exponent α = 1.323 ± 0.003. As Δ increased, variance still scaled with cluster size (R2 > 0.999), but the scaling exponent increased to α = 1.51 ± 0.01 (Figures 2A, D). Similarly, the quality of free energy scaling remained high (R2 > 0.99) as Δ increased (Figure 2D). At simulation resolution we found scaling exponent β = 0.838 ± 0.002, and at the maximum Δ = 10 s we found β = 0.714 ± 0.010 (Figure 2B).

Figure 2
Graphs A to D show how scaling changes with changing timebin size 1.
A depicts variance scaling (α) increasing with 1.
B represents free energy scaling (β) decreasing with 1.
C illustrates eigenvalue scaling (μ) for different values of K, each gradually
increasing with 1.
D shows scaling R2 comparison, with α and β remaining near ceiling with increasing \Delta while μ at various K values decreasing substantially.

Figure 2. Scaling exponents depend on temporal binning. (A) Variance scaling exponent α increases with bin size Δ. (B) Free energy scaling exponent β decreases with bin size Δ. (C) Eigenvalue scaling exponent μ increases with bin size Δ. (D) Quality of scaling relationship for variance (α, teal) and free energy (β, plum) is high regardless of bin size while eigenvalue scaling quality degrades with larger bin sizes (solid lines). Error bars are ±1s.d. across fifths of the simulation.

At the simulation timescale (minimum Δ = dt), the quality of eigenvalue scaling was high (R2 > 0.99 across all K) and values of μ changed little with K, ranging from μ = 0.66 ± 0.02 at K = 32 to μ = 0.59 ± 0.01 at K = 256. As Δ increased, the quality of eigenvalue scaling degraded, with R2 = 0.98 for K = 32 decreasing to R2 = 0.94 for K = 256, and increasing spread with cluster size K (Figure 2D). For instance, at Δ = 1 s, the scaling exponent μ ranged from μ = 1.41 ± 0.01 for K = 32 to μ = 1.13 ± 0.01 for K = 256. While the overall range across K for μ was only 0.07 for Δ = 0.01 s, it was four times larger (0.28) for Δ = 1 s (Figure 2C), suggesting a breakdown in the quality of scaling.

In summary, variance and free energy scaling is very robust: the quality of scaling is high at all bin sizes. By contrast, increasing bin size degrades the quality of scaling for eigenvalue vs. rank. In all cases, the value of the scaling exponent depends strongly on bin size.

3 Scaling in simulated imaged populations

The pRG exponents reported from imaged activity (Meshulam et al., 2019) and electrophysiologically recorded activity (Morales et al., 2023) in the hippocampus were similar, within small experimental error bars. Given the results of the previous section that exponents depend on the temporal resolution of recorded activity, this similarity is surprising, the more so because temporal resolution is only one difference between imaging and electrophysiology. Usually, neural activity is inferred by deconvolving imaged fluorescence, which reflects calcium dynamics as well as the indicator response function, both of which are heterogeneous across cells. Thus, having established that temporal coarse-graining alone can distort apparent scaling exponents, we next asked what distortions arise from the calcium-imaging process itself. In this section, we implement a forward model for imaged calcium dynamics and apply a standard deconvolution algorithm to infer spiking activity. We then compare the results of a coarse-graining analysis applied to spiking activity, simulated imaged calcium, and deconvolved activity.

3.1 Simulation: calcium dynamics, fluourescence, and deconvolution

Obtaining simultaneous large-scale electrophysiology and calcium imaging from the same neuronal population is not currently feasible. However, paired electrophysiology-imaging recordings at the single-cell level constrain forward models from spikes to fluorescence (Wei et al., 2020), enabling a controlled simulation-based approach to examine how calcium indicator dynamics and noise influence the pRG exponents. Starting from the simulated spiking activity analyzed in Section 2 (single cell in Figure 3Ai; population subset, Figure 3C), we generated fluorescence traces using the forward model of Wei et al. (2020) and sampled the resulting signals at 10 Hz (Δ = 100 ms), matching typical imaging frame rates (Figure 3Bii; population, Figure 3D). To set model parameters (τd, τr, Fmax, c1/2, k) across the population (Figure 3B, Table 2), we approximated distributions shown for GCamp6f-TG (Wei et al., 2020; see Figure 4 therein). We deliberately include this heterogeneity in model parameters to reflect heterogeneity observed across real populations of neurons. We simulated fluorescence traces with three levels of external noise, σext = 0.02, 0.2, and 0.4.

Figure 3
Graphs illustrates calcium imaging simulation. Panel A shows simulated spiking, convovled calcium signal traces and deconvolution for one neuron. Panel B displays probability density functions of various parameters used in convolution. Panels C to E compare spike counts, simulated fluorescence, and deconvolved activity over time for 20 neurons. Panel F depicts scatterplot comparing recovered and true decay times. Panel G shows cumulative fraction curves for deconvolution accuracy with varying noise. Panel H is a scatterplot comparing correlation of imaged and spike count data.

Figure 3. Forward model for calcium imaging and recovery of population activity. (A) Example single cell showing (i) simulated spike train, (ii) synthetic fluorescence generated using a GCaMP6f-like forward model, and (iii) deconvolved activity obtained with the constrained foopsi algorithm (Giovannucci et al., 2019). Spikes are simulated at 10 ms resolution and fluorescence is sampled at 10 Hz (Δ = 100 ms). (B) Distributions of calcium-model parameters across the population (τr, τd, c1/2, Fm), chosen to approximate those reported for GCaMP6f-TG (Wei et al., 2020). (C–E) Population activity for 20 example cells over a 30 s window: (C) binned spike counts, (D) simulated fluorescence traces, and (E) deconvolved activity. In (D), each row is normalized by its range for visibility. (F) Recovered decay time constants are strongly correlated with ground-truth values in the low-noise condition. (G) Deconvolution accuracy (correlation between inferred and true spike counts) declines systematically with increasing external noise. (H) Pairwise correlations from fluorescence (Ca(t)) and deconvolved activity (d-Ca(t)) closely track spike-count correlations, with smaller deviations after deconvolution.

Table 2
www.frontiersin.org

Table 2. Simulated calcium imaging parameters.

We then deconvolved the simulated fluorescence using standard methods (Giovannucci et al., 2019), producing inferred spike trains for each neuron (Figure 3Aiii; population, Figure 3E). Calcium decay time constants were inferred from simulated traces using short segments of data (10 segments of 100 s), and we confirmed that recovered time constants were close to the ground-truth values for most cells (Figure 3F), with larger errors for some cells with τd ≈ 1 s. The latent dynamics of the population model have an autocorrelation time τF = 1 s, and these errors may arise from cells that are strongly coupled to the latent dynamics.

In our our controlled, simulation-based environment, deconvolution was highly accurate, but decreased at higher noise levels. For the low-noise simulation, deconvolution ran successfully for 1,017 of 1,024 cells, and the median accuracy, quantified by the correlation between ground-truth spike counts and the deconvolved signal, was 0.86 (inter-quartile range (IQR), 0.05) (Figure 3E). Accuracy degraded with increasing noise level, with fewer activity traces successfully deconvolved (medium noise: 908 of 1024 successful; high noise: 705 of 1024 successful) and lower accuracy following successful deconvolution (medium noise: 0.72 (IQR, 0.18); high noise: 0.63 (IQR, 0.22)).

Using our simulations, we next asked how the transformations from spiking activity to fluorescence to deconvolved activity affected the cell-cell correlations that form the basis of the pRG coarse-graining (Figure 3H). Correlations were highly preserved: Spearman's ρ was 0.98 between the spike-count and simulated fluorescence correlations, and likewise 0.98 between spike-count and deconvolved-trace correlations. Correlation magnitudes were slightly increased in the simulated fluorescence (mean change +0.004 ± 0.018 relative to spike counts) and slightly decreased after deconvolution (−0.002 ± 0.019 relative to spike counts).

For fluorescence traces, we analyzed scaling relationships at the native imaging resolution (Δ = 100 ms). Unlike spike counts, where temporal coarse-graining can be performed over sub-millisecond resolution, calcium imaging already imposes a coarse temporal sampling interval, and additional binning would substantially reduce the number of independent time samples. Consistent with previous pRG analyses of calcium imaging (Meshulam et al., 2019; Munn et al., 2024), we therefore analyzed the fluorescence traces at a single value of Δ.

3.2 Scaling of variance and free-energy in simulated imaged activity

As in the earlier analyses, activity was coarse-grained into clusters of size K, now for each of seven conditions: spike counts (Δ = 0.1 s), simulated fluorescence traces (dF/F, for low, medium, and high noise) and deconvolved fluorescence traces (dCa(t)) at the same noise levels. Variance scaling fit quality remained high (R2 > 0.99) across all levels of noise, for both simulated fluorescence and the deconvolved signals (Figures 4A, B). Increasing noise levels decreases α, following expectations: fully uncorrelated data would show α = 1. We found the scaling exponent α was substantially lower for the simulated fluorescence (α = 1.334 ± 0.004 for the low-noise case) than for spike counts (α = 1.46 ± 0.01) or for deconvolved spike counts (α = 1.43 ± 0.01). One reason for this may be that the simulated fluorescence reflects the timing of the activity within the spike count window Δ, which adds variability to the signal that is not present in the summed spike counts in a window of the same size, thereby decreasing α. Deconvolution, which infers the spike count, appears to largely remove this variability, as α is similar between spike count and deconvolved data. This suggests that variance scaling can be recovered across multiple recording modalities.

Figure 4
Data visualization consisting of four panels: A) A line graph showing variance versus K with different colored lines indicating spiking data and transformed spiking using different convolution variables and resulting deconvolutions. B) A scatter plot labeled alpha on the y-axis with different data points for the different transformations. C) A line graph labeled -ln Psilence versus K showing different deconvolution transformations along with original spiking. D) A bar graph labeled beta, with data points for spikes and different deconvolution transformations. A legend provides color coding for spikes, dF/F, and d-Ca with varying noise levels.

Figure 4. Variance and free-energy scaling exponents derived from spikes, simulated fluorescence, and deconvolved calcium signals. (A) Variance scales with cluster size as Kα for spike counts (Δ = 100 ms), simulated fluorescence traces (dF/F) at three noise levels, and deconvolved fluorescence traces (dCa(t)). All conditions exhibit a clear power-law relationship. (B) The exponent α decreases with increasing external noise in the fluorescence simulation. Deconvolution partially restores α toward the spike-based value. (C) Free-energy, quantified by −ln Psilence, scales with cluster size as Kβ for spike counts and deconvolved fluorescence traces. (D) The exponent β is higher in the deconvolved signal and increases with added noise. Because the simulated fluorescence traces include baseline fluctuations and lack a true zero-activity state, β cannot be defined reliably for raw dF/F; we therefore report β only for spike-count and deconvolved data. Error bars denote ±1 s.d. across fifths of the data.

We also computed the free-energy scaling exponent β (Figures 4C, D). Because the simulated fluorescence includes baseline noise and lacks a true zero-activity state, β cannot be reliably defined for the raw dF/F signals, and we therefore report β only for the spike-count and deconvolved data. As for α, the fit quality remained high for all conditions (R2 > 0.99). Increasing noise levels increases β: noise contributes spurious non-zero activity, which causes the population silence probability to decrease faster (i.e., higher β) with increasing cluster size than it would otherwise. These results indicate that both α and β are less sensitive to transformations: the quality of scaling fits remains high, while the exponents change predictably with increasing levels of noise.

3.3 Scaling of covariance eigenvalues in simulated imaged activity

Next, we examined how eigenvalue scaling transformed, from the covariance obtained from spike counts (Δ = 0.1 s, Figure 5A), imaged fluorescence (Figures 5BD, teal), and deconvolved fluorescence (Figures 5BD, red, ‘d-Ca') at low (B), medium (C), and high (D) noise levels. In each case, the plot of cluster eigenvalues against scaled rank collapsed onto a curve across cluster sizes K = 32, 64, 128, 256. Figure 5E shows this curve for the K = 256 cluster eigenvalues across conditions. We quantified the eigenvalue scaling by extracting μ for each condition at each cluster size. As observed previously, μ varies somewhat with K, even for spike counts.

Figure 5
Seven graphs labeled A to G display eigenvalues and other statistical measures for different noise levels analyzed with various methods. Graph A shows eigenvalues for spikes, graph B with low noise, C with medium noise, and D with high noise for methods dF/F and d-Ca. Graph E compares eigenvalues at K = 256. Graph F plots mu values across changing K values, while graph G shows R² for mu fits. Legends distinguish dF/F and d-Ca using different noise levels and colors. Axes show logarithmic scales for eigenvalues and varying K values.

Figure 5. Comparison of eigenvalue scaling exponents derived from spike counts, calcium traces, and deconvolved data (Δ = 100 ms). (A–D) Eigenvalues vs. rank, for clusters of size 32, 64, 128 and 256 for (A) the spike counts, and the calcium traces (teal) and deconvolved data (red-orange) for (B) low noise, (C) medium noise, and (D) high noise conditions (see Section 5.3). (E) Eigenvalues vs. rank, for clusters of size K = 256 across different experimental and noise conditions. Colors as in (A–D). (F) Eigenvalue scaling exponent μ vs. K for each signal. (G) R2 of fit to power-law decay (see Section 5.2) for each signal at each cluster size. Error bars on α and μ are the standard deviation across fifths of the data.

At very low noise levels (darkest red/teal lines, Figure 5F), the calcium dynamics shift μ from the spike-count value significantly: at all values of K, there is a positive shift in μ of 0.15 to 0.2 from spiking to simulated fluorescence. Deconvolution decreases this shift slightly, but does not recover the scaling exponents obtained from spike counts directly.

Deconvolution had a larger impact on exponents under noisier conditions (medium, high noise cases; light red/teal lines), and as expected, increasing noise in the simulated fluorescence decreased μ. The spectra plotted in Figures 5CE suggest that the reason for this decrease in μ is that the white noise adds an offset across all eigenvalues, flattening the eigenvalue vs. rank curve. The effect of deconvolution on noisy calcium traces is to increase μ, by 0.02 to 0.14 for medium noise (σ = 0.2) and 0.11 to 0.28 for high noise (σ = 0.4).

Finally, we examined fit quality across conditions. The quality of fit assessed by R2 was comparable to that found in experimental observations (> 0.95) (Figure 5G). However, the highest quality fits were observed in noisy simulated calcium traces (teal curves, Figure 5G), largely due to the flattening of the curve in the tail. Fit quality was lower for deconvolved signals (red curves, Figure 5G) relative to simulated calcium traces, but still was higher for higher noise than lower noise. Thus, it is possible for a white-noise background to inflate fit quality for the eigenvalue scaling observed under pRG analysis. Moreover, the scaling exponent is highly sensitive to noise level and method of recording.

4 Discussion

The phenomenological renormalization group (pRG) offers a novel approach to probing criticality in large-scale neural recordings (Meshulam et al., 2019), emphasizing how population statistics transform across coarse-grained scales rather than relying on avalanche detection or reduction to a one-dimensional activity signal, which can obscure signatures of criticality (Dahmen et al., 2019; Fontenele et al., 2025). Given this promise, which is also reflected in its rapid adoption across electrophysiology, calcium imaging, and fMRI datasets (Morales et al., 2023; Castro et al., 2024; Munn et al., 2024; Calvo et al., 2025), it becomes increasingly important to understand the methodological limits of pRG analyses. In particular, comparisons across modalities require understanding how temporal sampling, nonlinear measurement dynamics, and processing choices shape the scaling relationships used to infer criticality. In this work, we examined how temporal coarse-graining and the nonlinear transformations inherent in calcium imaging and deconvolution impact the scaling relationships extracted via pRG, and we found that both factors can substantially alter the apparent scaling exponents.

4.1 Trends for each exponent

For overall context, trivial scaling (α = 1, β = 1, μ = 0) is observed in a system of independent variables, and if scaling persists with increasing levels of correlational structure, these exponents shift systematically: α and μ increase with correlation and β decreases. The trends we observed broadly reflect this intuition: more noise drove the system closer to randomness while averaging over this noise had the opposite effect.

4.1.1 Variance scaling (α)

None of the forward-modeling steps substantially degraded the quality of the variance-cluster size relationship, and in our controlled setting, the exponent α changed in predictable ways with added noise or temporal filtering. For example, α decreased when external noise is increased and increased when noise is reduced, such as when we used larger time bins or after deconvolution. The trend with bin size was consistent with recent supplemental analyses (Morales et al., 2023, Figure S6; Castro et al., 2024, Figure S1). The most substantial departure arose when moving from spike counts to raw fluorescence: the simulated fluorescence introduced additional temporal variability tied to precise spike timing, which is not present when spikes are directly summed within a bin. This additional variability lowered the apparent α. Deconvolution appeared to remove much of this variability, with scaling exponents approaching those obtained from spike counts directly.

4.1.2 Free energy scaling (β)

Similar to variance scaling, the quality of free energy scaling remained high under the transformations examined here. When noise was reduced via larger time bins, β decreased, which is also in agreement with recent results (Morales et al., 2023, Figure S6; Castro et al., 2024, Figure S1). We did not calculate β for simulated fluorescence, because free energy calculations require defining the zero-activity state (silence), which would have required an additional thresholding step on the simulated fluorescence, which is at any rate accomplished in the deconvolution step. The addition of high levels of noise in the simulated imaging pipeline drove β toward 1 (trivial scaling, independent neurons), as deconvolution had limited success at high noise levels.

4.1.3 Eigenvalue scaling (μ)

Eigenvalue scaling was more sensitive than variance or free energy scaling to temporal binning, measurement nonlinearities, and deconvolution. This sensitivity appeared in two forms: the inferred exponent (μ) shifted systematically with noise and the quality of scaling fits changed with these transformations. Similar to variance scaling trends, μ increased when noise levels decreased. This was seen in the increase of μ with bin size Δ and in the decrease in μ with increasing external noise, from μ ≈ 1.4 under low noise to μ ≈ 0.8 under high noise. The dependence of μ on bin size has been examined previously in supplementary analyses supporting other studies (Morales et al., 2023; Castro et al., 2024), but the trends reported therein were discrepant with each other. Whether this discrepancy arises from differences in underlying dynamics [Castro et al. (2024) recorded only from the visual cortex while Morales et al. (2023) compiled recordings from across many brain areas], or a difference in processing is not clear, but our results were consistent with the supplementary analysis found by Morales et al. (2023) (see Figure S6).

One counterintuitive finding was that apparent scaling quality (as measured by R2 or collapse across cluster sizes) could improve when noise was added to the simulated fluorescence. White noise introduces a flat component across all covariance modes, masking deviations from a power-law tail and yielding apparently better scaling.

4.1.4 Future work: dynamical scaling (z~)

The pRG framework can also be applied to temporal structure by examining the scaling of autocorrelation times τc of cluster activity, yielding a dynamic exponent z~ (Meshulam et al., 2019). In the context of the dynamic latent variable model, dynamic scaling behavior depends on the intrinsic time constants of the underlying latent variables (see, e.g., Figure S7 ( C ) of Morrell et al. (2021)). Depending on the choice of time constants in the model, the quality of dynamical scaling is often substantially lower than that of the static exponents α, β, and μ. In the present work, we introduced additional temporal structure through temporal binning, calcium indicator dynamics, and deconvolution. Future work incorporating z~ will require exploring the space of intrinsic and observational time constants.

In summary, while α, β and μ ideally reflect intrinsic properties of the underlying neural dynamics, they differ in their susceptibility to measurement effects. Under identical conditions, μ, which captures structure across the full covariance eigenvalue spectrum, was more sensitive to observation noise and nonlinearities, whereas α and β were easier to measure. Even when the quality of scaling remained high, the size of the shifts we observed were on the order of the variations in exponents seen across brain areas (Morales et al., 2023). Overall, our results highlight how measurement effects can drive large differences in apparent scaling exponents.

4.2 Critical dynamics and the dynamic latent variable model

Our finding that α increases with time bin size (Δ) aligns with recent modeling work examining how temporal coarse-graining shifts α in a critically tuned random recurrent network model (Calvo et al., 2025). Of note, the exponent ranges for α (1.2 to 1.5) in a network with symmetric connections are similar to those found here. The model we use is, on the surface, quite different from a critically tuned dynamical network. However, under some circumstances, they may be equivalent, much in the same way that a restricted Boltzmann machine, with hidden variables controlling observed patterns, can approximate a Hopfield network, in which attractors arise from interactions (Smart and Zilman, 2021). In this view, our statistical network model, where population patterns are assumed to arise from coupling to a few latent variables, could be mapped onto a network with symmetric interactions, where activity samples population states close to attractors. Both frameworks produce low-dimensional manifolds of population activity governed by a small number of slow variables; the latent-variable model enforces this statistically, whereas the recurrent model realizes it dynamically. Establishing the conditions under which the statistical and dynamical models of neural activity dynamics leading to robust pRG scaling are equivalent is left to future work.

4.3 Connections to experimental work

Our simulations represent a best-case scenario for inference. While we ran the deconvolution procedure blind to underlying ground-truth parameters, the forward model used to generate fluorescence matches the model class assumed by the deconvolution algorithm, and this is unlikely to be true for real imaged population activity. Under these ideal conditions, pairwise structure is remarkably well preserved: the Spearman correlation between cell-cell correlation matrices derived from ground-truth spikes, simulated fluorescence, and deconvolved activity is high (ρ≈0.97). Even so, eigenvalue exponents and their fit qualities shifted substantially. In real data, where indicator kinetics are heterogeneous and partially misspecified, baselines drift, shot noise and motion artifacts are present, and deconvolution priors are only approximate, we expect larger departures than those reported here.

Our model-based approach comes with two main limitations. While our fluorescence simulations were guided by paired electrophysiology-imaging data (Wei et al., 2020), scaling the model to populations of thousands of neurons required additional assumptions. We focused on a single calcium indicator (GCaMP6f-TG), but extending the forward-model approach to other indicators and modalities could reveal the conditions under which scaling relationships are most reliably measured. More broadly, similar forward-modeling approaches could be used for other high-dimensional recording modalities for which pRG analysis may provide useful insight, such as fMRI. A second limitation of our model is that it could not be used to examine how differences in the geometry of sampling across modalities may influence coarse-graining. Imaging can more easily access volumes of cells, whereas multielectrode probes typically sample neurons from different cortical layers and substructures, with potentially distinct correlation statistics. Determining how geometry affects pRG analyses will require empirical investigation, for instance by sub-sampling volumetric imaging data to clarify how layer-specific or topographic organization affects apparent scaling.

The consequences of the many choices in how neural data is recorded and processed are challenging to address solely in experimental data, in which there will always be constraints in the population size and recording duration. As a testbed to examine the impact of these choices, the dynamic latent variable model is an excellent tool: it replicates pRG scaling and exponents and can be scaled up to arbitrary population size and recording duration (Morrell et al., 2021), and can be combined, as done here, with indicator-specific forward models (Wei et al., 2020) and deconvolution methods. Ultimately, an understanding of how the measurement process shapes observed scaling could lead to the development of methods to calibrate or subtract these effects. Our results indicated that some, but not all, signatures are robust to the method in which activity is measured. Thus, the present work provides both a caution and a path forward: scaling analyses remain powerful, provided the effect of transformations between neural activity and observed data are understood.

5 Methods

5.1 Dynamic latent variable model for population spiking activity.

We use the dynamic latent variable model from Morrell et al. (2021). This is a model in which each neuron is coupled to independent latent variables that determine the time evolution of the neuron activity. Each latent variable is an Orstein-Uhlenbeck random walk. The interaction between each neuron and each latent variables is drawn from a normal distribution with mean 0 and variance NF-1 as shown in Supplementary Figure 1. Aside from shared drive arising from the common set of latent variables, neurons are independent. Activity si in neuron i is determined by

H=ηi,mN,NFJimhm(t)si+ϵisi    (1)
P(si|{hm(t)})=1Z(si,{hm(t)})e-H.    (2)

In Equation 1, N is the total number of neurons, NF is the number of latent fields, hm(t) represents a given latent field, S is the state of the neuron, and η and ϵ are parameters that control the coupling strength and bias. In Equation 2, H sets the probability of activity of each neuron, given the latent variables.

The model was simulated for N = 1024 neurons for T = 2·104 s (2·106 time steps, or 5.5 h). At the largest value of Δ = 10 s, this translates to T/Δ = 2000 samples in time. For pRG analysis, scaling relationships were fit on fifths of the data, and error bars are estimated as the standard deviation across fifths of the data.

5.2 Phenomenological renormalization group (pRG) analysis

We applied the real-space phenomenological renormalization group (pRG) procedure introduced in Meshulam et al. (2019) and subsequently used in multiple large-scale neural datasets (Morales et al., 2023; Munn et al., 2024; Morrell et al., 2021; Castro et al., 2024). The goal of pRG is to examine how statistical structure changes under systematic coarse-graining, thereby revealing scale-invariant relationships in population activity. At each step, units that share the strongest correlations are merged, producing a hierarchy of coarse-grained descriptions that preserves the dominant shared fluctuations in the population. We refer the reader to the original work for a complete description, but walk through the procedure briefly here.

On iteration k = 1, the pair of neurons with the highest correlation coefficient are identified and merged into a cluster. The activity of a cluster is defined as the sum of its constituent neurons' activities, normalized by the mean activity in the non-silent state so that activity remains comparable across cluster sizes. This procedure is repeated until all neurons are assigned to clusters of size K = 2k = 2. On iteration k+1, correlations are recomputed among the clusters of size 2k, and the merging procedure is repeated to form clusters of size 2k+1. Iterating this process yields coarse-grained variables spanning scales from single neurons to large clusters.

Below, we describe the three scaling relationships analyzed in this work.

5.2.1 Variance scaling and the exponent α

For each cluster size K, we computed the variance σK2 of the coarse-grained activity. The variance-scaling exponent α was extracted as the slope of the best-fit line of logσK2 vs. log K. This exponent captures how shared variability grows with cluster size, with α = 1 corresponding to uncorrelated activity and α = 2 corresponding to perfectly correlated activity among all neurons in the cluster.

5.2.2 Free-energy scaling and the exponent β

Following Morrell et al. (2021), we computed the free-energy scaling exponent β from the growth of −ln Psilence with cluster size K, where Psilence is the probability that all neurons in a cluster are inactive in a time bin. For spike-count data, inactivity is naturally defined as a zero spike count. For deconvolved calcium traces, inactivity was defined as a bin with no inferred spikes. Because the fluorescence traces do not possess a true zero-activity state, β cannot be reliably defined for the raw dF/F signals, and we therefore report β only for spike-count and deconvolved data.

5.2.3 Covariance-spectrum scaling and the exponent μ

To characterize how the structure of correlations evolves under coarse-graining, we computed the covariance matrix for each cluster of size K, calculated its eigenvalues, and averaged the spectra across clusters. Following previous work (Meshulam et al., 2019; Morales et al., 2023; Castro et al., 2024), we plotted the eigenvalues λ(r) against their scaled rank r/K and estimated the exponent μ as the slope of the best-fit line of logλ(r) vs. log(r/K) over the range r = 1 to r = K/2, sampled logarithmically (e.g., at r = 1, 2, 4, ...K/2). This exponent quantifies how rapidly covariance modes decay at each coarse-grained scale.

We used two indicators of the quality of scaling. First, we computed R2 over logarithmically spaced ranks up to half of the maximum cluster size (1, 2, 4, 8, ...K/2). For reference, values of R2 reported for experimental data are typically larger than 0.95 (see Table S2 of Morales et al., 2023). Second, we extracted slopes μ from fits for K = 32, 64, 128, and 256; smaller differences across K indicated higher quality of scaling. Uncertainties in exponents are reported as the standard deviation across fifths of the data.

5.3 Simulated imaging and deconvolution

5.3.1 Simulating fluorescence

We generated synthetic calcium fluorescence signals by applying the forward model of Wei et al. (2020) to the simulated spike trains. Spikes at times {tk} were first converted into an intracellular calcium concentration trace c(t) using a double-exponential kernel,

c(t)=t>tkexp(-t-tkτd)[1-exp(-t-tkτr)],    (3)

where τd and τr are decay and rise constants. Internal noise was omitted for simplicity. The calcium concentration was then transformed into a synthetic fluorescence signal through a sigmoidal nonlinearity,

ΔF/Fsynth(t)=Fmax1+exp[-k(c(t)-c1/2)]+ne(t)    (4)

where Fmax is the maximum fluorescence change, k controls the steepness of the nonlinearity, c1/2 is the half-activation concentration, and ne(t)~N(0,σext2) is additive external noise. We simulated three noise levels, σext = 0.02, 0.2, 0.4.

Model parameters were drawn from distributions chosen to approximate those reported for GCaMP6f-TG in Wei et al. (2020) (their Figure 4B and Supplementary Figure 3A). Fmax was drawn from a normal distribution with mean 8 and standard deviation 2. To avoid saturation at high firing rates, the rise and decay constants were made weakly dependent on each cell's average firing rate ri: τr was drawn from a lognormal distribution in which logτr had mean log(0.05)+0.2log(ri) (s.d. 0.3), and τd was set to τr+c0/(τ0+ri), with c0 = 0.5 and τ0 = 0.4. For each neuron, c1/2 was set to half of the maximum calcium concentration. Because c1/2 and k are tightly coupled, we used the empirical relationship

c1/2=327k+1,    (5)

which was obtained from the parameters extracted from Wei et al. (2020). This relationship should be viewed as a heuristic.

Fluorescence signals were then sampled at 10 Hz (Δ = 100 ms), matching a typical imaging frame rate. In all pRG analyses below, fluorescence traces were evaluated at this native temporal resolution.

5.3.2 Deconvolution of simulated fluorescence

Fluorescence traces were deconvolved using the constrained foopsi algorithm (Giovannucci et al., 2019; Pnevmatikakis et al., 2016) with the “lars” optimization method, as implemented in the software package of Wei (2020). For each neuron, the baseline fluorescence was initialized as the median of the full fluorescence trace, and baseline drift was not modeled. Calcium decay time constants were estimated by fitting AR(2) models to ten nonoverlapping 100 s segments of the fluorescence trace, and the slow pole of the AR model was taken as the estimate of τd. In cases where no stable AR(2) model was found, we increased model order up to 5; if no stable solution was obtained, deconvolution for that neuron was excluded. This occurred most frequently at higher levels of noise, and accounts for the missing fraction of cells in Figure 3G.

Estimated decay constants were used to initialize the full deconvolution, which was run in batches (batch size 2000) over the full 200,000-frame dataset. The output was an inferred spike train for each neuron, which we used for subsequent pRG analysis.

5.3.3 Correlation analysis

Because the pRG procedure relies on the ordering of pairwise correlations to form clusters, we computed correlation matrices for spike counts, simulated fluorescence, and deconvolved traces at the same temporal resolution (Δ = 100 ms). For each noise condition, we computed Pearson correlations across time for all neuron pairs, as well as the Spearman correlation between the ranking of spike-count correlations and the rankings obtained for fluorescence and deconvolved signals. These analyses verified that pairwise correlations were highly preserved across observational transformations, providing a consistent substrate for pRG coarse-graining.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary material, further inquiries can be directed to the corresponding author.

Author contributions

ML: Investigation, Software, Writing – original draft, Writing – review & editing, Methodology, Visualization, Validation, Conceptualization. BL: Software, Validation, Methodology, Conceptualization, Investigation, Writing – review & editing. KA: Methodology, Validation, Writing – review & editing, Conceptualization. EA: Validation, Writing – original draft, Conceptualization, Writing – review & editing, Investigation, Software, Methodology. AS: Software, Investigation, Funding acquisition, Writing – original draft, Validation, Conceptualization, Resources, Writing – review & editing, Supervision, Project administration, Visualization, Methodology.

Funding

The author(s) declared that financial support was received for this work and/or its publication. AS, BL, KA, EA, and ML were supported by the National Institute of Mental Health (NIH BRAIN Grant RF1MH130413 to AS).

Conflict of interest

The author(s) declared that this work was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Generative AI statement

The author(s) declared that generative AI was used in the creation of this manuscript. Generative AI was used to solve technical formatting issues in LaTeX, to optimize memory handling and parallel processing in simulation software (MatLab), and as a tool to aid in the revision of text. The author(s) take full responsibility for the content of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fncom.2025.1724190/full#supplementary-material

References

Aitchison, L., Corradi, N., and Latham, P. E. (2016). Zipf's law arises naturally when there are underlying, unobserved variables. PLoS Comput. Biol. 12, 1–32. doi: 10.1371/journal.pcbi.1005110

PubMed Abstract | Crossref Full Text | Google Scholar

Angelaki, D., Benson, B., Benson, J., Birman, D., Bonacchi, N., Bougrova, K., et al. (2025). A brain-wide map of neural activity during complex behaviour. Nature 645, 177–191.

PubMed Abstract | Google Scholar

Beggs, J. M., and Plenz, D. (2003). Neuronal avalanches in neocortical circuits. J. Neurosci. 23, 11167–11177. doi: 10.1523/JNEUROSCI.23-35-11167.2003

PubMed Abstract | Crossref Full Text | Google Scholar

Calvo, R., Martorell, C., Roig, A., and Mu noz, M. A. (2025). Robust scaling in human brain dynamics despite latent variables and limited sampling distortions. arXiv [preprint] arXiv:2506.03640. doi: 10.48550/arXiv.2506.03640

Crossref Full Text | Google Scholar

Cambrainha, G. G., Castro, D. M., de Vasconcelos, N. A. P., Carelli, P. V., and Copelli, M. (2025). Criticality at work: Scaling in the mouse cortex enhances performance. PRX Life 3:033026. doi: 10.1103/w49n-2vz8

Crossref Full Text | Google Scholar

Castro, D. M., Feliciano, T., de Vasconcelos, N. A. P., Soares-Cunha, C., Coimbra, B., Rodrigues, A. J., et al. (2024). In and out of criticality? State-dependent scaling in the rat visual cortex. PRX Life 2:023008. doi: 10.1103/PRXLife.2.023008

Crossref Full Text | Google Scholar

Chialvo, D. R. (2010). Emergent complex neural dynamics. Nat. Phys. 6, 744–750. doi: 10.1038/nphys1803

Crossref Full Text | Google Scholar

Chung, J. E., Joo, H. R., Fan, J. L., Liu, D. F., Barnett, A. H., Chen, S., et al. (2019). High-density, long-lasting, and multi-region electrophysiological recordings using polymer electrode arrays. Neuron, 101, 21–31.e5. doi: 10.1016/j.neuron.2018.11.002

PubMed Abstract | Crossref Full Text | Google Scholar

Cohen, M. R., and Kohn, A. (2011). Measuring and interpreting neuronal correlations. Nat. Neurosci. 14, 811–819. doi: 10.1038/nn.2842

PubMed Abstract | Crossref Full Text | Google Scholar

Dahmen, D., Grün, S., Diesmann, M., and Helias, M. (2019). Second type of criticality in the brain uncovers rich multiple-neuron dynamics. Proc. Nat. Acad. Sci. 116, 13051–13060. doi: 10.1073/pnas.1818972116

PubMed Abstract | Crossref Full Text | Google Scholar

Demas, J., Manley, J., Tejera, F., Barber, K., Kim, H., Martínez Traub, F., et al. (2021). High-speed, cortex-wide volumetric recording of neuroactivity at cellular resolution using light beads microscopy. Nat. Methods 18, 1103–1111. doi: 10.1038/s41592-021-01239-8

PubMed Abstract | Crossref Full Text | Google Scholar

Fontenele, A. J., de Vasconcelos, N. A. P., Feliciano, T., Aguiar, L. A. A., Soares-Cunha, C., Coimbra, B., et al. (2019). Criticality between cortical states. Phys. Rev. Lett. 122:208101. doi: 10.1103/PhysRevLett.122.208101

PubMed Abstract | Crossref Full Text | Google Scholar

Fontenele, A. J., Sooter, J. S., Ziarati, E., Barreiro, A. K., Ly, C., and Shew, W. L. (2025). Is critical brain dynamics more prevalent than previously thought? bioRxiv. doi: 10.1101/2025.09.02.673722

Crossref Full Text | Google Scholar

Giovannucci, A., Friedrich, J., Gunn, P., Kalfon, J., Brown, B. L., Koay, S. A., et al. (2019). Caiman: an open source tool for scalable calcium imaging data analysis. Elife 8:e38173. doi: 10.7554/eLife.38173

Crossref Full Text | Google Scholar

Hengen, K. B., and Shew, W. L. (2025). Is criticality a unified setpoint of brain function? Neuron 113, 2582–2598.e2. doi: 10.1016/j.neuron.2025.05.020

PubMed Abstract | Crossref Full Text | Google Scholar

Kass, R. E., and Ventura, V. (2006). Spike count correlation increases with length of time interval in the presence of trial-to-trial variation. Neural Comput. 18, 2583–2591. doi: 10.1162/neco.2006.18.11.2583

PubMed Abstract | Crossref Full Text | Google Scholar

Kinouchi, O., and Copelli, M. (2006). Optimal dynamical range of excitable networks at criticality. Nat. Phys. 2, 348–351. doi: 10.1038/nphys289

Crossref Full Text | Google Scholar

Legenstein, R., and Maass, W. (2007). Edge of chaos and prediction of computational performance for neural circuit models. Neural Netw. 20, 323–334. doi: 10.1016/j.neunet.2007.04.017

PubMed Abstract | Crossref Full Text | Google Scholar

Levina, A., Priesemann, V., and Zierenberg, J. (2022). Tackling the subsampling problem to infer collective properties from limited data. Nat. Rev. Phys. 4, 770–784. doi: 10.1038/s42254-022-00532-5

Crossref Full Text | Google Scholar

McGregor, J. N., Farris, C. A., Ensley, S., Schneider, A., Fosque, L. J., Wang, C., et al. (2024). Failure in a population: tauopathy disrupts homeostatic set-points in emergent dynamics despite stability in the constituent neurons. Neuron 112, 3567–3584.e5. doi: 10.1016/j.neuron.2024.08.006

PubMed Abstract | Crossref Full Text | Google Scholar

Meshulam, L., Gauthier, J. L., Brody, C. D., Tank, D. W., and Bialek, W. (2019). Coarse graining, fixed points, and scaling in a large population of neurons. Phys. Rev. Lett., 123:178103. doi: 10.1103/PhysRevLett.123.178103

PubMed Abstract | Crossref Full Text | Google Scholar

Meshulam, L., Gauthier, J. L., Brody, C. D., Tank, D. W., and Bialek, W. (2023). Successes and failures of simple statistical physics models for a network of real neurons. arXiv [preprint] arXiv:2112.14735. doi: 10.48550/arXiv.2112.14735

Crossref Full Text | Google Scholar

Morales, G. B., di Santo, S., and Mu noz, M. A. (2023). Quasiuniversal scaling in mouse-brain neuronal activity stems from edge-of-instability critical dynamics. Proc. Nat. Acad. Sci. 120:e2208998120. doi: 10.1073/pnas.2208998120

PubMed Abstract | Crossref Full Text | Google Scholar

Morrell, M. C., Nemenman, I., and Sederberg, A. J. (2024). Neural criticality from effective latent variables. eLife 12:RP89337. doi: 10.7554/eLife.89337

PubMed Abstract | Crossref Full Text | Google Scholar

Morrell, M. C., Sederberg, A. J., and Nemenman, I. (2021). Latent dynamical variables produce signatures of spatiotemporal criticality in large biological systems. Phys. Rev. Lett., 126:118302. doi: 10.1103/PhysRevLett.126.118302

PubMed Abstract | Crossref Full Text | Google Scholar

Munn, B. R., Müller, E. J., Favre-Bulle, I., Scott, E., Lizier, J. T., Breakspear, M., et al. (2024). Multiscale organization of neuronal activity unifies scale-dependent theories of brain function. Cell 187, 7303–7313.e15. doi: 10.1016/j.cell.2024.10.004

PubMed Abstract | Crossref Full Text | Google Scholar

Petermann, T., Thiagarajan, T. C., Lebedev, M. A., Nicolelis, M. A. L., Chialvo, D. R., and Plenz, D. (2009). Spontaneous cortical activity in awake monkeys composed of neuronal avalanches. Proc. Nat. Acad. Sci. 106, 15921–15926. doi: 10.1073/pnas.0904089106

PubMed Abstract | Crossref Full Text | Google Scholar

Pnevmatikakis, E. A., Soudry, D., Gao, Y., Machado, T. A., Merel, J., Pfau, D., et al. (2016). Simultaneous denoising, deconvolution, and demixing of calcium imaging data. Neuron 89, 285–299. doi: 10.1016/j.neuron.2015.11.037

PubMed Abstract | Crossref Full Text | Google Scholar

Ponce-Alvarez, A., Jouary, A., Privat, M., Deco, G., and Sumbre, G. (2018). Whole-brain neuronal activity displays crackling noise dynamics. Neuron 100, 1446–1459.e6. doi: 10.1016/j.neuron.2018.10.045

PubMed Abstract | Crossref Full Text | Google Scholar

Priesemann, V., and Shriki, O. (2018). Can a time varying external drive give rise to apparent criticality in neural systems? PLoS Comput. Biol. 14, 1–29. doi: 10.1371/journal.pcbi.1006081

PubMed Abstract | Crossref Full Text | Google Scholar

Rocha, R. P., Koçillari, L., Suweis, S., De Filippo De Grazia, M., Thiebaut de Schotten, M., Zorzi, M., et al. (2022). Recovery of neural dynamics criticality in personalized whole-brain models of stroke. Nat. Commun. 13:3683. doi: 10.1038/s41467-022-30892-6

PubMed Abstract | Crossref Full Text | Google Scholar

Safavi, S., Chalk, M., Logothetis, N. K., and Levina, A. (2024). Signatures of criticality in efficient coding networks. Proc. Nat. Acad. Sci. 121:e2302730121. doi: 10.1073/pnas.2302730121

PubMed Abstract | Crossref Full Text | Google Scholar

Shew, W. L., Clawson, W. P., Pobst, J., Karimipanah, Y., Wright, N. C., and Wessel, R. (2015). Adaptation to sensory input tunes visual cortex to criticality. Nat. Phys. 11, 659–663. doi: 10.1038/nphys3370

Crossref Full Text | Google Scholar

Shew, W. L., and Plenz, D. (2013). The functional benefits of criticality in the cortex. The Neuroscientist 19, 88–100. doi: 10.1177/1073858412445487

PubMed Abstract | Crossref Full Text | Google Scholar

Smart, M., and Zilman, A. (2021). “On the mapping between Hopfield networks and restricted Boltzmann machines,” in International Conference on Learning Representations (Amherst, MA: OpenReview.net). Available online at: https://openreview.net/forum?id=RGJbergVIoO (Accessed October 10, 2025).

Google Scholar

Steinmetz, N. A., Aydin, C., Lebedeva, A., Okun, M., Pachitariu, M., Bauza, M., et al. (2021). Neuropixels 2.0: a miniaturized high-density probe for stable, long-term brain recordings. Science 372:eabf4588. doi: 10.1101/2020.10.27.358291

PubMed Abstract | Crossref Full Text | Google Scholar

Stringer, C., Pachitariu, M., Steinmetz, N., Reddy, C. B., Carandini, M., and Harris, K. D. (2019). Spontaneous behaviors drive multidimensional, brainwide activity. Science 364:eaav7893. doi: 10.1126/science.aav7893

PubMed Abstract | Crossref Full Text | Google Scholar

Touboul, J., and Destexhe, A. (2017). Power-law statistics and universal scaling in the absence of criticality. Phys. Rev. E. 95:012413. doi: 10.1103/PhysRevE.95.012413

PubMed Abstract | Crossref Full Text | Google Scholar

Wei, Z. (2020). Forward and deconvolution models for spike and fluorescence transform. Zenodo, 1.0.1. Software.

Google Scholar

Wei, Z., Lin, B. J., Chen, T., Daie, K., Svoboda, K., and Druckmann, S. (2020). A comparison of neuronal population dynamics measured with calcium imaging and electrophysiology. PLoS Comput. Biol. 16:e1008198. doi: 10.1371/journal.pcbi.1008198

PubMed Abstract | Crossref Full Text | Google Scholar

Xu, Y., Schneider, A., Wessel, R., and Hengen, K. (2024). Sleep restores an optimal computational regime in cortical networks. Nat. Neurosci. 27, 328–338. doi: 10.1038/s41593-023-01536-9

PubMed Abstract | Crossref Full Text | Google Scholar

Keywords: calcium imaging, coarse graining, electrophysiology, neural criticality, population dynamics, renormalization group

Citation: Li MS, Liu B, van Antwerp KW, Abdelaleem E and Sederberg AJ (2026) Measurement effects on critical scaling in neural systems. Front. Comput. Neurosci. 19:1724190. doi: 10.3389/fncom.2025.1724190

Received: 13 October 2025; Revised: 10 December 2025;
Accepted: 29 December 2025; Published: 23 January 2026.

Edited by:

Tiago Ribeiro, National Institutes of Health (NIH), United States

Reviewed by:

Antonio Jorge Fontenele, University of Arkansas, United States
Yin-Jui Chang, The University of Texas at Austin, United States

Copyright © 2026 Li, Liu, van Antwerp, Abdelaleem and Sederberg. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Audrey J. Sederberg, YXVkcmV5LnNlZGVyYmVyZ0BnYXRlY2guZWR1

Present Address: M. Shane Li, Department of Philosophy, University of Cincinnati, Cincinnati, OH, United States

Disclaimer: All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.