## Abstract

What information single neurons receive about general neural circuit activity is a fundamental question for neuroscience. Somatic membrane potential (*V*_{m}) fluctuations are driven by the convergence of synaptic inputs from a diverse cross-section of upstream neurons. Furthermore, neural activity is often scale-free, implying that some measurements should be the same, whether taken at large or small scales. Together, convergence and scale-freeness support the hypothesis that single *V*_{m} recordings carry useful information about high-dimensional cortical activity. Conveniently, the theory of “critical branching networks” (one purported explanation for scale-freeness) provides testable predictions about scale-free measurements that are readily applied to *V*_{m} fluctuations. To investigate, we obtained whole-cell current-clamp recordings of pyramidal neurons in visual cortex of turtles with unknown genders. We isolated fluctuations in *V*_{m} below the firing threshold and analyzed them by adapting the definition of “neuronal avalanches” (i.e., spurts of population spiking). The *V*_{m} fluctuations which we analyzed were scale-free and consistent with critical branching. These findings recapitulated results from large-scale cortical population data obtained separately in complementary experiments using microelectrode arrays described previously (Shew et al., 2015). Simultaneously recorded single-unit local field potential did not provide a good match, demonstrating the specific utility of *V*_{m}. Modeling shows that estimation of dynamical network properties from neuronal inputs is most accurate when networks are structured as critical branching networks. In conclusion, these findings extend evidence of critical phenomena while also establishing subthreshold pyramidal neuron *V*_{m} fluctuations as an informative gauge of high-dimensional cortical population activity.

**SIGNIFICANCE STATEMENT** The relationship between membrane potential (*V*_{m}) dynamics of single neurons and population dynamics is indispensable to understanding cortical circuits. Just as important to the biophysics of computation are emergent properties such as scale-freeness, where critical branching networks offer insight. This report makes progress on both fronts by comparing statistics from single-neuron whole-cell recordings with population statistics obtained with microelectrode arrays. Not only are fluctuations of somatic *V*_{m} scale-free, they match fluctuations of population activity. Thus, our results demonstrate appropriation of the brain's own subsampling method (convergence of synaptic inputs) while extending the range of fundamental evidence for critical phenomena in neural systems from the previously observed mesoscale (fMRI, LFP, population spiking) to the microscale, namely, *V*_{m} fluctuations.

- balanced networks
- membrane potential
- neural computation
- neuronal avalanches
- renormalization group
- scale-free

## Introduction

How do cortical population dynamics impact single neurons? What can we learn about cortical population dynamics from single neurons? These questions are central to neuroscience. Uncovering the functional significance of multiscale organization within cerebral cortex requires knowing the relationship between the dynamics of networks and individual neurons within them (Nunez et al., 2013).

For pyramidal neurons in the visual cortex, somatic spike generation is ambiguously related to presynaptic firing (Tsodyks and Markram, 1997; Brunel et al., 2014; Gatys et al., 2015; Stuart and Spruston, 2015; Moore et al., 2017). Such neurons pass spiking information to many postsynaptic neurons (Lee et al., 2016). However, a presynaptic pool with multifarious neighboring and distant neurons (Hellwig, 2000; Wertz et al., 2015) provides excitatory and inhibitory synaptic inputs throughout the soma and complex dendritic architecture (Magee, 2000; Larkum et al., 2008; Moore et al., 2017). Input propagation to the axon hillock has both active and passive features (London and Häusser, 2005) and the membrane potential (*V*_{m}) response is increasingly nonlinear near the action potential threshold. Thus, such details of network propagation give *V*_{m} more utility than focusing solely on spiking.

Most computational neuroscientists use spiking data because spikes are “the currency of the brain” (Wolfe et al., 2010) and extracellular recording is straightforward compared with whole-cell recording. Yet, the paucity of single-neuron spiking (Shoham et al., 2006) and limited foreknowledge about connections (Helmstaedter, 2013) make extracellular single-unit observation an impoverished means of studying neuronal circuits. In contrast, subthreshold *V*_{m} fluctuations contain rich information about the circuits containing each neuron (Sachidhanandam et al., 2013; Petersen, 2017). Integral to gaining a neuron's view of the brain is uncovering relationships between the statistics of *V*_{m} fluctuations and fluctuations of local spiking and then contrasting against other plausible one-dimensional signals.

We look for such relationships in the strict predictions and rigorous measurements of scale-freeness used to identify a fragile network connectivity pattern known as “critical branching.” This pattern exhibits emergent properties valuable for information processing, such as higher susceptibility and dynamic range (Haldeman and Beggs, 2005; Beggs, 2008; Shew and Plenz, 2013; Shriki and Yellin, 2016; Timme et al., 2016), but omits some neuronal dynamics (Poil et al., 2008, 2012) without extension (Porta and Copelli, 2018). The pattern is as follows: on average over all neuronal avalanches (spiking above baseline; Friedman et al., 2012), one spike leads to exactly one other spike. In most arbitrary networks there is less or more than one; these are “subcritical” and “supercritical” respectively. Among the dazzling emergent properties of “criticality” are universality, self-similarity, and scale-free correlations (Stanley, 1999).

These are as follows: A “universality class” is a set of incongruous systems exhibiting identical statistics only at their “critical points.” “Self-similarity” includes fractal patterns and power laws in geometrical analysis of avalanches (power laws are “scale-invariant,” popularly called “scale-free”). Avalanches of any duration have identical average shapes after normalization (Shaukat and Thivierge, 2016). Avalanche areas grow with duration as another power law (Sethna et al., 2001). However, observation methods must be consistent with event propagation (Priesemann et al., 2009; Yu et al., 2014; Levina and Priesemann, 2017). Additionally, pairwise correlation versus length or time are also power laws (Chialvo, 2010), meaning any input has a nonzero chance of propagating forever or anywhere.

In summary, the theory of critical branching networks offers superb standards of comparison for three reasons: neuronal avalanche analysis applies to *V*_{m}, offers promising insights, and makes precise predictions about fluctuation geometry. We study both *V*_{m} fluctuations and criticality with one simple question: Do *V*_{m} fluctuations match the scale-free statistics of cortical populations (see Fig. 1)?

To address this question, we simultaneously recorded somatic *V*_{m} from pyramidal neurons and local field potential (LFP) in visual cortex and performed avalanche analysis on fluctuations. We found that subthreshold *V*_{m} fluctuation statistics match published microelectrode array (MEA) data. We used surrogate testing to show why negative LFP fluctuations don't match and modeling to demonstrate dependence on critical branching.

## Materials and Methods

### Surgery and visual cortex

All procedures were approved by Washington University's Institutional Animal Care and Use Committees and conform to the National Institutes of Health's Guide for the Care and Use of Laboratory Animals. Fourteen adult red-eared sliders (*Trachemys scripta elegans*, 150–1000 g) were used for this study; their genders were not recorded. Turtles were anesthetized with propofol (2 mg of propofol/kg) then decapitated. Dissection proceeded as described previously (Saha et al., 2011; Crockett et al., 2015; Wright et al., 2017a).

To summarize, immediately after decapitation, the brain was excised from the skull with the right eye intact and bathed in cold extracellular saline containing the following (in mm): 85 NaCl, 2 KCl, 2 MgCl_{2}*6H_{2}O, 20 dextrose, 3 CaCl_{2}-2H_{2}O, 45 NaHCO_{3}. The dura was removed from the left cortex and right optic nerve and the right eye hemisected to expose the retina. The rostral tip of the olfactory bulb was removed, exposing the ventricle that spans the olfactory bulb and cortex. A cut was made along the midline from the rostral end of the remaining olfactory bulb to the caudal end of the cortex. The preparation was then transferred to a perfusion chamber (Warner RC-27LD recording chamber mounted to PM-7D platform) and placed directly on a glass coverslip surrounded by Sylgard. A final cut was made to the cortex (orthogonal to the previous and stopping short of the border between medial and lateral cortex), allowing the cortex to be pinned flat with the ventricular surface exposed. Multiple perfusion lines delivered extracellular saline to the brain and retina in the recording chamber (adjusted to pH 7.4 at room temperature).

We used a phenomenological approach to identify the visual cortex, described previously (Shew et al., 2015). In brief, this region was centered on the anterior lateral cortex, in agreement with voltage-sensitive dye studies (Senseman and Robbins, 1999, 2002). Anatomical studies identify this as a region of cortex receiving projections from lateral geniculate nucleus (Mulligan and Ulinski, 1990). We further identified a region of neurons as belonging to the visual cortex when the average LFP response to visual stimulation crossed a given threshold and patched within that neighborhood (radius of ∼300 μm).

### Intracellular recordings

For whole-cell current-clamp recordings, patch pipettes (4–8 MΩ) were pulled from borosilicate glass and filled with a standard electrode solution containing the following (in mm): 124 KMeSO_{4}, 2.3 CaCl_{2}-2H_{2}O, 1.2 MgCl_{2}, 10 HEPES, and 5 EGTA adjusted to pH 7.4 at room temperature. Cells were targeted for patching using a differential interference contrast microscope (Olympus). *V*_{m} recordings were collected using an Axoclamp 900A amplifier, digitized by a data acquisition panel (National Instruments PCIe-6321), and recorded using a custom LabVIEW program (National Instruments), sampling at 10 kHz. As described in (Crockett et al., 2015; Wright and Wessel, 2017; Wright et al., 2017a,b), before recording from a cell after initial patching current was injected to elicit spiking. This current injection test was also repeated intermittently between recording trials. Recording did not proceed if a cell spiked inconsistently (e.g., failure to spike, insufficient spike amplitude) in response to injected current or exhibited extreme depolarization in response to small current injection amplitudes. If a clog or loss of seal was suggested by unusually erratic *V*_{m} short timescales current, the current injection test was performed and, upon failure, the affected recording was marked for exclusion from analysis. We excluded cells that did not display stable resting *V*_{m}s for long enough to gather enough avalanches. Up to three whole-cell recordings were made simultaneously. In total, we obtained recordings from 51 neurons from 14 turtles.

Recorded *V*_{m} fluctuations taken in the dark (no visual stimulation) were interpreted as ongoing activity. Such ongoing cortical activity was interrupted by visual stimulation of the retina with whole-field flashes and naturalistic movies as described previously (Wright and Wessel, 2017; Wright et al., 2017a,b). An uninterrupted recording of ongoing activity lasted for 2–5 min. Periods of visual stimulation were too short and were too frequently interrupted by action potentials to yield the great number of avalanches required for rigorous power-law fitting.

A sine-wave removal algorithm was used to remove 60 Hz line noise. Action potentials in turtle cortical pyramidal neurons are relatively rare. An algorithm was used to detect spikes and the *V*_{m} recordings between spikes were extracted and filtered from 0 to 100 Hz. *V*_{m} recordings were detrended by subtracting the fifth percentile in a sliding 2 s window. The resulting signal was then shifted to have the same mean value as before subtraction. Detrending did not affect the size of *V*_{m} fluctuations (data not shown).

### Extracellular recordings

Extracellular recordings were achieved with tungsten microelectrodes (microprobes heat-treated tapered tip) with ∼0.5 MΩ impedance. Electrodes were slowly advanced through tissue under visual guidance using a manipulator (Narishige) while monitoring for activity using custom acquisition software (National Instruments). The extracellular recording electrode was located within ∼300 μm of patched neurons. Extracellular activity was collected using an AM Systems Model 1800 amplifier, band-pass filtered between 1 Hz and 20,000 Hz, digitized (NI PCIe-6231), and processed using custom software (National Instruments). Extracellular recordings were downsampled to 10,000 Hz and then filtered (100 Hz low-pass), yielding the LFP. The LFP was filtered and detrended as described above (see “Intracellular recordings” section) except that the mean of the entire signal was subtracted and the signal was multiplied by −1 before it was detrended. This final inverted signal is commonly featured in the literature as negative LFP or nLFP (Kelly et al., 2010; Kajikawa and Schroeder, 2011; Okun et al., 2015; Ness et al., 2016).

### Experimental design and statistical analysis

#### Setwise comparisons

To measure differences between sets of statistics, we relied on three nonparametric measures. We used the MATLAB Statistics and Machine Learning Toolbox implementation of Fisher's exact test (Hammond et al., 2015). This allowed us to measure the effect size (odds ratio *r*_{OR}) and statistical significance (*p*-value) of finding that consistency-with-criticality is more frequent or less frequent in an experimental group than a control group.

To quantify the similarity between the exponents measured in different sets of data, we used the MATLAB Statistics and Machine Learning Toolbox implementations of the exact Wilcoxon rank-sum test (Hammond et al., 2015) and the exact Wilcoxon signed-rank test. In both cases effect size, *r*_{SDF} is measured by the simple difference formula (Kerby, 2014). The rank-sum test is used when comparing nonsimultaneous recordings, such as comparing MEA data with *V*_{m} data. The signed-rank test is used when comparing data that can be paired, such as *V*_{m} data to concurrent LFP. When comparing whether a dataset differs from a specific value, we can use the sign test.

The significance level was set at *p* = 0.05 for all tests. Each setwise comparison test stands alone as its own conclusion. They were not combined to assess the significance of any effect across sets-of-sets. Thus, we are not making multiple comparisons and no corrections are warranted (Bender and Lange, 2001).

#### Random surrogate testing

It is possible that scale-free observations have an origin in independent random processes of a type demonstrated previously (Touboul and Destexhe, 2017). To control for this, we phase-shuffled the *V*_{m} fluctuations using the amplitude-adjusted Fourier transform (AAFT) algorithm (Theiler et al., 1992). This tests against the null hypothesis that a measure on a time series can be reproduced by performing a nonlinear rescaling of a linear Gaussian process with the same autocorrelation (same Fourier amplitudes) as the original process. Phase information is randomized, which removes higher-order correlations but preserves the scale-free power-spectrum.

The AAFT tests only higher-order correlations, but a simpler algorithm tests against the null hypothesis that an un-rescaled linear Gaussian process with the same autocorrelation as the original process can produce the same results (Theiler et al., 1992). This is known as the unwindowed Fourier transform (UFT). Once we see what measures depend on the higher-order correlations with the AAFT, we can use the UFT to see how measures depend on the non-Gaussianity (nonlinear rescaling), which is inherent to excitable membranes. Using the UFT alone would make it difficult to attribute whether statistically significant differences are due to the rescaling or to the higher-order correlations (Rapp et al., 1994).

We performed AAFT and UFT on each *V*_{m} time series once and then compared how the two datasets performed on every metric used in this study. The datasets were compared with a matched Wilcoxon sign-rank test implemented via MATLAB's statistics tool box. Doing the comparison at a dataset level allowed us to obtain a discrimination statistic for every metric that we used without repeating the computationally expensive analysis procedure hundreds or thousands of times on every *V*_{m} trace. With enough individual recordings in each dataset, the matched Wilcoxon sign-rank test is a reliable measure that empowered us to efficiently compare all important metrics.

#### Neuronal avalanche analysis

Neuronal avalanches were defined by methods analogous to those described previously (Poil et al., 2012), which are used for uninterrupted ongoing signals; conversely, methods based on event detection (Beggs and Plenz, 2003) require periods of nonactivity. A threshold is defined and an avalanche starts when the signal crosses the threshold from below and ends when the signal crosses the threshold from above. The choice of threshold is a free parameter and we set it to the 25th percentile before conducting the complete analysis. In similar situations (continuous nonzero signals), researchers chose one-half the median (Poil et al., 2012; Larremore et al., 2014). However, one-half the median cannot work for negative signals or signals with high mean but low variance. Before analysis, threshold choices between the 15^{th} to 50^{th} percentile were tested on data from the five cells with the most recordings to see how threshold may affect the number of avalanches. The 25^{th} percentile was consistent with the existing literature and gave many avalanches compared with alternatives. Having a large number of avalanches is important because it gives the best statistical resolution. An analysis with a choice of threshold that yields fewer avalanches (or changing the threshold for each recording) would be suspect for selecting serendipitous results. After the analysis was conducted, eight percentiles between the 15^{th} to 50^{th} percentile were tested and gave similar power-law exponents.

We quantified each neuronal avalanche by its size *A* and its duration *D*. The avalanche size is the area between the processed *V*_{m} recording and the baseline. The baseline is another free parameter that was set at the second percentile of the processed *V*_{m} recording. The second percentile was chosen because its value is more stable than the absolute minimum. The avalanche duration *D* is the time between threshold crossings.

The lower limit of avalanche duration is defined by the membrane time constant, which has been reported to be between 50 and 140 ms for the turtle brain at room temperature (Ulinski, 1990; Larkum et al., 2008). We took a conservative approach by setting the limit at less than half the lower bound on membrane time constant, which was significantly less than the lower cutoff from power-law fits. Only avalanches with a duration longer than 20 ms were included in the analysis. Thus, we avoided artificially retaining only the events most likely to be power-law distributed.

Following the procedure described above, each processed *V*_{m} recording of uninterrupted ongoing activity (i.e., a recording of 2–5 min duration) yielded 327 ± 148 (mean ± SD) avalanches. This is insufficient for rigorous statistical fitting on recordings individually (Clauset et al., 2009). Therefore, we grouped avalanches from multiple recordings of ongoing activity of the same cells. Each cell produced between 3 and 19 recordings of ongoing activity (2–5 min duration each recording), with trials recorded intermittently over a period of 10–60 min. We grouped recordings based on whether they occurred in the first or second 20 min period since the beginning of recording from that neuron. Then, all the avalanches from the first or second 20 min period were grouped together with one data object (the group) storing the size and duration of each avalanche. It is rare for neurons to have recordings in the third 20 min periods, so these data were not included. Since there was a slow drift in the mean *V*_{m} over a period of several minutes, we scaled the avalanche sizes from each recording to have the same median as other recordings from the same group. *Z*-scoring was not useful for accounting for trial to trial variability because it does affect whether a specific time point is above or below a certain percentile threshold. Therefore, it is not useful for removing variability in avalanche duration. Windowed *z*-scoring introduces artifacts near action potentials. On average, four recordings were possible in each 20 min period. There were 51 neurons with multiple recordings of ongoing activity in the first 20 min of experimentation (thus 51 recording groups); of these, 18 neurons had an additional 20 min period with more than one recording. This produced a total of 69 groups with 1346 ± 1018 (mean ± SD) avalanches for each group. Of these 69 groups, 57% had >1000 avalanches. The largest number of avalanches was 7495 and the smallest was 313. Only five groups had <500 avalanches. We report on the 51 groups from the first 20 min period separately from the 18 groups with recordings from the second 20 min period of experimentation.

For each group, we evaluated the avalanche size and duration distributions with respect to power laws. To test whether a distribution followed a power law, we applied the rigorous statistical fitting routine described previously (Clauset et al., 2009). We tested three power-law forms: *P*(*x*) ∝ *x*^{−α} (with and without truncation) (Deluca and Corral, 2013), as well as a power law with exponential cutoff, *P*(*x*) ∝ *x*^{−α}*e*^{−x/r}. We compared these against log normal and exponential alternative (non-power-law) hypotheses. Distribution parameters were estimated using maximum likelihood estimation (MLE) and the best model out of those fitted to the data was chosen using the Akaike information criterion (Bozdogan, 1987). It should be acknowledged that a small power-law region in the truncated form would be suspect for false positives, likewise for a strong exponential cutoff (Deluca and Corral, 2013). Finally, to decide whether a fitted model was plausible, pseudorandom datasets were drawn from a distribution with the estimated parameters and then the fraction that had a lower fit quality (Kolmogorov–Smirnov distance) than the experimental data was calculated. If this fraction, called the comparison quotient *q*, was >0.10, then the best-fit model (according to the Akaike information criterion) was accepted as the best candidate. Otherwise, the next best model was considered.

We applied several additional steps and strict criteria to control for false positives. One such step was assessing whether the scaling relation was obeyed over the whole avalanche distribution for each group (not just the portion above the apparent onset of power-law behavior). The scaling relation is another power-law, 〈*A*〉(*D*) ∝ *D*^{γ}, predicting how the measured size of avalanches increases geometrically with increasing duration (on average). For any dataset with three power laws, 〈*A*〉(*D*) ∝ *D*^{γ} (scaling relation), *P*(*A*) ∝ *A*^{−τ} (size distribution), and *P*(*D*) ∝ *D*^{−β} (duration distribution), the scaling relation exponent is predicted by the other two exponents by γ ≈ γ* _{p}* = (Scarpetta et al., 2018). Note that γ

*= 1 is a trivial value because it implies 〈*

_{p}*A*〉(

*D*) ∝

*D*and that would suggest that individual avalanches were just noise symmetric about a constant value. This would mean that the average avalanche shape is just a flat line at some constant of proportionality, =

*a*, where is a function describing the shape of an avalanche of duration

*D*,

*t*

_{0}is the beginning of the avalanche, and

*a*is a constant.

##### Standards for consistency with critical point behavior.

We applied four standardized criteria to provide a transparent and systematic way to produce a binary classification; either “no inconsistencies with activity near a critical point were detected” or “some inconsistencies with activity near a critical point were detected.”

First, a collection of avalanches must be power-law distributed in both its size and duration distributions.

Second, the collection of avalanches must have a power-law scaling relation as determined by *R*^{2} > 0.95 (coefficient of determination) for linear least-squares regression to a log-log plot of average size vs durations: log (*A*(*D*)) ∼ γ log (*D*) + *b*. This *R*^{2} represents the best that any linear fit can achieve and must include all the avalanches, not a subset. We denote the scaling exponent (slope from linear regression) from this fit as γ_{f}.

Third, the scaling relation exponent predicted by theory (denoted as γ_{p}) must correspond to a trendline on a log-log scatter plot of 〈*A*〉(*D*) with an *R*^{2} that is within 90% of the best-case fitted trendline from the second criterion. Again, the *R*^{2} for the predicted scaling relation is calculated across all avalanches, not just the subset above the inferred lower cutoff of power-law behavior (which was found for the first criterion). This cross-validates agreement with theory.

Fourth, the fitted scaling relation exponent must be significantly greater than 1: (γ* _{f}* − 1) > σ

_{γf}where σ

_{γf}is the SD. This last requirement eliminates scaling that might be trivial in origin. It is measured after getting the fitted scaling relation exponent for all of the data so that a dataset SD can be determined. It is necessary to also check that the set of scaling relation exponents from the power-law fits to all avalanche sets is significantly different from 1 at a dataset level. A scaling relation exponent equal to one suggests a linear relationship between mean size and duration that is not consistent with criticality in neural systems (Haldeman and Beggs, 2005).

Our four-criterion test cannot measure distance from a critical point nor eliminate all risk of false positives. To complete our analysis, we also looked at three additional factors: (1) whether exponent values match exponent values from other experiments as expected from the universality prediction of theory, (2) whether all the exponents within our dataset have similar scaling relation predictions, and (3) whether the avalanches within our dataset exhibit shape collapse across all the recordings.

##### Obtaining shape collapse and analyzing it quantitatively and qualitatively.

Shape collapse is a very literal manifestation of scale-invariance (also called “self-similarity”) (Sethna et al., 2001; Beggs and Plenz, 2003; Friedman et al., 2012; Pruessner, 2012; Timme et al., 2016). Avalanches of different durations should rise and fall in the same way on average. This average avalanche profile is called a scaling function. The average avalanche profile for avalanches of duration *D* is predicted to be *A* (*t*,*D*) = *D*^{(γ−1)} where *D*^{(γ−1)} is the power-law scaling coefficient that modulates the height of the profile and is the universal scaling function itself (normalized in time). Shape collapse analysis provides an independent estimate of the scaling relation exponent γ_{SC}, which is only expected to be accurate at criticality (Sethna et al., 2001; Scarpetta and de Candia, 2013; Shaukat and Thivierge, 2016), and a visual test of conformation to an empirical scaling function.

Exponent estimation is very sensitive to the unrelated, intermediate rescaling steps involved in combining the avalanches from multiple recordings into one group. To get an estimate of the scaling relation exponent for each group, γ_{SC}, we averaged the scaling exponents, γ_{i}, found individually for each recording in that group where *i* denotes the *i*^{th} recording and SC is shape collapse).

Naturally, individual avalanche profiles are vectors of variable length *D*. We must first “rescale in time” to make them vectors of equal length without losing track of what each vector's original duration was. We do that by linearly interpolation with 20 evenly spaced points. So, the *j*^{th} avalanche profile of the *i*^{th} recording is denoted as a 20-element vector where the top arrow denotes a vector.

Next, the set of all profiles from recording *i* with the exact same duration *D*, denoted as **Γ_{Di}** where bold indicates a set, were averaged and divided by a test scaling factor

*D*

^{(γ′i−1)}. We define this as (γ′) = 〉

*D*

^{−(γ′i−1)}. The prime indicates a test rescaling. The average is over all vectors in the set

**Γ**. The choice of γ

_{Di}_{i}was optimized using MATLAB's fminsearch function to minimize the mean relative error between the average over all durations 〈(γ′)〉 and the set members (γ′) so that for recording

*i*: This error minimization and applying the rescaling is the “collapse” in “shape collapse.”

Once we have the γ_{i} for the avalanches in each individual recording of ongoing activity, we compare the average, γ* _{SC}* = 〈γ

*〉, with the predicted and fitted scaling relation exponents for the group of recordings, γ*

_{i}_{p}and γ

_{f}(statistical comparison tests are described in a previous section). Thus, quantitative analysis of shape collapse was done by comparing γ

_{SC}, γ

_{p}, and γ

_{f}for each of the 69 groups individually.

Visual assessment of how well avalanche profiles can be described by one universal scaling function, supports the quantitative exponent estimation. This was performed by averaging all the profiles within specific duration bins (regardless of trial or group) and plotting them on top of one another. Shape collapse always requires a very large number of avalanches, so we had to combine avalanches from all 69 groups. However, the resting *V*_{m} differs from recording to recording and cell to cell. Therefore, avalanche profiles from different recordings are vertically misaligned. To combine avalanche profiles from different recordings, we divided all the profiles by a scalar value unique to each recording: the time average over all the collapsed profiles. This produced rescaled and mean-shifted profiles, denoted by a double prime, = /〈**Γ′_{ijk}**〉, where

*k*∈ [1,20] denotes the interpolated time point. The set of avalanches from each recording were thus aligned, but individual variability was preserved and thus profiles from different recordings could be averaged without introducing artifacts. This set,

**Γ″**, contained a total of 106,220 shifted and rescaled profiles for the

_{ij}*V*

_{m}data.

The set of shifted and rescaled profiles falling into a duration bin is denoted **Γ″_{D}**. Each duration bin then provides its own estimate of the scaling function ∼ . For each bin,

*D*was defined as the average duration of all constituent profiles. If <700 avalanches had a particular duration, we included the next longest duration iteratively until we met or exceeded 700 avalanches. This only applied to long durations. The choice of 700 was made because it allowed smooth averaging without excessively wide duration bin widths.

We also assessed the mean curvature of avalanche profiles from the rescaled profile for a particular duration . This allowed us to plot how curvature depends on duration. Mean curvature 〈κ〉 is defined as follows, where *k* still denotes time points:

### Model simulations

We simulated a model network consisting of *N* = 10^{4} binary probabilistic model neurons. The model neurons form a directed random network (Erdős–Rényi random graph) where the probability that neuron *j* connects to neuron *i* is *c*. In a network of *N* neurons, this results in a mean in-degree and out-degree of *cN*. We tested nine not quite evenly distributed values of connection probabilities, *c* ∈ [0.5, 1, 3, 5, 7.5, 10, 15, 20, 25] × 10^{−2}. As discussed previously (Kinouchi and Copelli, 2006; Larremore et al., 2011a, 2014), the impact of connectivity on network dynamics is nonlinear, so we took a finer look at smaller connection probabilities while maintaining thorough coverage of intermediate connection probabilities.

The strength of the connection from neuron *j* to neuron *i* is quantified in terms of the network adjacency or weight matrix, *W*, with the fortune of having a simple and intuitive meaning. For each existing connection from neuron *j* to neuron *i*, *W*_{ij} is the direct change in the probability that neuron *i* will fire at the next time step if neuron *j* spikes in the current time step.

The dynamics of this network are well characterized by the largest eigenvalue λ of the network weight matrix *W* with criticality occurring at λ = 1 (Kinouchi and Copelli, 2006; Larremore et al., 2011a,b, 2012, 2014). The physical interpretation of λ is a “branching parameter”(Haldeman and Beggs, 2005) that governs expected number of spikes immediately caused by the firing of one neuron. If λ = 1, then one spike causes one other spike on average, while, if λ > 1,then one spike causes more than one on average and vice versa.

We tested five different values of largest eigenvalue at, near and far from criticality λ ∈ [0.9, 0.95, 1, 1.015, 1.03]. A fraction, χ, of the neurons is designated as inhibitory. This is done by multiplying all outgoing connections of an inhibitory neuron by −1. We tested nine different values of the fraction of inhibitory neurons in the range from 0 to 0.25, thus including the value 0.2, corresponding to the fraction of inhibitory neurons in the mammalian cortex (Meinecke and Peters, 1987). The magnitudes of nonzero weights are independently drawn from a distribution of positive numbers with mean η, where the distribution is uniform on [0,2η] and η is given by η = λ/(*cN*(1 − 2χ)). The maximum eigenvalue is then fine-tuned by dividing *W* by the current maximum eigenvalue and set to the exactly desired value *W* = λ*W′*/λ′ where *W′* and λ′ are the matrices and eigenvalues before correction.

The binary state *S*_{i}(*t*) of neuron *i* at time *t* denotes whether the model neuron spikes (*S _{i}*(

*t*) = 1) or does not spike (

*S*(

_{i}*t*) = 0) at time

*t*. At each time step, the states of all neurons are updated synchronously according to the following update rule: Where ξ

_{i}(

*t*) is a random number on [0,1] drawn from a uniform distribution, and Θ is the Heaviside step function. In addition to this update rule, a refractory period of one time step (translated to ∼2 ms) was imposed for certain parameter conditions. A simulation begins with initiating the activity of one randomly chosen excitatory neuron and continuing the simulation until overall network activity had ceased. The process was then repeated.

From the simulated binary states of 10^{4} model neurons, we extracted three measures of simulated activity. First, the network activity *F*(*t*) = ∑_{i=1}^{N}*S _{i}*(

*t*) ̸

*N*is the fraction of neurons spiking at time

*t*. Second, the input to model neuron

*i*at time

*t*is

*P*(

_{i}*t*) = ∑

*(*

_{j}^{N}W_{ij}S_{j}*t*− 1), which is almost always positive for our parameters. Note that

*P*(

_{i}′*t*) =

*P*(

_{i}*t*) × Θ(

*P*(

_{i}*t*)) directly represents the probability for the neruon i to spike at time

*t*. Third, we constructed a proxy for the

*V*

_{m}signal, Φ

*(*

_{i}*t*) = (α

***

_{h}*P*)(

_{i}*t*), by convolving the input

*P*

_{i}(

*t*) with an alpha-function: α

*(*

_{h}*t*) = exp with

*h*

_{m}= 2 time steps (assumed to be ∼4 ms).

A total of 405 different parameter combinations (connection density, inhibition, maximum eigenvalue) were simulated. Each combination was simulated 10 times. Based on the connection probability *c* and the fraction of inhibition χ, we distinguish four regions in parameter space classified according to the behavior of the critical model; that is, λ = 1.

The first region is the “positive weights” region. Without inhibition activity increases or dies out in accordance with the branching parameter. This region is defined by χ = 0. With moderate inhibition and dense connectivity, there is a region of parameter space we call “quiet”; activity lasts only slightly longer than in a system with no inhibition. This region is defined by the *ex post facto* boundaries *c* ≥ *e*^{11χ}/25 and χ > 0. Further increasing inhibition relative to connection density produces a behavior like “up and down” states (or “telegraph noise”) (Sachdev et al., 2004; Millman et al., 2010). We call this the “switching” regime because network activity switches between a low mean and a high mean. This region is defined by *c* < *e*^{11χ}/25 and *c* ≥ (10*e*^{12χ} − 13)/100 and χ > 0. When inhibition is high relative to connection density, the system enters the “ceaseless” region where stimulating one neuron causes activity that effectively never dies out. An especially attractive feature of this model is that the “ceaseless” and “switching” regimes exhibit sustained self-generated activity. This provides a way to model spontaneous neural activity without externally imposed firing patterns.

Refractoriness was studied in the network without inhibition and it was found that dynamic range was inversely proportional to refractory period (Larremore et al., 2011a), but the empirical branching parameter (criticality) displayed no dependence on refractory period (Kinouchi and Copelli, 2006). In studies that featured inhibition and introduced ceaselessness, no refractory was used (Larremore et al., 2014). However, we found that, for some networks in the switching regime, the maximum eigenvalue was a better predictor of the empirical branching ratio if the refractory period was one time step. Because this relationship is central to our understanding of criticality in this model, we ran an initial testing cycle before each simulation begins to decide whether to set the refractory period to one time step or to zero. Doing so ensures that the network displays critical-like phenomena in all regimes (the maximum eigenvalue of connectivity), but also ensures that the model is consistent with prior studies.

We performed avalanche analysis on each of the simulated signals using the methods described above for *V*_{m} recordings. If the network was in the switching regime, then we only performed analysis on the periods when the network was in the mode (high or low mean) in which it spent the majority of its time. As before, the 25th percentile defined the avalanche threshold. If the signal had negative values, as in the case of single neuron *V*_{m} proxies in networks with inhibition, then the signal was shifted by subtracting the second percentile. To obtain good statistics, we continued stimulating and extracting avalanches until a simulation either reached 10^{4} avalanches or 5 × 10^{3} avalanches and a very large file size or a very long computational time. This ensured that there were between 2000 and 10,000 avalanches per trial.

### Data and software accessibility

All raw data are available at https://github.com/jojker/continuous_signal_avalanche_analysis and the code developed for this analysis is available upon request to the corresponding author.

## Results

Single-neuron *V*_{m} fluctuations are thought to be dominated by synaptic inputs from multitudes of presynaptic neurons (Stepanyants et al., 2002; Brunel et al., 2014; Petersen, 2017). Because the way neurons integrate their diverse inputs is central to information processing in the brain, it is important that neuroscience gain a thorough understanding of the relationship between subthreshold *V*_{m} fluctuations and population activity. A basic step is to compare statistical analyses, especially analyses in which a meaningful relationship is expected. We investigated whether an avalanche analysis on *V*_{m} fluctuations would reveal the same signatures of scale-freeness and critical network dynamics found in measures of population activity (Fig. 1) (Friedman et al., 2012; Shew et al., 2015; Marshall et al., 2016). To address this comparison across organizational levels, we recorded *V*_{m} fluctuations from 51 pyramidal neurons in visual cortex of 14 turtles and assessed evidence for critical network dynamics from these recordings.

In a model investigation we corroborated results evaluated the conditions needed to enable inferring dynamical network properties from the inputs to single neurons. Finally, we extended the analysis to other commonly recorded time series of neural activity for comparison with the information content of *V*_{m} fluctuations about the dynamical network properties.

*V*_{m} fluctuations reveal signatures of critical point dynamics

We obtained whole-cell recordings from pyramidal neurons in the visual cortex of the turtle *ex vivo* eye-attached whole-brain preparation (Fig. 2*A*). Recorded *V*_{m} fluctuations taken in the dark (no visual stimulation) were interpreted as ongoing activity. We analyzed the recorded ongoing *V*_{m} fluctuations employing the concept of “neuronal avalanches” (Beggs and Plenz, 2003; Poil et al., 2012; Shew et al., 2015), which are positive fluctuations of network activity. For continuous time series such as the *V*_{m} recording, one selects a threshold and a baseline. We defined a neuronal avalanche based on the positive threshold crossing followed by a negative threshold crossing of the *V*_{m} time series (Poil et al., 2012; Hartley et al., 2014; Larremore et al., 2014; Karimipanah et al., 2017a). We quantified each neuronal avalanche by its size, *A*, defined as the area between the curve and the baseline, and its duration *D*, defined as the time between threshold crossings (Fig. 2*B*).

To quantify the statistics of avalanche properties, we applied concepts and notations from the field of “critical phenomena” in statistical physics (Nishimori and Ortiz, 2011; Pruessner, 2012). Because the critical point is such a small target for any naturally occurring self-organization (Pruessner, 2012; Hesse and Gross, 2014; Cocchi et al., 2017) and there is considerable risk of false positives (Taylor et al., 2013; Hartley et al., 2014; Touboul and Destexhe, 2017; Priesemann and Shriki, 2018), asserting criticality in a new system or with a new tool requires extraordinary evidence. Because this is a new tool, we created four criteria and set quantifiable standards for concluding a system is consistent with criticality based on avalanche power laws and completed this exhaustive battery of tests with shape collapse, a geometrical analysis of self-similarity in the avalanche profiles (see “Experimental design and statistical analysis” section).

In brief, we found that both the size and duration distributions of the fluctuations treated as avalanches were consistent with power laws (Fig. 2*C*), *P*(*A*) ∝ *A*^{−τ} and *P*(*D*) ∝ *D*^{−β} matching widely reported exponents (Beggs and Plenz, 2003; Priesemann et al., 2009, 2014; Hahn et al., 2010; Klaus et al., 2011; Friedman et al., 2012; Shriki et al., 2013; Arviv et al., 2015; Shew et al., 2015; Karimipanah et al., 2017a,b), obeyed the scaling relation (Fig. 2*D*), and exhibited shape collapse over an expansive set of durations (Fig. 2*E*).

Specifically, of the 51 recording groups featuring data from the first 20 min period of recording from one cell, 98% had power laws in both size and duration distributions. The exponent values for the size distribution were τ = 1.91 ± 0.38 (median ±SD). Exponent values for the duration distribution were β = 2.06 ± 0.48. Of the 51 neurons with a recording group from the first 20 min, 18 had an additional 20 min period spanning multiple recordings. All these 18 groups had power laws in both size and duration; the exponent values for the size distribution were τ = 1.87 ± 0.29 and the exponent values for the duration distribution were β = 2.21 ± 0.39.

It is also important to confirm that power-law behavior extends across several orders of magnitude of avalanche durations. We demonstrate a power-law distribution over 2.45 ± 0.39 orders of magnitude of duration. For the scaling relation, we found a larger span with 2.62 ± 0.23 orders of magnitude across our whole avalanche duration range.

Another statistic crucial to signatures of criticality measures the relationship between the power laws describing size and duration of avalanches (Sethna et al., 2001; Beggs and Timme, 2012; Friedman et al., 2012). If the average avalanche size also scales with duration according to 〈*A*〉(*D*) ∝ *D*^{γ}, then the exponent γ is not independent, but rather depends on the exponents τ and β according to γ = (β − 1)̸(τ − 1) regardless of criticality (Scarpetta et al., 2018). For critical systems this condition is enforced because avalanche profiles follows the same shape for all durations which means that this prediction is believed to be more precise than for non-critical systems and the exact values are important (Sethna et al., 2001; Nishimori and Ortiz, 2011). We found that average avalanche size scaled with duration 〈*A*〉(*D*) ∝ *D*^{γ} according to a power law and that the observed values of τ and β provided a good prediction γ = (β − 1)̸(τ − 1) of the fitted γ (Fig. 2*D*).

Specifically, of the 51 recording groups from the first 20 min period, the fitted scaling relation exponents were γ_{f} = 1.19 ± 0.05 and the predicted scaling relation exponents were γ_{p} = 1.17 ± 0.35. For the additional second 20 min period (18 groups/neurons), the fitted scaling relation exponents were γ_{f} = 1.21 ± 0.05 and the predicted scaling relation exponents were γ_{p} = 1.28 ± 0.21.

To affect a more convincing analysis, we defined four stringent criteria that must be independently satisfied before any set of avalanches can be deemed consistent with network dynamics near a critical point (see “Experimental design and statistical analysis” section). Overall, of the 69 groups of recordings (which includes 18 out of 51 cells twice), 98.6% had power laws in both the size and duration distributions of avalanches and 92.8% had scaling relations that were well fit by power laws (*R*^{2} > 0.95). All were deemed nontrivial by the test (γ* _{f}* − 1) > σ

_{γf}where σ

_{γf}is the dataset SD and σ

_{γf}= 0.051. The smallest value was γ

_{f}= 1.094. The fourth constraint, that the

*R*of the predicted scaling relation was within 10% of the best-fit scaling relation, was satisfied 85.6% of the time. Together, this set of criteria cannot measure distance from a critical point nor eliminate false positives. However, the take away is that 81% of all recording groups examined were judged to be consistent with network activity near a critical point.

^{2}Separating out results: 76% of the 51 recording groups from the first 20 min period and 94% of the recording groups from the second 20 min period were judged consistent with criticality. The general pattern is that the first 20 min period and the second are both consistent with criticality, but the second group meets our criteria much more frequently. This could be an effect related to the length of time we are able to maintain a patch or it could be that a better patching results in both longer stable recording ability and better inference of dynamical network properties.

To further discount the possibility of false positives we investigated whether the avalanches within our dataset exhibited “shape collapse” (Fig. 2*E*). The scaling relation is a consequence of self-similarity (Sethna et al., 2001; Papanikolaou et al., 2011; Friedman et al., 2012; Marshall et al., 2016; Shaukat and Thivierge, 2016; Cocchi et al., 2017). In other words, avalanches all have the same “hump shape” no matter how long they last; this shape is called the scaling function or avalanche profile. The shape collapse also provides an independent estimate of the scaling relation exponent γ. If estimated exponent, γ_{SC}, matches the fitted exponent, γ_{f}, then it is considered strong evidence of critical point behavior. For critical systems, the average avalanche profile of an avalanche of duration *D* is given as *A*(*t*,*D*) = *D*^{(γ−1)} where *D*^{(γ−1)} is a coefficient governing the scaling of height with duration, and is the scaling-function that describes the universal shape of an avalanche at any duration. The similarity of avalanche profiles of different durations is qualitatively judged (Sethna et al., 2001; Beggs and Plenz, 2003; Friedman et al., 2012; Pruessner, 2012; Timme et al., 2016) by plotting empirically estimated scaling functions for several durations on top of one another after they have been rescaled as part of the process of estimating γ_{SC}.

We obtained shape collapse across more than one order of magnitude (between ∼50 and 700 ms) of avalanche durations. Below 50 ms, distinct peaks arose. Above 700 ms, the profile height grew faster than the power-law scaling that worked for shorter duration avalanches. This is observed as an apparent outlier in Figure 2*E*. This likely marks the point where avalanches become so long and so large that they begin to weakly activate the nonlinear action potential mechanism of the neuron. When comparing to plausible alternatives to *V*_{m} in later sections, we included analysis of mean curvature and avalanche profile peak height along with visual inspection of shape collapse quality (Fig. 2*E*). The shape collapse plots begin with short avalanches (20 ms) that are below the median lower cutoff for power-law behavior (which was 256 ms) but are well predicted by the scaling relation.

The exponents estimated from the shape collapse were a good match for both the predicted and fitted scaling relation exponents. The groups of recordings from the first 20 min yielded γ_{SC} = 1.1868 ± 0.042. The average matched absolute percentage error was 1.3% with respect to γ_{f}. A matched signed-rank difference of median test revealed that γ_{f} was not significantly different from γ_{SC}, simple difference effect size *r*_{SDF} = 0.089, *p* = 0.063 (where *p* < 0.05 indicates that they are different).

This stage of the analysis showed that, when fluctuations of *V*_{m} are treated like neuronal avalanches, they are consistent with criticality by the standards of power laws governing size and duration. We also showed that *V*_{m} avalanches exhibit geometrical self-similarity across more than one order of magnitude. These factors showed that the cortical circuits driving fluctuations of *V*_{m} are consistent with activity near a critical point according to standards of self-similarity. In our next investigation, we compared with population data from microelectrode arrays and other results from the literature to determine whether *V*_{m} fluctuations are consistent with the universality requirement of behavior near critical points and if they can be used to measure dynamical network properties.

*V*_{m} fluctuations are consistent with avalanches from previously obtained microelectrode array LFP recordings

We sought to interpret our results from the analysis of single-neuron *V*_{m} fluctuations in the context of the more commonly used analysis of multiunit spiking activity (Friedman et al., 2012; Shew et al., 2015; Marshall et al., 2016; Karimipanah et al., 2017a) or multisite LFP event detection from MEA data (also known as “multielectrode array”) (Beggs and Plenz, 2003; Shew et al., 2015).

In a previous study, avalanche analysis was performed on LFP multisite MEA recordings from the visual cortex of a different set of 13 *ex vivo* eye-attached whole-brain preparations in turtles (Shew et al., 2015). Avalanches were inferred from the steady-state (after on response transients but before off response transients) of responses to visual presentation of naturalistic movies as opposed to the resting-state activity between presentations (which is where the *V*_{m} data come from). Avalanche size and duration distributions followed power laws.

The median exponents were τ = 1.94 ± 0.27 for the avalanche size distributions and β = 2.14 ± 0.32 for the avalanche duration distributions (Fig. 3*A*). A scaling relation existed with average exponent γ_{f} = 1.20 ± 0.06 fitted to the data and γ_{p} = 1.19 ± 0.07 from the average of the predicted scaling based on theory. The scaling power-law extended over one to two orders of magnitude. Critical branching was more firmly established in Shew et al. (2015) by analyzing the branching ratio. The branching ratio is the average ratio of events (i.e., spikes) from one moment in time to the next, but only during identified avalanches. A critical branching network has a branching ratio of one, but empirically estimating it requires discrete events and an assiduous choice of time binning for analysis. Shew et al. (2015) found that a branching ratio near one that was robust to reasonable choices of time bin and varied with choice of time bin in expectation with critical branching. We are not aware of methods for estimating a branching ratio in continuous signals like *V*_{m}.

The set of avalanche size, duration, and scaling relation exponents obtained from *V*_{m} fluctuations (Fig. 3*B*) were not distinguishable from the MEA obtained set. The fitted scaling relation exponent γ_{f} had the least variability of all three kinds of exponents, so it is the most likely to show a difference. Thus, if a difference is not significant, then this suggests universality more strongly than for the avalanche size τ or duration β distribution exponents.

When we limited our analysis to the first 20 min period that contained multiple recordings (51 cells), neither the fitted scaling relation exponent nor the predicted scaling relation exponent were significantly different from the MEA results. The Wilcoxon rank-sum difference of medians test against the MEA data yielded (*r*_{SDF} = 0.164, *p* = 0.37), and (*r*_{SDF} = 0.08, *p* = 0.67), respectively. The median exponent values for the size and duration distributions were not significantly different from the median of the MEA data (*r*_{SDF} = 0.164, *p* = 0.37) and (*r*_{SDF} = 204, *p* = 0.265), respectively.

These results establish *V*_{m} fluctuations as an informative gauge of high-dimensional information while also demonstrating that the power-law characteristics are universal properties of the brain by showing a close match between data at different scales and under different conditions. Further underscoring universality, our results are also similar to the critical exponents measured from other animals such as the τ = 1.8 result from *in vivo* anesthetized cats (Hahn et al., 2010). Although an exhaustive literature search was not conducted here, others have conducted incomplete surveys (Ribeiro et al., 2010; Priesemann et al., 2014).

### Single-neuron estimate of network dynamics is optimized at the network critical point

To gain a deeper insight into the relation between single-neuron input and network activity, we investigated a model network of probabilistic integrate and fire model neurons (Kinouchi and Copelli, 2006; Larremore et al., 2011a,b, 2012, 2014; Karimipanah et al., 2017a,b). This model network contains fundamental features of cortical populations, such as low connectivity, inhibition, and spiking while being sufficiently tractable for mathematical analysis (see “Model simulations” section).

In brief, the model network consists of *N* = 10^{4} binary probabilistic model neurons (Fig. 4*A*). The connection probability *c* results in a mean in-degree and out-degree of *cN*. The connection strength from neuron *j* to neuron *i* is quantified in terms of the network adjacency matrix *W*. Each connection strength *W*_{ij} is drawn from a distribution of (initially) positive numbers with mean η, where the distribution is uniform on [0,2η]. A fraction χ of the neurons are designated as inhibitory; that is, their outgoing connections are made negative. The binary state *S*_{i}(*t*) of neuron *i* is updated according to *S _{i}*(

*t*) = Θ(∑

*(*

_{j}^{N}W_{ij}S_{j}*t*− 1) − ξ

*(*

_{i}*t*)) where ξ

_{i}(

*t*) is a random number between 0 and 1 drawn from a uniform distribution and Θ is the Heaviside step function.

The largest eigenvalue, λ = η*cN*(1 − 2χ), of the network adjacency matrix *W* characterizes the network dynamics, with critical network dynamics occurring at λ = 1. This tuning parameter λ controls the degree to which spike propagation “branches”: λ = 1 means that one spike creates one other spike on average, λ > 1 implies that one spike creates more than one other spike, whereas λ < 1 means that one spike creates less than one other spike (Haldeman and Beggs, 2005; Kinouchi and Copelli, 2006; Levina et al., 2007; Larremore et al., 2011b, 2012; Kello, 2013).

The input to model neuron *i* is *P _{i}*(

*t*) = ∑

*(*

_{j}^{N}W_{ij}S_{j}*t*− 1) and provides the link between network activity and single-neuron activity. From this we can derive a simple mathematical result characterizing how estimation of network properties is optimized at criticality.

If we let *K _{i}*(

*t*− 1) denote the number of active neurons in the presynaptic population of neuron

*i*, then we can rewrite the input to a model neuron as a sum of independent and identically distributed random variables drawn from the nonzero entries of

*W*:

*P*(

_{i}*t*) = ∑

_{k}

^{Ki(t − 1)}W

*. After implementing inhibition by inverting some elements of*

_{ijk}*W*, the distribution of weights is not uniform but piecewise uniform. Weights are drawn uniformly from the interval [−2η,0] with probability χ and from the interval [0,2η] with probability 1 − χ. The mean of the nonzero entries of

*W*are denoted with a prime so that the mean is 〈

*W*〉 = η(1 − 2χ) and the SD is . Now we can find the mean behavior of the input integration function as it relates to the presynaptic population: We learn three things by examining the mean behavior of the input integration function. First, the mean grows as

_{ij}^{′}*O*(

*K*) but the SD grows as the root O(), so the function becomes a more precise estimator of network activity with increasing activity in the presynaptic population (increasing

_{i}*K*

_{i}). Second, the input integration function

*P*

_{i}(

*t*) is rarely negative. At the parameter combination

*c*= 0.005 and χ = 0.25 (which has the largest variance relative to the mean), the mean becomes >1 SD larger than 0 when

*K*

_{i}> 5. Third, and most importantly, the input integration function is an averaging operator and the tuning parameter λ biases that averaging operation. To show this, we only need two observations: the instantaneous firing rate averaged over the presynaptic population is the number of active neurons divided by the expected total number of presynaptic neurons, ω

*(*

_{i}*t*) =

*K*(

_{i}*t*)/

*cN*. Next, we rearrange the definition of λ to get λ̸

*cN*= η(1 − 2χ). Substituting these two observations into the mean behavior of our input integration function we get the key mathematical result: Note that

*P*(

_{i}^{′}*t*) =

*P*(

_{i}*t*) × Θ(

*P*(

_{i}*t*)) directly represents the probability for the neuron to spike at time

*t*.

These results demonstrate that the inputs to a neuron *P*_{i} and the instantaneous firing rate of that neuron are the result of an averaging operator acting on the presynaptic population, which is a subsample of the network. Furthermore, the tuning parameter λ not only modulates the relationship of single neuron firing to downstream events (also known as branching), but also governs how the input to a neuron relates to the presynaptic population. It biases the averaging operator to either amplify firing rate (λ > 1) or dampen it (λ < 1). Therefore, our model implements both critical branching and the inverse of the critical branching condition, a critical coarse-graining condition. The model is a network of subsampling operators that only capture whole-system statistics when λ = 1 and the operators reflect an unbiased stochastic estimate of mean firing rate among the subsample (the presynaptic population). This averaging operation may exist in many kinds of networks, including those with structure and those that are not critical branching networks, so this result helps to establish plausible generalizability.

To further evaluate the relation between single-neuron input and network activity under different conditions, we simulated the described network of 10^{4} model neurons for a total of 405 different parameter combinations, including connection probability, inhibition, and maximum eigenvalue (Fig. 4*A*), each parameter combination was repeated 10 times. We then compared the avalanche analysis results of simulated network activity *F*(*t*) = ∑_{i=1}^{N}*S _{i}*(

*t*) with the input to a single neuron (the input integration function). However,

*P*(

_{i}*t*) = ∑

*(*

_{j}^{N}W_{ij}S_{j}*t*) is the probability that neuron

*i*will fire at time

*t*, also known as the instantaneous firing rate of neuron

*i*.

*V*_{m} is not a direct representation of firing rate, but rather the firing rate is related to synaptic input through the *F–I* curve which is nonlinearly related to *V*_{m}. This nonlinearity could destroy the correspondence between the simulated single neuron signal and network activity. To better facilitate comparison of the simulated input integration function with the experimentally recorded *V*_{m}, we constructed a proxy for the subthreshold *V*_{m}, Φ_{i}(*t*), of a model neuron by convolving the simulated input *P*_{i}(*t*) with an alpha-function (see “Model Simulations” section).

The parameter space has four distinct patterns of critical network behavior (Fig. 4*B*). Qualitatively, these were reflected in the network activity. The presence of these paradoxical behaviors may indicate the presence of second phase transition tuned by the balance of excitation to inhibition (Shew et al., 2011; Poil et al., 2012; Kello, 2013; Hesse and Gross, 2014; Larremore et al., 2014; Scarpetta et al., 2018). Several key results differ strongly and thus are reported separately for these regions of parameter space.

These regions are defined in terms of the connection density and inhibition and are shown in Figure 4*B*. First is the “positive weights” region, where there is no inhibition (χ = 0) and the network is a standard critical branching network. The second region, “quiet,” has a small increase in the fraction of inhibitory neurons. Activity lasts slightly longer than for the classically critical network. The third region is called the “switching” regime because network activity switches between a low mean and a high mean (like “up and down states”; Destexhe et al., 2003; Millman et al., 2010; Larremore et al., 2014; Scarpetta et al., 2018). This occurred in the middle portion of the values of connectivity and inhibition. Last, we have the “ceaseless” region, with a large fraction of inhibition relative to connection density, where activity never dies out. This region is defined by *c* < (10*e*^{12χ} − 13)/100 and χ > 0. Three of these regimes are displayed in Figure 5*A*; the “quiet” region is mostly redundant to the “positive weights” region. The “ceaseless” and “switching” regimes exhibit sustained self-generated activity and are included with the intention to model ongoing spontaneous activity dynamics without contamination by externally imposed firing patterns (Mao et al., 2001).

We looked at the magnitude of relative error between estimated exponents for the avalanche size distribution (Fig. 4*C*) to determine how well our proxy neural inputs, φ_{i}(*t*), reflected network activity, *F*(*t*), in different parameter regions, and with different values for the tuning parameter, λ. Importantly, the least error occurred for λ = 1 with and without the presence of inhibitory nodes. This insensitivity to parameter differences supports the claim (Larremore et al., 2014) that the system becomes critical when λ = 1 even in the presence of inhibition.

However, the four regions of parameter space perform differently according to our four standardized criteria for consistency with criticality. In the “positive weights” region, 90% of 90 trials (nine points in parameter space with 10 trials per point) have network activity that meets all four criteria when the tuning parameter is set at criticality (λ = 1) (Fig. 4*C*). A total of 39% meet the criteria in the “ceaseless” region, 19% in the “quiet” region, and 67% in the “switching” region, which may indicate the location of a second phase transition and shows that evidence for precise criticality in this model is limited once inhibition is included.

As we vary the tuning parameter, we can clearly distinguish critical from noncritical systems. Overall, 47% of trials met all four criteria when λ = 1, whereas 3% did when λ = 0.95, 18% when λ = 1.015, 1% when λ = 0.9, and 1% when λ = 1.03 (Fig. 4*D*).

The estimated power-law exponents show that the avalanche size distributions for *F*(*t*), *P _{i}*(

*t*), and Φ

*(*

_{i}*t*) are most alike at criticality. Note that estimated exponents serves as the “scaling index,” a measure of the heavy tail even when a power-law is not the statistical model that fits best (Jeżewski, 2004). The fact that matching between network activity and the input integration function was best at criticality is important because it underscores the scale-free nature of critical phenomena and contrasts with the results obtained when testing subsampling methods with a different relationship to network structure (Priesemann et al., 2009; Yu et al., 2014; Levina and Priesemann, 2017).

While the system was both critical (λ = 1) and in the positive weights region, our *V*_{m} proxy Φ* _{i}*(

*t*) met all four criteria for consistency with criticality 74% of the time for 90 trials (Fig. 4

*D*), whereas

*P*(

_{i}*t*) met all four only 1% of the time. The network activity had avalanche size and duration exponent values τ

_{F}= 1.43 ± 0.04 and β

_{F}= 1.87 ± 0.09 (Fig. 4

*D*) and had a fitted scaling relation exponent, γ

*= 1.83 ± 0.02 and a predicted exponent γ*

_{Ff}*= 1.99 ± 0.23. The*

_{Fp}*V*

_{m}proxy, Φ

*(*

_{i}*t*), had slightly lower avalanche size and duration exponent values that fluctuated around the paired network values, τ

_{Φ}= 1.40 ± 0.06 and β

_{Φ}= 1.73 ± 0.17 (Fig. 4

*D*), and exclusively lower scaling relation exponents γ

_{Φf}= 1.68 ± 0.02. Although the unsmoothed

*P*(

_{i}*t*) varied considerably more it had size and duration exponents that were almost exclusively higher than the paired network values, τ

*= 1.87 ± 0.50 and β*

_{P}*= 1.87 ± 0.50 with a fitted scaling relation exponent that was exclusively lower, γ*

_{P}*= 1.68 ± 0.02.*

_{Pf}In Figure 5, we compared different population dynamics estimation techniques by looking at avalanches inferred from *P _{i}*(

*t*) (the inputs to neuron

*i*), and the

*V*

_{m}proxy Φ

*(*

_{i}*t*). Both

*P*(

_{i}*t*) and Φ

*(*

_{i}*t*) fluctuate about

*F*(

*t*), but

*P*(

_{i}*t*) is much noisier (Fig. 5

*A*); in the ceaseless regime,

*P*(

_{i}*t*) and Φ

*(*

_{i}*t*) are systematically offset. Avalanches inferred from Φ

*(*

_{i}*t*) had average sizes that scaled with duration (Fig. 5

*B*). Avalanches from Φ

*(*

_{i}*t*) consistently had duration and size distribution exponents that were closer to network avalanches than avalanches from

*P*(

_{i}*t*). However,

*P*(

_{i}*t*) performed satisfactorily in the sense that its error was systematically offset and best at criticality (Fig. 5

*C*).

Including inhibition introduced several important differences. For the ceaseless region with λ = 1, far fewer trails met our criteria; however, *P _{i}*(

*t*) followed

*F*(

*t*) much more closely. The network activity had avalanche size and duration exponent values τ

*= 1.48 ± 0.09, and β*

_{F}*= 1.53 ± 0.09 and had a fitted scaling relation exponent, γ*

_{F}*= 1.23 ± 0.11. The*

_{Ff}*V*

_{m}proxy, Φ

*(*

_{i}*t*), had slightly higher avalanche size and duration exponent values that fluctuated around the paired network values, τ

_{Φ}= 1.51 ± 0.19 and β

_{Φ}= 1.57 ± 0.17, but nearly identical scaling relation exponents γ

_{Φf}= 1.23 ± 0.11. Although the unsmoothed

*P*(

_{i}*t*) varied considerably more, it had size and duration exponents that were almost exclusively higher than the paired network values, τ

*= 1.88 ± 0.20 and β*

_{P}*= 2.18 ± 0.34, with a fitted scaling relation exponent that was slightly lower, γ*

_{P}*= 1.19 ± 0.07.*

_{Pf}When λ ≠ 1, both Φ* _{i}*(

*t*) and

*P*(

_{i}*t*) failed to meet all four criteria for criticality at the same high rate as

*F*(

*t*) (to within 1%). This lack of false positives confirms that these signals are useful for characterizing critical branching. In Figure 4

*B*, we calculated the absolute magnitude of relative error between the size exponent from avalanche analysis performed on

*F*(

*t*) and Φ

*(*

_{i}*t*). As expected, the avalanches were usually not power laws according to our standards; in this case, the exponent is known as the “scaling index” and describes the decay of the distribution's heavy tail (Jeżewski, 2004).

When we set λ = 0.95, we see a moderate deterioration in the ability of either Φ* _{i}*(

*t*) or

*P*(

_{i}*t*) to recapitulate network exponent values. The error is no longer systematic; so they cannot be used to predict network values. The variability of the exponents increases greatly for Φ

*(*

_{i}*t*), whereas it decreases for

*P*(

_{i}*t*). The exponent error increases slightly over the λ = 1 case and the base of the distribution is much broader.

Reducing λ further, to λ = 0.90, the input integration function, *P _{i}*(

*t*) ∼ λω

*(*

_{i}*t*− 1), rapidly dampens impulses (ω

_{i}is the instantaneous firing rate over the presynaptic population for neuron

*i*). Variability continues to increase and a systematic offset does not return. Exponent error is now much broader. With branching this low, events often are not able to propagate to the randomly selected neuron; an exception is the “ceaseless” regime where activity is still long lived.

When we set λ = 1.015, we see a dramatic deterioration in the ability of either Φ* _{i}*(

*t*) or

*P*(

_{i}*t*) to recapitulate network values. Variability in exponent estimation increases for both Φ

*(*

_{i}*t*) and

*P*(

_{i}*t*). Exponent error increases rapidly, underscoring the inability to estimate network activity from neuron inputs.

Increasing λ further to λ = 1.03 produces an input integration function, *P _{i}*(

*t*) ∼ λω

*(*

_{i}*t*− 1), which rapidly amplifies all impulses and the network saturates. The effect is that variability in the estimated exponents decreases and a systematic offset returns, with both Φ

*(*

_{i}*t*) and

*P*(

_{i}*t*), producing exponents that are exclusively and considerably higher than network values. Exponent error reveals that estimating network properties from the inputs to a neuron is probably not possible for supercriticality in this model.

The results here show that the *V*_{m} proxy represents an effective way of subsampling network flow. This is a hallmark of the near-critical region in the PIF model and a manifestation of scale-freeness. Criticality in our model corresponds to the point when the inputs to a neuron represent an average of the activity of the presynaptic population. Importantly we explored why it works, as well as showing that it does work in experimental data. This analysis, presented in forthcoming sections, uncovered that proper temporal and spatial aggregation is important as is the role of inhibition in *V*_{m} dynamics. This supports both the criticality hypothesis and tight balance (Barrett et al., 2013; Boerlin et al., 2013; Denève and Machens, 2016). Additionally, it has specific implications for the information content of *V*_{m}.

### Predicted scaling relation exponent is more stable than avalanche size or duration exponents

A key part of the study of criticality in neural systems is the assumption that biological systems must self-organize to a critical point. The precise critical point is a very small target for a self-organizing mechanism in any natural system. Therefore, a key question is whether the self-organizing mechanism of the brain prioritizes efficiently achieving information processing advantages of scale-free covariance at the expense of being slightly sub or supercritical (which is a larger target) (Priesemann et al., 2014; Tomen et al., 2014; Williams-García et al., 2014; Gautam et al., 2015; Clawson et al., 2017).

Our data offered unexpected insight. It is known that so long as three requirements are met the scaling relation will be marginally obeyed: Avalanche size and durations must be power-law distributed (with exponents τ and β, respectively) and average size must scale with duration according to a power-law with exponent γ. Given those three requirements one can derive a prediction for the scaling exponent, γ* _{p}* = (β − 1)/(τ − 1) without needing to assume criticality (Scarpetta et al., 2018). However, without any other assumptions, one expects β and τ to be independent, so plotting one against the other should make a point-cloud that is symmetrical, not stretched along a trendline (Fig. 3).

We analyzed the independence of τ, β, and γ measured from experimental data (where self-organization is hypothesized) and compared it with model data (where self-organization is impossible, but criticality is guaranteed). We found that β and τ are more independent and the predicted scaling relation is more variable for the model than for experimental data in which β and τ covary, apparently to maintain a fixed scaling relation prediction.

The previous multisite LFP recordings displayed a range of values for the avalanche size τ and duration β distribution exponents across the tested brain preparations. Interestingly, the exponent values were not independent; rather, the duration exponent varied linearly with the size exponent (Shew et al., 2015) (Fig. 3*A*). The single-neuron *V*_{m} fluctuations reported here produced a similar linear relationship between size and duration exponent (Fig. 3*B*). Algebraic manipulation of the predicted scaling exponent γ* _{p}* = (β − 1)/(τ − 1) provides a clue. If the scaling relation (β − 1) = γ(τ − 1) is obeyed and if γ

_{p}is a fixed universal property, then the linear relationship β

*∝ 〈γ*

_{j}*〉τ*

_{p}*holds across different cells and animals.*

_{j}To demonstrate this important result, variability in the predicted scaling relation is much less than expected, we propagate errors and assume independent β and τ. We would expect the SD of γ_{p} to be σ_{γp}^{*} ∼ ∼ 0.72, which is approximately twice the real value in *V*_{m} data, σ_{γp} ∼ 0.35.

The Pearson correlation, ρ, confirms strong dependence between τ and β, ρ_{τβ} = 0.61, *p* = 2.57 × 10^{−6} for the *V*_{m} data, whereas for the MEA data ρ_{τβ} = 0.96, *p* =1.01 × 10^{−7}. From this, we confirm what Figure 3 shows: the variability in τ and β are not independent and this implies the existence of an organizing principle connecting τ to β. Whatever the principle may turn out to be, one of its effects is the maintenance of low variability in γ_{p} at the expense of greater variability in τ and β.

A principle reason to suspect self-organization is that this trend is not seen in the model results. Importantly, τ and β are independent of the scaling-relation exponent function, although still weakly correlated. In this model, there is no adaptive organizing principle driving this network to criticality, instead the structure is fixed and set to be at the critical point. This shows how systems behave in the absence of self-organization. No parameter is being maintained at low variability at the expense of other parameters.

Limiting ourselves to simulated network activity for the λ = 1 case without inhibition (Fig. 4*C*), propagation of errors leads us to expect the SD of the scaling-relation prediction to be σ*_{γFp} ∼ 0.27, which is very close to real value of σ_{γFp} ∼ 0.23. The correlation is statistically significant at the 5% level, but much smaller ρ_{τβ} = 0.23, *p* = 0.027. This was noted in the original study (Shew et al., 2015), where the authors were able to reproduce the linear trend between avalanche size and duration exponents by simulating a network with synaptic depression to adaptively restore critical behavior after an increase in network drive. They showed that the trendline is produced by corrupting their simulated data via randomly deleting 70–90% of spiking events and then changing the way that they group events in time (adaptive time binning). Our *V*_{m} fluctuations have no counterpart to the adaptive time binning other than the intrinsic membrane time constant, which is not manipulated experimentally.

In conclusion, the linear trend between avalanche size and duration exponents is not a universal property of critical systems because it was not found in the model. This suggests that the linear trend is enforced by an organizing principle at work in the brain but absent in the model. This principle prioritizes maintaining stability in either the scaling of avalanche size with duration, or the power-law scaling of autocorrelation which is closely related to the scaling relation and scale-free covariance via the power-law governing autocorrelation (Bak et al., 1987; Sethna et al., 2001).

### Nonlinearity and temporal characteristics such as high-order correlation, proper combination of synaptic events, and signal timescale are required to reproduce network measures from single-electrode recordings

To demonstrate that subthreshold *V*_{m} fluctuations can be used as an informative gauge of cortical population activity, it is necessary to compare against alternative signals that have either been used by experimentalists as a measure of population activity or that share some key features of *V*_{m} but are missing others. By making these comparisons, we can illuminate which features of the *V*_{m} signal are responsible for its ability to preserve properties of cortical network activity. Additionally, it is necessary to determine whether the statistical properties of avalanches can be explained by random processes unrelated to criticality. To address these points of the investigation, we analyzed five surrogate signals: single-site LFP recorded concurrently with the *V*_{m} recordings, two phase-shuffled versions of *V*_{m} recordings, computationally inferred excitatory current, and the same inferred excitatory current further transformed to match *V*_{m} autocorrelation (which tests the role of IPSPs by making a *V*_{m}-like signal that lacks them).

### Negative fluctuations of LFP disagree with *V*_{m} and MEA results and are inconsistent with avalanches in critical systems

The first alternative hypothesis to test is whether the LFP could yield the same results. We used low-pass filtered and inverted single site LFP, which is commonly believed to measure local population activity. However, in our analysis, it did not recapitulate the results from either MEA or *V*_{m} avalanche analysis. We obtained viable single-site LFP recordings (see “Extracellular recordings” section), simultaneous and adjacent with whole-cell recordings, for 38 of the 51 neurons reported above. We performed avalanche analyses on the LFP recordings using a procedure like the one described for the *V*_{m} recordings (see “Intracellular recordings” section) (Fig. 6). LFP recordings were grouped in the same way that *V*_{m} recordings were to match them for comparison. However, the numbers of recordings are not the same because there were two or three cells being patched alongside (within 300 μm) one extracellular electrode and there was not always a simultaneous LFP recording. LFP also produced more avalanches per 2–5 min recording, *N _{AV}* = 1128 ± 348. There are 23 20 min periods spanning multiple LFP recordings. These recordings were gathered into groups and matched against 49

*V*

_{m}recording groups (38 from the first 20 min period, 11 from the second). Additionally, there were 16 20 min periods spanning only one LFP recording but with >500 avalanches. The concurrent

*V*

_{m}recordings did not have enough avalanches. This gives us 39 LFP avalanche datasets.

The LFP recording groups performed poorly according to our four criteria for consistency with criticality. Of the 39 LFP recording groups, only 41% had acceptable scaling relation predictions and only 36% met all four standard criteria for criticality (Fig. 7*A*). The additional criterion of shape collapse was not observed (Fig. 6*C*); there was no linear trend among the exponents governed by the scaling relation and the exponents did not match MEA data (Fig. 3*A*). However, 85% produced power-law fits for size and duration, 92% had scaling relations well fit by power laws, and all were nontrivial. We expect from previous data (Touboul and Destexhe, 2017) that some fraction of noncritical data will pass the four standard criteria by chance as long as the data have a 1/*f* power spectrum.

To emphasize that these results are chance, we can limit ourselves to just those with the best chance of meeting the scaling relation criteria by picking those that have power laws in the size and duration distributions. This is enough to expect the scaling relation to be obeyed if mean size scales geometrically with duration (Scarpetta et al., 2018). It is still the case that only 42% of recording groups meet the three remaining standard criteria for consistency with criticality. Therefore, having power laws is statistically independent of meeting the other criterion for consistency with criticality.

Not only does the single-site LFP data differ from MEA and *V*_{m} data because it fails to demonstrate consistency with criticality, it is also the case that the scale-free properties that do exist are not representative of the MEA data or the simultaneous *V*_{m} recordings. The failure was not because LFP recordings cooccurred with decreased consistency with criticality more generally. Eighty-one percent of the matched *V*_{m} recordings met all the criteria, whereas 58% of the LFP recordings did, a statistically significant dissimilarity (*r _{OR}* = 7.65,

*p*= 1.1 × 10

^{−5}).

The estimated exponents from all 39 LFP recording groups were highly variable. The duration distribution and scaling relation were most dissimilar to *V*_{m} and MEA data. Of the 33 LFP groups that were power-law distributed, the avalanche size exponent had a median value of τ = 1.90 ± 0.63, whereas the duration exponent was β = 1.41 ± 0.9 (very low) (Fig. 7*A*) and the fitted exponent was γ* _{f}* = 1.11 ± 0.02. The predicted scaling-relation exponents were inaccurate, with γ

*= 0.89 ± 0.76 for the subset of recording groups that had power laws.*

_{p}The extreme variability makes it difficult to determine whether the size and duration exponents match other data, but the fitted scaling relation exponent was much less variable and more clearly separated from MEA or *V*_{m} results. The matched difference of median test (Wilcoxon signed-rank) between 49 recording groups found that the best fit τ, (τ = 1.90 ± 0.63), was not significantly distinguishable from the *V*_{m} data (*r _{SDF}* = 0.15,

*p*= 0.33), but β, (β = 1.41 ± 0.9), was dissimilar with a comparable effect size (

*r*= 0.17,

_{SDF}*p*= 0.028); γ

_{f}, (γ

*= 1.11 ± 0.02), was also dissimilar (*

_{f}*r*= 0.25,

_{SDF}*p*= 7.1 × 10

^{−15}).

When comparing against the 13 samples of MEA data γ_{f} was significantly different from the MEA data (*r _{SDF}* = 0.88, and

*p*= 9.21 × 10

^{−08}). This contrasts with our comparison between

*V*

_{m}and MEA data. In that case the scaling relation was not distinguishable even with 51 points of comparison and very low variability making a difference easier to detect. However, because of their extreme variability the size and duration exponents fail a 5% significance threshold for distinguishing from the MEA data by a Wilcoxon rank-sum result (

*r*= 0.06,

_{SDF}*p*= 0.766 for τ and

*r*= 0.29,

_{SDF}*p*= 0.123 for β). This failure of inverted LFP to show the same statistical properties as multiunit activity may add a caveat to the assumptions behind the use of inverted LFP as a proxy for population activity (Kelly et al., 2010; Einevoll et al., 2013; Okun et al., 2015). Specifically, the amplitude of single-electrode negative LFP excursions is ambiguously related to the number of spiking neurons, whereas the use of electrode arrays as described previously (Beggs and Plenz, 2003; Shew et al., 2015) is more appropriate.

To summarize, single-site LFP fluctuations result from the superposition of local spiking and extracellular synaptic current from juxtaposed network elements (Kajikawa and Schroeder, 2011; Einevoll et al., 2013; Pettersen et al., 2014; Ness et al., 2016). These fluctuations were found to be less informative about the network dynamics than single-neuron *V*_{m} fluctuations. *V*_{m} fluctuations result from the superposition of EPSPs and IPSPs, indicating neuronal responses propagating in a manner consistent with the true neural network architecture. In other words, synaptic and spiking events driving fluctuations at single extracellular electrodes may be too badly out of sequence and distorted to faithfully represent neuronal avalanches, whereas the sequence of synaptic and spiking events driving somatic *V*_{m} fluctuations is functionally relevant by definition.

### Stochastic surrogates are distinguishable from *V*_{m} or MEA results, revealing the importance of nonlinear filtering

After eliminating inverted LFP as an alternative single-electrode signal, it was important to establish whether our results could have been created from a linear combination of independent random processes (Touboul and Destexhe, 2017; Priesemann and Shriki, 2018) similar to those used when contesting evidence for critical brain dynamics (Bédard et al., 2006; Touboul and Destexhe, 2010, 2017). We also wanted to learn what effects nonlinearity (non-Gaussianity) has in signals such as *V*_{m}.

To address these questions, we used both the AAFT and UFT phase shuffling algorithms (see “Experimental design and statistical analysis” section). AAFT (Fig. 6) preserves both the exact power spectrum (autocorrelation) of the signal and nonlinear skew of signal values, but randomizes the phase (higher-order temporal correlations). UFT is the same but forces the distribution of signal values to be Gaussian. Using both allows us to attribute some characteristics to nonlinear rescaling and others to precise temporal correlation structure.

Phase shuffling tends to preserve power laws since it explicitly preserves the 1/*f* trend of the power spectrum. However, the matched signed-rank test reveals that the values of the exponents change in both methods. Under UFT transformation the scaling relation and shape collapse became more trivial and like the LFP. This suggests that both the nonlinear rescaling of input currents by membrane properties and the way that input populations interact throughout the intricate dendritic arborization are important.

For the 51 recording groups from the first 20 min, the AAFT reshuffled data yield a median size exponent of τ = 1.74 ± 0.29, whereas the duration exponent was β = 2.0 ± 0.34 (Fig. 7*D*). The fitted scaling relation exponent was γ* _{f}* = 1.19 ± 0.06 and the predicted scaling relation exponent was γ

*= 1.21 ± 0.49.*

_{p}Pairing the surrogates to the original *V*_{m} data and performing the Wilcoxon signed-rank test for difference of medians gives (*r _{SDF}* = 0.053,

*p*= 2 × 10

^{−4}), (

*r*= 0.091,

_{SDF}*p*= 0.08), and (

*r*= 0.207,

_{SDF}*p*= 3 × 10

^{−5}) for τ, β, and γ

_{f}, respectively. Therefore, τ and γ

_{f}are both significantly different; this is supported by the fact that only 55% of the groups meet all four standard criteria for criticality, whereas 76% of meet them for the original

*V*

_{m}time series. This difference between success rates is significant by Fisher's exact test (r

_{OR}= 2.67,

*p*= 0.0363).

The failure mode for AAFT shuffled data was almost entirely in reduced goodness of fit (*R*^{2}) for a power-law fit to its scaling relation, 17% fewer recording groups met the criterion *R*^{2} > 0.95 than for *V*_{m} (*r _{OR}* = 4.18,

*p*= 0.0093). When the shape collapse is examined, we see another clear, if qualitative, difference in the symmetry of any presumed scaling function (Fig. 6

*C*). The AAFT shuffled dataset is not consistent with critical point behavior. Thus, we show that the exponent values and evidence for criticality, especially scaling and shape collapse that we inferred from

*V*

_{m}are not likely to come from random processes and are dependent on nonlinear temporal correlation structure.

The key feature of the UFT result is that the fitted scaling relation exponent is much lower, γ* _{f}* = 1.05 ± 0.049, which is significantly less than for AAFT (

*r*= 0.25,

_{SDF}*p*= 1 × 10

^{−13}) and less than the LFP (

*r*= 0.228,

_{SDF}*p*= 3 × 10

^{−6}). It is very close to trivial scaling but is still distinguishable from γ

_{f}= 1 at a population level via the sign test (

*r*= 0.843,

_{SDF}*p*= 2 × 10

^{−10}). Because the fitted scaling relation exponent and shape collapse were similar in both the UFT and LFP data, this suggests that lack of nonlinear rescaling (nonlinear filtering) may be a key feature of LFP that explains its failure to accurately reflect critical point behavior.

The UFT was universally poorer performing, 39% do pass the criticality test; but given that the scaling relation exponent is so low, this is simply random chance and significantly worse than the *V*_{m} results (*r _{OR}* = 5.04,

*p*= 3 × 10

^{−4}). The UFT phase shuffling results obtain a median size exponent of τ = 1.69 ± 0.45, while the duration exponent was β = 1.81 ± 0.49. The predicted scaling relation exponent was γ

*= 1.01 ± 0.72. All are significantly different from the*

_{p}*V*

_{m}results, (

*r*= 0.183,

_{SDF}*p*= 0.005), (

*r*= 0.199,

_{SDF}*p*= 2 × 10

^{−4}), and (

*r*= 0.249,

_{SDF}*p*= 2 × 10

^{−13}) for τ, β, and γ

_{f}, respectively. These results are redundant with the AAFT, confirming that our results do not have a trivial explanation.

When the scaling relation was examined, we saw another clear, if qualitative, difference in the symmetry of any presumed scaling function (Fig. 6*C*). When taken together, our four standardized criteria followed by shape collapse analysis let us distinguish phase-shuffled *V*_{m} fluctuations from the original *V*_{m} fluctuations, even limiting ourselves to data that meets the four criteria. Therefore, the phase-shuffled data showed that the evidence for criticality in the original *V*_{m} fluctuations is dependent on nonlinear temporal correlations.

### Excitatory and inhibitory synaptic activity are both required for *V*_{m} fluctuations to match MEA avalanches

Having learned that single-site LFP recordings cannot be used to accurately infer the statistics of population activity and knowing that low-pass filtered and inverted LFP is believed to reflect excitatory synaptic activity (Kajikawa and Schroeder, 2011; Buzsáki et al., 2012; Einevoll et al., 2013; Ness et al., 2016), this begs the question: to what extent do excitatory synaptic events contain evidence for network criticality?

Somatic *V*_{m} fluctuations are the complex result of spatially and temporally distributed excitatory and inhibitory synaptic inputs further mangled by active and passive membrane properties in dendrites and soma. There is reason to believe that these features conspire to enforce the condition that *V*_{m} faithfully represents inputs to the presynaptic network (Barrett et al., 2013; Boerlin et al., 2013; Denève and Machens, 2016), similar to how input signals relate to presynaptic populations in our model. To address the stated question, we estimated the excitatory synaptic conductance changes *g*_{exc}^{*} from the *V*_{m} recordings using a previously developed inverse modeling algorithm (Yaşar et al., 2016) and applied the avalanche analysis on the inferred *g*_{exc}^{*} time series (Fig. 6).

The inferred excitatory conductance is plausibly related to the presynaptic population, however it failed to be a reliable measure of network dynamics (Fig. 7*B*). We cannot know whether the failure is because excitatory current does not contain enough information or because the signal's time constant is too short. Power laws in the avalanche size and duration distributions were observed in only 12% of the 51 groups from the first 20 min of recording. Comparing with *V*_{m}, this was very different (*r _{OR}* = 375,

*p*= 6 × 10

^{−14}). Shape collapse was absent from the inferred excitatory conductance (Fig. 6

*C*) and none passed all four criteria for criticality. From this, we conclude that inferred excitatory conductances are not a good network measure.

One of many potential reasons for this failure could be the much shorter time constant of the inferred *g*_{exc}^{*} signal compared with the *V*_{m} signal. We saw exactly that situation when examining model results: *P _{i}*(

*t*) failed to reproduce network values as well as its smoothed version, Φ

*(*

_{i}*t*). Therefore, we smoothed the

*g*

_{exc}

^{*}signal with an alpha-function chosen because it should impose a similar non-Gaussian distribution as the

*V*

_{m}signal. The time constant of the alpha-function was tuned to minimize the error between the autocorrelation of the smoothed

*g*

_{exc}

^{*}signal and the original

*V*

_{m}signal. By doing so, we create a signal with a 1/

*f*power spectrum that should exhibit power laws and reproduce many

*V*

_{m}statistical features (Fig. 6).

Reinstating the autocorrelation does not summon the return of scale-freeness. The smoothed signal did demonstrate power laws (94%) and one serendipitously met the standardized criteria for consistency with critical point behavior (Fig. 6*D*). However, this is chance. The average coefficient of determination for a fitted scaling relation on a log-log plot was *R*^{2} = 0.84 ± 0.14, so overall average avalanche sizes did not scale with duration as a power law. Nonetheless, this is a substantial improvement on the unsmoothed version *R*^{2} = 0.68 ± 0.17. This is a statistically significant difference (*r _{SDF}* = 0.054,

*p*= 3 × 10

^{−4}).

The smoothed inferred *g*_{exc}^{*} signal (Fig. 6*A*) is visually more like the original *V*_{m} (Fig. 2*B*) than the AAFT shuffled *V*_{m} surrogate (Fig. 6*A*); however, it was a worse match. This shows that signals dependent only on excitation, even ones with very similar non-Gaussian distribution and power-spectrum do not reflect the statistics of population activity. Interactions between EPSPs and IPSPs may be needed.

In conclusion, the single-site LFP, the phase-shuffled recorded *V*_{m}, and the inferred excitatory conductance *g*_{exc}^{*}, including its smoothed version, all failed to reveal the critical network dynamics. However, there are either similarities between the signals or some remaining scale-free signatures which reveal the importance of signal aspects. To faithfully represent population activity statistics, a candidate signal must: have the right non-Gaussian distribution, the right 1/*f* power-spectrum characteristics, and be sensitively dependent on higher-order temporal correlations such as may result from the complex interplay of excitation and inhibition within the dendritic arborization of a pyramidal neuron in the visual cortex.

## Discussion

Leveraging *V*_{m} and LFP recordings with modeling and MEA data yielded two principle findings: subthreshold *V*_{m}s are a useful indicator of network activity and this correspondence is inherent to critical coarse-graining. Scrutiny revealed that avalanche size and duration distribution parameters covary to maintain similar geometrical scaling across different experiments, a noteworthy observation. The following discussion emphasizes possible significance and research intersections, such as explaining disagreement with theory via subsampling effects or quasicriticality or relating neural computation to a mathematical apparatus within critical systems theory.

Although “appropriating the brain's own subsampling method” is a novel description of whole-cell recordings, it was inspired by examples. Whole-cell recordings contain information about the network (Gasparini and Magee, 2006; Mokeichev et al., 2007; Poulet and Petersen, 2008; El Boustani et al., 2009; Okun et al., 2015; Cohen-Kashi Malina et al., 2016; Hulse et al., 2017; Lee and Brecht, 2018) and stimulus (Anderson et al., 2000; Sachidhanandam et al., 2013). Usually, the focus is using neural inputs to predict outputs, not to measure population dynamics (Destexhe and Paré, 1999; Carandini and Ferster, 2000; Isaacson and Scanziani, 2011; Okun et al., 2015). Additionally, long-time or large-population statistics, like our avalanche analysis, are useful for understanding neural code (Sachdev et al., 2004; Churchland et al., 2010; Crochet et al., 2011; Graupner and Reyes, 2013; McGinley et al., 2015; Gao et al., 2016) and are robust to noise. Our finding that single *V*_{m} recordings reflect scale-free network activity is significant as recording stability in behaving animals improves (Poulet and Petersen, 2008; Kodandaramaiah et al., 2012; Lee and Brecht, 2018). We open the door to using *V*_{m} fluctuations as windows into network dynamics.

Rigorous analysis supports our experimental conclusion: subthreshold *V*_{m} fluctuations mimic neuronal avalanches and evince critical phenomena, but negative LFP deflections do not despite being purported network indicators (Bédard et al., 2006; Liu and Newsome, 2006; Kelly et al., 2010; Einevoll et al., 2013; Okun et al., 2015). We invoke network not single-neuron criticality (Gal and Marom, 2013; Taillefumier and Magnasco, 2013) because the trend between size and duration exponents agrees with MEA data. Our findings originate from spontaneous activity of ex-vivo turtle visual cortex which shares many connectivity and functional features with mammalian cortex (Ulinski, 1990; Larkum et al., 2008). Last, the results are not serendipitous noise because the *V*_{m} dataset significantly differed from a dataset of phase-shuffled and rescaled surrogates (Theiler et al., 1992).

Readers keen on critical phenomena may notice our exponents differ from the exact theoretical predictions (τ = 1.5, β = 2) (Haldeman and Beggs, 2005). Others observing this mismatch have suggested the brain operates slightly off-critical (Hahn et al., 2010; Priesemann et al., 2014; Tomen et al., 2014).

An extension of this suggestion, quasicriticality (Williams-García et al., 2014), also explains the highly stable scaling relation: biological systems blocked from precise critically may optimize properties that are maximized only for critical systems, becoming “quasicritical.” Correlation time and length are maximized only at criticality and closely related to avalanche geometrical scaling (Tang and Bak, 1988; Sethna et al., 2001). If brains optimize correlation length, then a highly stable scaling relation may result. Furthermore, including inhibition (Larremore et al., 2014) makes our otherwise critical model less consistent with criticality except that population statistics can still be inferred from input fluctuations. The stable scaling was not in the model, which lacks any plasticity mechanisms. Stable scaling may be a rare observation of self-organization principles such as quasicriticality. A contributing explanation is subsampling effects (Priesemann et al., 2009; Levina and Priesemann, 2017), but it does not explain the stable scaling relation unless quasicriticality is also invoked.

### Neuronal avalanche and neural input fluctuation similarity are captured by a critical recurrent coarse-graining network

Our main modeling finding, inputs to a neuron reflect network activity best for critical branching networks, is supported by a parameter sweep and detailed analysis. Our network had no structure, but structure exists at all scales of brain networks (Song et al., 2005; Perin et al., 2011; Shimono and Beggs, 2015) and can have profound impacts on network dynamics (Litwin-Kumar and Doiron, 2012; Mastrogiuseppe and Ostojic, 2018). We derived a relationship showing that the findings may be transferrable to networks where neural inputs fluctuate about proportionality to some subsample's activity. We tune proportionality to be one, but that can also emerge from plasticity (Shew et al., 2015; Del Papa et al., 2017). Tight balance suggests a biological mechanism causing subthreshold *V*_{m} to track excitation into a presynaptic population because IPSPs can have their timing and strength “balanced” to truncate EPSPs that would otherwise last longer than spurts of presynaptic excitation (Barrett et al., 2013; Boerlin et al., 2013; Gatys et al., 2015; Denève and Machens, 2016). We use the *V*_{m} proxy Φ* _{i}*(

*t*), an alpha-function convolved with a point process,

*P*(

_{i}*t*). Φ

*(*

_{i}*t*) is more like

*V*

_{m}than

*P*(

_{i}*t*) and reproduces our experimental findings. Last, we investigate quasicriticality by including inhibition but tuning the maximum eigenvalue to what would be the critical point without inhibition.

Our model provides insights on network subsampling and renormalization group. Usually, subsampling means selecting neurons at random or modeling an MEA with an arbitrary grid (Priesemann et al., 2009). Our “subsample” is the presynaptic population represented by summing weighted inputs from active neurons. This is the first analysis intersecting network convergence (i.e., postsynaptic soma).

Subsampling distorts avalanche size and duration, likely creating differences between experimental results and theoretical predictions (Priesemann et al., 2009; Ribeiro et al., 2014; Levina and Priesemann, 2017; Wilting and Priesemann, 2018). Subsampling may explain disagreement between avalanche analysis on simulated network activity, *F*(*t*), the *V*_{m} proxy Φ* _{i}*(

*t*), and the single-neuron firing rate

*P*(

_{i}*t*). However,

*V*

_{m}and MEA results are off theory but match each other. Either their subsampling errors are alike enough to produce similar distortions or subsampling cooccurs with quasicriticality (Priesemann et al., 2014; Williams-García et al., 2014).

Intriguingly, the restricted Boltzmann machine (RBM) (Aggarwal, 2018) (a related model) was exactly mapped to a “renormalization group” (RG) operator (Mehta and Schwab, 2014; Koch-Janusz and Ringel, 2018). RG is a mathematical apparatus relating bulk properties to minute interactions (Maris and Kadanoff, 1978; Nishimori and Ortiz, 2011; Sfondrini, 2012). It characterizes critical points of phase transitions (Stanley, 1999; Sethna et al., 2001) and helps to derive neuronal avalanche analysis predictions (Sethna et al., 2001; Le Doussal and Wiese, 2009; Papanikolaou et al., 2011; Cowan et al., 2013). RG operators coarse-grain and then rescale, like resizing a digital image. Crucially, iterating an appropriate operator on a critical system produces statistically identical “copies,” but on noncritical systems, the iterations diverge. Our model averages (coarse-grains) presynaptic pools to get an instantaneous firing probability for each neuron; then, a logical operation (rescaling) sets the spiking states for the next iteration, demonstrating an RG-like operation that reproduces our experimental findings. Denève and Machens (2016) proposed a similar relationship between real *V*_{m} and presynaptic pools. The finding that a similar neural operation emerges in RBMs underscores the relevance of RG and the extension of our findings to structured or nonbranching networks. The importance is that a recurrent coarse-graining network may be like a scale-free ouroboros, displaying widespread scale-freeness if any component is critical or briefly driven by critical or scale-free inputs (Mehta and Schwab, 2014; Schwab et al., 2014; Aoki and Kobayashi, 2017; Koch-Janusz and Ringel, 2018).

Significantly, associating neuronal processing with critical branching may induce an organizing principle, the “information bottleneck principle.” This balances dimensionality reduction (compression) against information loss (Tishby and Zaslavsky, 2015) and is reminiscent of efficient coding (Friston, 2010; Denève and Machens, 2016) and origins of tuning curves (Wilson et al., 2016; Heeger, 2017). Koch-Janusz and Ringel (2018) trained their network by maximizing mutual information between many inputs and few outputs. This produced nodes with receptive fields matching popular RG operators. They derived correct power laws by iterating the network. Applications of RG to neural computation are known: image processing (Gidas, 1989; Mehta and Schwab, 2014; Saremi and Sejnowski, 2016), brain and behavior (Freeman and Cao, 2008), emergent consciousness (Werner, 2012; Fingelkurts et al., 2013; Laughlin, 2014), and hierarchical modular networks (Lee et al., 1986; Willcox, 1991) important for criticality (Moretti and Muñoz, 2013). Furthermore, our model's RG-like features are crucial to reproducing experimental results. It follows that elegant RG operators as in the RBM might also capture biological neuronal processing, fulfilling the demand for beautiful neuroscience models (Roberts, 2018) while offering insights into organizing principles and scale-freeness.

### Conclusion

We have established that subthreshold fluctuations of *V*_{m} in single neurons agree with neuronal avalanche statistics and with critical branching, but fluctuations in other single-electrode signals do not. Computational modeling showed that accurate inference requires critical-branching-like connectivity. Fluctuation size scales with duration more self-consistently in experimental than model results, hinting at self-organization. These findings are consistent with a nascent reduction of neural computation to coarse-graining operations that may explain the prevalence of critical-like behavior during spontaneous neural activity. Fully articulating the implications requires more investigation, but we have substantially extended the evidence for critical phenomena in neural systems while rigorously demonstrating that subthreshold *V*_{m} fluctuations of single neurons contain useful information about dynamical network properties.

## Footnotes

This work was supported by the Whitehall Foundation (Grant 20121221 to R.W.) and the National Science Foundation (Collaborative Research in Computational Neuroscience Grant 1308159 to R.W.). We thank Woodrow Shew at the University of Arkansas for helping with the comparison with MEA data.

The authors declare no competing financial interests.

- Correspondence should be addressed to James K. Johnson at jkjohnson{at}wustl.edu