Deconstructing scale-free neuronal avalanches: behavioral transitions and neuronal response

Observations of neurons in a resting brain and neurons in cultures often display spontaneous scale-free (SF) collective dynamics in the form of information cascades, also called ‘neuronal avalanches’. This has motivated the so called critical brain hypothesis which posits that the brain is self-tuned to a critical point or regime, separating exponentially-growing dynamics from quiescent states, to achieve optimality. Yet, how such optimality of information transmission is related to behavior and whether it persists under behavioral transitions has remained a fundamental knowledge gap. Here, we aim to tackle this challenge by studying behavioral transitions in mice using two-photon calcium imaging of the retrosplenial cortex (RSC)—an area of the brain well positioned to integrate sensory, mnemonic, and cognitive information by virtue of its strong connectivity with the hippocampus, medial prefrontal cortex, and primary sensory cortices. Our work shows that the response of the underlying neural population to behavioral transitions can vary significantly between different sub-populations such that one needs to take the structural and functional network properties of these sub-populations into account to understand the properties at the total population level. Specifically, we show that the RSC contains at least one sub-population capable of switching between two different SF regimes, indicating an intricate relationship between behavior and the optimality of neuronal response at the subgroup level. This asks for a potential reinterpretation of the emergence of self-organized criticality in neuronal systems.


Introduction
Ever since the initial findings of Beggs and Plenz [1], the scale-free (SF) statistics of the temporal organization of neuronal activity (dubbed neuronal avalanches) has been used as evidence toward the critical brain hypothesis. It is conjectured that, in an effort to maximize information processing capabilities, the brain is self-tuned to a point, or regime, perched between an ordered phase and a disordered phase [2][3][4]. From statistical physics it is known that a necessary condition of systems near a phase transition are SF statistics for various system properties, e.g., the size and duration of an avalanche [5]. Since those initial findings, substantial experimental and theoretical efforts have been dedicated toward studying the presence of SF statistics in neuronal dynamics (for a recent review, see [6]).
SF statistics have no characteristic size, and as such, one expects to observe the same statistical behavior (i.e., the same critical exponents), regardless of the scale at which the system is observed. Studies of the brain from the neuronal level, up to the whole-brain level in fMRI have all displayed signatures of SF statistics during ongoing, spontaneous activity [7][8][9][10]. This phenomena has also been studied in multiple species [8,11,12], as well as in vitro [1,13], using different experimental setups [14,15]. In many of these studies, regardless of scale and species, researchers have found similar critical exponents for the size distribution of spontaneous neuronal avalanches, often with an avalanche size exponent near −3/2, though not always. Deviations away from criticality have been tied to changes in the underlying state of the subject, such as changes in consciousness [15,16], or disruptions to homeostatic firing rates [17]. The value of −3/2 is consistent with, and has been used to motivate, the mean-field branching process (MFBP) as a candidate universality class for neural dynamics [18]. This suggests that normal neuronal dynamics belong to the same universality class as a stochastic branching process, while deviations from 'normal' neuronal dynamics may yield different statistics, potentially indicating different universality classes [13], or a loss of criticality altogether [7,15]. More recent efforts have been devoted to understanding why different experiments demonstrate different exponents, such as how random sub-sampling may result in different avalanche statistics [19], even if the underlying dynamics are fixed [20]. Other recent studies have focused on the role of time-scale separation and spontaneous activity may contribute to the observed SF statistics [21], and by extension the interpreted universality class [22].
Neuronal activity can be broadly classified into two sub-types: ongoing spontaneous activity, and evoked transient activity. The statistics of neuronal avalanches of spontaneous activity has largely consistently shown to follow SF statistics, at least in the awake state [7,15]. In anesthetized subjects there is an observed loss and subsequent recovery of SF statistics, suggesting that SF statistics are a hallmark of the dynamics necessary for day-to-day waking activity [7,15]. A more controversial topic is the presence of SF statistics during transient periods, such as during task based experiments. While some studies have suggested that avalanche distributions during transient non-spontaneous activity do not display SF dynamics [8,[23][24][25], others have found they do [26]. One specific case study using LFP data in the monkey cortex has shown that using an adaptive time binning and thresholding approaches recover SF dynamics that would otherwise not be present [27]. The authors argue that the adaptive methods are necessary to account for the non-stationary levels of neuronal activity associated with behavior.
Except for in the case of whole brain studies using fMRI, neuronal avalanche statistics are traditionally generated from a localized sub-sample of neurons. Recent research has shown how the avalanche statistics of a sub-sampled population may change with respect to the total population [19,20]. However, what is typically not considered in task-based avalanche studies is the diverse neuronal response, even within a localized subpopulation of neurons. As such, how the diverse response of neurons to behavioral transitions contributes to the generation of neuronal avalanche statistics and reflects the optimality of the information processing is an open question.
In this study we demonstrate that neurons of the mouse retrosplenial cortex (RSC) can be categorized with respect to their response to behavioral transitions from resting to moving states: we find two distinct subpopulations; cells whose firing rate is suppressed during periods of moving, and cells which increase activity during periods of moving. In order to understand how avalanche statistics are generated or affected by these sub-populations we take a top-down approach-we initially consider the total population, and then proceed to deconstruct these avalanches. First temporally by conditioning on behavior, and then further deconstructing the population into dynamical sub-groups, along with their corresponding conditional and unconditional avalanche statistics. We find that the unconditional statistics of avalanches generated from neurons belonging to different sub-populations largely coincides with the behaviorally conditioned statistics of the total population, indicating two fundamentally different classes of approximately SF responses. In particular, the subgroup of neurons that increase activity in the moving state can switch between these two dynamical population responses associated with behavior. Without access to the underlying network structure to establish causality between firing events, these diverse dynamics are typically mixed if the total population is considered to generate effective avalanche statistics. Our findings indicate that in general one needs to take into account such superposition effects in order to correctly interpret the presence (or absence) of SF behavior and the value of the observed, potentially effective exponents during cued or self-initiated behavior.

Data
Adult (9-11 months old) Thy1-GCaMP6m mice (n = 7), were housed in standard rodent cages, and maintained at 24 • C under a 12 h light/dark cycle, with experiments performed during the light cycle. Procedures were in accordance with the Canadian Council on Animal Care guidelines, and with protocols approved by the Animal Welfare Committee of the University of Lethbridge (protocol #1512).
Mice received a 5 mm bilateral craniotomy (AP: +1 to −4; ML: −2.5 to +2.5), which was then covered with three layers of coverslips attached to the skull using Vetbond. A titanium head plate was fixed to the skull using dental acrylic. Head-fixed mice were trained to run on a treadmill and received a drop of 10% sucrose solution on every lap of the treadmill belt. Animals were water restricted during the training and testing (with 30 min per day free access to water) and their body weight was carefully monitored throughout the experiment to ensure the weight loss did not exceed 15% of their baseline value. Training continued in daily sessions until mice performed at least 20 trials in 30 min. The treadmill belt was 150 cm long and 4 cm wide with four tactile cues in different locations. The reward was delivered in the same location on the belt for all trials by a solenoid valve. An optical encoder attached to the wheel shaft was used to monitor belt movement. The behavior was acquired with Clampex software (Axon Instruments). Neural activity was imaged in daily sessions of 20-30 min while the mice were performing the task. Automatic image pre-processing was performed using the Suite-2P algorithm [28] for MATLAB (MathWorks), as described in [29]. The regions of interest (ROIs) detected were inspected manually and labeled as cells or non-cells by experienced users.
We measured continuous neuronal activity via two-photon (2P) imaging of endogenously expressed calcium indicators (GCaMP6) in the RSC. The field of view of the recording was 835 μm × 835 μm, with a depth between 135-170 μm (layers II/III). The 2P scanning rate was 19 Hz. The average number of neurons over all recordings was 418 ± 162 neurons. We studied both unconditional and behavior-conditioned avalanche distributions. Each mouse has four 20 min long recordings obtained on consecutive days.
During the course of the experiment the mouse is free to move (or stand) on a non-motorized belt, while head-fixed to the 2P imaging system. A liquid reward is administered upon completion of one lap of the belt (150 cm). The track speed is recorded and later thresholded to indicate behavioral transitions (see section 2.2.1). We consider only the resting-to-moving (R-to-M) behavioral transition in our analysis because (a) the moving-to-resting (M-to-R) transition is confounded with reward consumption and thus may necessitate more than just the binary 'resting' and 'moving' states to adequately describe, and (b) the M-to-R transition is also typically not as sharp as the R-to-M transition and thus also may not be adequately captured by a two state approximation.

Neuronal avalanches and scale-free properties
The self-organized patterns of neuronal activity are typically referred to as neuronal avalanches. The size S of the avalanche can be defined as the number of above threshold events (e.g., neuronal firings) that participated in that avalanche, and the duration, T, is the time over which the avalanche was active. Often the distributions (i.e., probability density functions) of these statistics have been found to be power-law, (also called SF): The deconvolved calcium signal obtained from the Suite-2P algorithm is not binarized. To generate a more straightforward interpretation of the total signal in terms of the number of active neurons, we apply a binarizing threshold to each neuron individually (see figures 1(a) and (b)). No association between the signal intensity of the deconvolved calcium signal and the cells correlation with behavior was observed. From the binarized signal, avalanches are obtained by summing across all thresholded ROIs as consecutive frames of non-zero activity, delineated by periods of inactivity. Distributions were obtained by concatenating statistics across all recording sessions of that animal. We typically do not concatenate distributions across mice however as we found that there was too much variation.
The number of quiescent frames, and by extension the number of avalanches, depends on the threshold. If the threshold is too large, few neuronal events are registered. If the threshold is too small, periods of quiescence become infrequent, merging together what might have otherwise been considered multiple neuronal avalanches. As such we choose the threshold, θ max , so that number of avalanches (or equivalently, the avalanche rate) is maximized [7] for each recording, as shown in figure 1(c). Defining the mean deconvolved calcium activity as η avg = N −1 N i=1 η i , where η i is the deconvolved calcium signal of neuron i averaged over the recording, the inset of figure 1(c) shows that larger population activity results in larger thresholds, as expected. Figure S1 (https://stacks.iop.org/JPCOMPLEX/2/045010/mmedia) of the supplementary material shows the range of thresholds over which robust power-law exponents were found.
The largest size of the avalanche depends on the number of neurons, N, and the average population (binarized) activity λ avg = N −1 N i=1 λ i , where λ i is the mean number of binary firings of neuron i, averaged over the duration of the recording. To compare across recordings and animals, we normalize the avalanche sizes by Λ = Nλ avg [7]. For simplicity we will call S/Λ the avalanche size. All avalanche statistics of the total population presented in the following are largely unchanged if N is reduced by performing random sub-sampling (see, e.g., figures S2(a) and (b)), partially due to the (scaling) behavior of Λ (see, e.g., figures S2(d) and (e)).
SF statistics in the size and duration of avalanches is a necessary condition for systems at a critical point [5]. Note however, unlike critical systems such as the Ising model, or the branching process which need to be tuned to the critical point, neuronal systems are out-of-equilibrium systems that are thought to self-organize toward a critical point [30][31][32]. For critical systems, the mean duration of an avalanche depends on the size S; with the exponents obeying the following scaling relationship: As both sides need to be satisfied independently, the equality of the scaling relation has been used to test for critically [33]. For simplicity, we refer to the right-hand side of equation (4) as rhs throughout this document. The set of exponents obtained from the power-law distributions of various statistics specify the universality class the system is said to belong to. An oft-cited potential universality class for spontaneous on-going neural activity is the MFBP class, which corresponds to the set of exponents τ = 1.5, α = 2, and γ −1 = 2 [1].

Behaviorally conditioned avalanches
Generating conditional avalanche distributions is achieved by taking five avalanches immediately prior to and following each R-to-M transition, for the pre and post-transition distributions respectively. Avalanches occurring during the transition are discarded to avoid ambiguity. We chose five avalanches so as to maximize the statistics, while at the same time remaining below the mean moving/resting duration so as to not sample too close to the either the previous or next moving-to-resting transition. Avalanches extending past another behavioral transition are discarded.
To identify behavioral transitions we use accompanying track speed data. We first process the track data by applying a sliding mean filter with a width of 40 frames (roughly 2 s). This is to avoid very brief transitions.
A mean filter elongates periods of activity, but does so consistently by an amount roughly determined by the width of the window. We compensate for this simply by applying a constant shift to each data set. In this case we found that ten frames was sufficient. We thresholded for each track recording (using the tenth percentile) to identify transitions. In some rare instances brief occurrences of spurious backward movement can occur. We consider any such track data to be below threshold. See figure S3 and the surrounding discussion an analysis of the quality of the automatic detection.

Adaptively thresholded avalanches
Adaptively thresholded avalanches were generated by taking a 41 frame (≈2 s) window centered on R-to-M transitions (i.e., first 20 frames are the resting state, last 20 frames are the running state). Within each 20 frames, a local threshold is calculated that maximizes the number of avalanches. This also includes calculating a local Λ for before and after each transition. We repeat this over all transitions within a recording.

Generating behaviorally-regulated subgroups
Cells are classified as either M (for 'moving'), or R (for 'resting'), according to the correlation between the corresponding deconvolved calcium signal and the track speed data. We define the 15 percent of cells with the most positive correlation belonging M group and the 15 percent with the most negative correlation as belonging R group. An example of the time series of an M and R cell, and how they correlate with the track, are shown in figure 1(d). We will denote the total population as the A (for 'all') group. We chose this approach to ensure that each group had the same number of cells. A significance based approach is shown in figure S4, but the result did not change. All cells appear to display similar activity in the standing phase.

Estimation of critical exponents
Power-law analysis is typically done over a range between a minimum and maximum avalanche size (S min and S max , respectively, see also section 2.4.2). Ranges over which the power-law analysis was done are displayed in figures as dashed lines. We calculate and cross-validate values of avalanche exponents using several independent methods. Avalanche exponents were calculated using maximum likelihood estimation (MLE) methods (mle, MATLAB) [34,35]. We also independently performed a rescaling analysis to extract the exponent. For a generic random variable X that follows a power law P(X) = CX −k (C is the normalization constant), rescaling the lefthand side by X k yields X k P(X) = C, meaning that the correct estimate for k is the one that minimizes the slope of the plot. This method is useful for cross-validating α as we found for distributions with larger slope and limited range MLE can over-estimate the exponent (figure S6(h)). However, in most cases (e.g., when the exponent is shallower and by extension the fitting range is sufficiently large) the MLE method and rescaling method agree very well (figure S6).
The scaling relation, equation (4), defines the relationship between the exponents τ and α that is expected to be obeyed at critically. We calculate 1/γ by first binning the avalanche sizes and finding the mean of the corresponding avalanche durations to get T (S) (see figure 2(g)). We then use least squares to fit a line between s min and s max (in log-log space) to T (S). The slope of this line corresponds to γ.
We calculate the right-hand side of equation (4) (rhs) using the exponents obtained from the rescaling method. This is because we found MLE often failed to give a reasonable estimate when the range is reduced, as is often the case for duration distributions. To cross-validate this result, we compare the α obtained from the rescaling method against the value one would expect to get from τ (obtained via MLE) and γ, and assuming the scaling relation (equation (4)) held. This second method is independent of rescaling and as the right column of figure S6 shows, the two methods are consistent with one another.

Uncertainties
We assume two sources of uncertainty for the MLE method. The first is the inherent uncertainty of the MLE, which was taken to be half of the 95 percent confidence interval. The second is the uncertainty in the fitting range. To estimate this, we vary the fitting interval in the lower range (where S min ∝ O(1)) by ±1, and in the upper range (where S max ∝ O(10 2 )) by ±10. The range of exponents obtained doing this was taken as the uncertainty. We then add these two uncertainties. Typically we found that the intrinsic uncertainty was greater than the fitting uncertainty.
We also assume two sources of uncertainty for the rescaling method. The first is the intrinsic uncertainty of estimating when the graph was 'flat'. We took ±0.1 as an estimate for this uncertainty. The second is that this method depends on the binning used as it is derived from P(S/Λ). We test for different bin sizes by varying the number of bins from N bins (1 − 0.2) to N bins (1 + 0.2), and take the half the interval from the minimum to the maximum value of α obtained.
The uncertainty on γ, denoted δγ, is the 95 percent confidence intervals of the fitting method (fit, MAT-LAB). The uncertainty for γ −1 is obtained using the standard error propagation formula, δ(γ −1 ) = γ −2 δγ. Similarly, uncertainties on rhs as well as α BC were also obtained via the standard error propagation formula. For example, the uncertainty on rhs is given by: where δα and δτ are the uncertainties on α and τ .

Evaluating goodness-of-fit
To test whether power-law behavior is a plausible fit we first establish a range of avalanche sizes, [S min , S max ] (using S as a short-hand for S/Λ) is identified so that p > 0.1 (pvcalc NCC toolbox [36]). S max was typically set by the (visually identified) cut-off of the distribution, and then S min was chosen so as to give the largest possible range. Next we tested against alternate hypothesis similar to a method similar to that shown in [20]; by first maximizing the likelihood function, denotedL, of the alternative distribution, then calculating the Akaike information criterion, AIC X = 2k − 2 lnL, where k is the number of parameters for the hypothesis X [37]. Common alternatives tested were the exponential, stretched-exponential, and log-normal distributions (see figure S7 and table S1). A value of AIC close to zero suggests a more parsimonious model. We use the metric ΔA = AIC X − AIC pl , where AIC pl corresponds to the power-law, as a measure of plausibility: ΔA > 0 suggests the power-law as a likely model, and ΔA < 0 suggests X as the likely model. Finally, we define the AIC score as the number of mice (out of seven) for which power-law was the favored distribution (see table S2). This score approach was taken due to significant variation in the actual value of ΔA from animal to animal, and the limited statistics which skewed means and uncertainties toward extreme outcomes.

Unconditional avalanche statistics
First we studied behaviorally unconditional avalanche statistics of the total population. This analysis serves as a baseline and reflects the results one would obtain without behavioral information. We find that P(S/Λ) and P(T) follow power-law like distributions over extended ranges with a slope of approximately τ = 2.0 ± 0.2, and α = 2.5 ± 0.2 (figures 2(a) and (b)), where uncertainties are given by the standard deviation across the n = 7 animals. Alternate hypothesis testing favored power-law over exponential and stretched-exponential, and slightly favored power-law over log-normal (see table S2). We see that the exponents derived from behaviorally unconditioned avalanches do not match MFBP expectations. We compared avalanche statistics against surrogate data generated by randomly cyclically permuting the deconvolved calcium signal of each cell, which are shown in gray. The surrogate avalanches are statistically different according to a two-sided Kolmogorov-Smirnov (KS, denoted D KS ) test (kstest2 in MATLAB, D KS ≈ 0.33, p 10 −10 for P(S/Λ), and D KS ≈ 0.07, p 10 −9 for P(T), where p indicates largest p-value obtained). Note that for the surrogate data we tested both calculating a new binarizing threshold for each surrogate (shown in figures) as well as using the same threshold as the real data. The resulting distributions were similar in each case (not shown).

Avalanche statistics conditioned on behavioral transitions
Next we focus on avalanche distributions conditioned on behavioral transitions (see section 2.2.1 for details). We denote behaviorally conditioned distributions by P(S|resting) and P(S|moving), but emphasize here that 'resting' and 'moving' in this context refer to the periods of neuronal activity strictly in a window around an Rto-M transition, meaning that not all times in the recording are designated as moving or resting. As mentioned in section 2.1, this was done so as to avoid including confounding factors such as reward administration upon completion of a lap. We did however test whether these windows were representative of the results one would find had we designated all times as either being running or standing, and indeed found this to be the case ( figure  S3). In other words, the following holds: P(S) ≈ P(S and resting) + P(S and moving), where the conditional distributions can be obtained from P(S and X) = P(S|X)P(X), where X is a stand-in for the 'moving' or 'resting' states. Figures 2(c) and (e) shows the conditional resting and moving size distributions, respectively, and (d) and (f) show the respective avalanche duration distributions. In the resting state we found that τ = 2.5 ± 0.5 and α = 3.1 ± 0.3, where as in the moving state we found τ = 1.5 ± 0.3 and α = 2.0 ± 0.2, which are close to the MFBP exponents. Alternate hypothesis testing favored power-law over exponential, stretched-exponential and log-normal in both the moving and resting states (see table S2).
As with the unconditional distributions we compared the distributions against surrogate data. While the differences between the real distribution and the surrogates are clear in figures 2(e) and (f), they are not as apparent in (c) and (d). Nevertheless, except for the duration distribution for a single animal, all were found to be statistically different according to a two-sided KS test (p = 0.002 was largest significant test among duration and size distribution comparisons). Figure 2(g) shows one example of the relationship between avalanche size and avalanche duration, in the moving state. We found that the scaling relations held for the unconditional and resting distributions (figure 2(h)). We found that in the moving case, however, that the difference between the two is much larger, but remains within uncertainty. For the unconditional case we find that γ −1 = 1.34 ± 0.08 and rhs = 1.51 ± 0.4. In the resting case γ −1 = 1.4 ± 0.1 and rhs = 1.36 ± 0.6. In the moving case γ −1 = 1.32 ± 0.08 and rhs = 2 ± 1. We note here that the larger uncertainty on running state value is due to τ being closer to unity, resulting in larger uncertainties due to the functional form of the rhs uncertainty formula shown in equation (5). The values of γ −1 differ from the MFBP expectation of γ −1 = 2, but do appear to be consistent with other reported values, as in [7,38].

Neuronal response to behavioral transition is diverse
Here we show that firing activity of the neurons is modulated by the resting/moving behavior, but not uniformly so across the entire population. Figure 3(a) shows a representative example of the thresholded population activity in a window averaged over all R-to-M behavioral transitions for a single recording (more examples in figure S4(a)). The cells are ordered in terms of the Pearson correlation between neuron activity and the track speed data. We see a subset of cells increases its activity when transitioning into a moving state, while another subset of cells suppress' its activity in the moving state. Cell groups are robust with recordings taken roughly an hour apart from one another, suggesting that these groups are not simply due to ordering random activity ( figure 3(b)). Comparisons between days could not be made as the cells are not always the same from recording to recording.

Unconditional avalanche statistics of sub-groups
As the activity of the M and R sub-groups might depend on the state of the animal we need to control for the discrepancy in the amount time spent in either the moving or resting state. To do so, we again consider only periods of neuronal activity strictly in a window around an R-to-M transition and consider the unconditional avalanche distribution to be the 'moving' and 'resting' conditional distributions joined together (i.e., P(S) = P(S and resting) + P(S and moving)). Figure 3(e) shows that the size of the unconditional avalanches for the A group followed a power-law like distribution with slope τ A = 1.8 ± 0.1. This is consistent with the result in  In figure 3(c) we plot the avalanche duration exponent α against the size exponent τ for the individual mice (n = 7). The aforementioned means of τ and α of each subgroup are represented by large markers. The star marker indicates values one would expect from MFBP. This suggests that sub-sets of cells with different responses to changes in behavior display power-law statistics that are different from each other and from the total population. Despite this, the M group is clearly represented much more by the A group statistics, suggesting this sub-group dominates the observed statistics. Finally figure 3(d) shows good agreement for the scaling relations.
Comparisons against cyclic permutation surrogates were qualitatively similar to that of the corresponding conditional distributions in figure 2-all were statistically different from the generated surrogates (p < 0.05, kstest2 MATLAB). For visual clarity we omitted these surrogates from figure 3. We also considered surrogate avalanches generated by randomly sub-sampling cells. We found that randomly sub-sampled distributions do not generate the statistics observed in the R and M groups (figures S2(b) and (c)).

Conditional avalanche statistics of sub-groups
One possible conceptual mechanism that could generate such behavior dependent statistics, is that the activity of the total population is switching between these two sub-populations. To test this we analyze the behaviorally conditioned avalanches of the sub-groups, which were generated in the same manner as the total population analysis. Figures 4(a)-(c) shows the conditional distributions in the moving (blue) and resting (red) states. While the resting state distribution was similar for all three groups, the moving state avalanche statistics of the R group diverge from the M and A groups substantially. Figures 4(e) and (f) show the avalanche duration exponents plotted against the size exponents in the resting and moving states, respectively. Large shapes indicate the mean positions and the large star represents MFBP expectations. In the resting state the size exponent τ was similar for all three groups, with a value of τ ≈ 2.4. In the moving state we found that the size exponents are τ A = 1.5 ± 0.3, τ R = 2.9 ± 0.7, and τ M = 1.4 ± 0.3. Similarly, the duration exponent in the resting case is similar for all groups with α ≈ 3. In the moving case we found the exponents α A = 2.0 ± 0.2, α R = 3.2 ± 0.4, and α M = 1.6 ± 0.3. In the moving state, alternate hypothesis testing favored power-law over other test distributions for both the M and R groups. In the resting phase the scores drop-while power-law is slightly favored for the M group the test suggests an exponential distribution as a more likely distribution in the R group (table S2).
Together this suggests that instead of a switching phenomena, one subset of cells (namely the M cells) dominates the avalanche statistics during the animals moving phases, whereas during the resting phase both sub-sets have the same dynamics (figures 4(d) and (e)). In terms of the dynamics observed at the total population level, our results suggest that exponents obtained from the total population in the resting state are representative of both sub-sets of neurons. Yet, this is not true in the moving state, where the exponents are largely representative of only one sub-set of dominating cells.
Finally, both sides of the scaling relations in equation (4) are plotted in figure 4(f). Just as in our previous cases, we found the moving case to deviate the most, but was still within uncertainty. Surrogate analysis was also performed. Using random cyclic permutations we found that, except for in the case of two trials of the R cells in the resting state (kstest2 MATLAB, p = 0.34 and p = 0.33 for those two trials), all distributions were found to be statistically different from their corresponding surrogates (largest p-value p = 10 −3 ).

Avalanche statistics with an adaptive threshold
A recent study [27] suggests that the non-stationary levels of neuronal activity found in task-based studies require adaptive methods in order to capture the intrinsic SF dynamics. We investigate if adaptive thresholding can account for the differences between pre and post-transition avalanches in the different cell groups. Figure 5(a) shows that applying the adaptive threshold generates a more homogeneous response across the transition, as further quantified in figures S5(a) and (b). Figures 5(b)-(d) shows the avalanche size distributions. In order to reduce statistical fluctuations between the avalanches obtained from the pre and post-transition states, we concatenated the avalanches of all animals. Uncertainties were obtained using 400 bootstrap samples. Red distributions correspond to the resting state, whereas blue corresponds to the moving. Because Λ is calculated locally, the smallest avalanche size, O(Λ −1 ), becomes a highly fluctuating variable. As such the figure domain shows only the sizes for which the these fluctuations become small. In the moving case only the A group and M group displayed potentially power-law statistics with τ A = 1.8 ± 0.2 and τ M = 1.9 ± 0.1, (p = 0.12 and p = 0.10, respectively, ranges chosen so as to maintain p > 0.1). In the resting state we found only the A group displayed SF statistics (τ = 1.7 ± 0.1 with p = 0.15), but the range of validity is severely reduced when compared to the other two. In all other cases we did not find SF statistics (p < 10 −4 ). Within the statistical uncertainties, the above exponents are similar to one another and bare resemblance to the exponents shown in figure 3 such that it is possible that they are indeed identical. We note that our findings are specific to the behavioral periods: neglecting any information about behavior (e.g., the track speed) and applying adaptive thresholding using a rolling adaptive threshold, the avalanche distributions take on a fundamentally different form (see figure S5(d)). This implies that access to or observation of the behavioral state is essential to understand the emerging avalanche dynamics.

Discussion
Using 2P calcium imaging of the RSC in the mouse brain, we analyzed behaviorally conditioned and unconditioned neuronal avalanches. Within this population of neurons, we identified robust sub-sets of cells with behaviorally dependent dynamics. In the resting phase of behavior all cells had similar dynamics. Indeed, choosing between them is akin to random sub-sampling ( figure S2(c)). However, the cells differentiate themselves in the moving phase where-in the 'resting' (R) group had its activity down-regulated and the 'moving' (M) group has its activity up-regulated. Because of this selective up/down-regulation, avalanches in the moving phase observed on the total population level are in fact almost entirely due to the smaller M sub-group. This preferential firing (or lack-there-of) of neurons in the moving phase results in avalanches that differ from what one would expect from random sampling [19,20]. As a consequence of these dynamics, we found the behaviorally unconditional avalanches derived from these sub-networks closely resemble the conditional avalanches of the total population.
Effective information processing in the brain requires a diverse portfolio of neuronal activity. This is particularly true of association cortices such as the RSC, which contains cells coding for a variety of features including velocity and head direction etc [39], sensory events [40], as well as cells which mimic the place cells of the hippocampus [29]. One potential explanation would be that the RSC acts as a hub between various parts of the brain, and therefore is reflecting the transfer of information from one region to the other [39]. Indeed, the RSC is well positioned to integrate sensory, mnemonic, and cognitive information by virtue of its strong connectivity with the hippocampus [41], medial prefrontal cortex [42], and primary sensory cortices [40]. However, even highly specialized areas such as the visual cortex, which is not typically interpreted as a hub, not only display a diverse response to visual stimuli, (e.g., tuning curves [43]), but can also directly or indirectly reflect motor activity [44], with the reverse (that other areas of the brain respond to visual stimuli) being true as well [45].
While the presence of movement-correlated neurons in the RSC is not in and of itself surprising, it is interesting that we also find robust groups of anti-correlated neurons. Indeed, to the best of our knowledge the existence of behaviorally anti-correlated cells within the RSC has not been established or emphasized in the literature-in contrast to other areas of the brain [46,47]). In particular, the presence of both movementcorrelated and anti-correlated neurons lends support to the hypothesis that the RSC is structurally a hub, with neuronal activity excited and suppressed appropriately so as to transfer information from one region to another [39,48], especially between the hippocampus and neo-cortex [49]. Why a hub would benefit from critical dynamics is not presently clear. However as computational models suggest, networks at criticality are able to learn complex tasks more easily as the possible response to a given input becomes more diverse [50].
Potentially a hub may benefit from these same principles-routing incoming information more efficiently due to the abundance of possible firing patterns provided by the critical state. Further work would be required to establish whether these R and M cells are indeed a manifestation of the functional networks within which they are embedded, or instead transient neural assemblies [51].
We re-emphasize that it is only in the moving state that the differences in the sub-network dynamics are observed (see figure 3(a) as an example). From a dynamical systems theory perspective, this would suggest that movement induces a bifurcation from ergodic network dynamics wherein all cells are accessible (in the resting state), to non-ergodic dynamics wherein a sub-set of cells (in this case R cells) become inaccessible [52]. Equivalently this could suggest M cells constitute a stable attracting sub-network in the running state [53]. We also found that behaviorally correlated M cells outnumber the anti-correlated R cells ( figure S4(b)). Why one group would contain more cells than the other is conceptually consistent with the notion of sparse orthogonal representation in RSC [29] (i.e., minimizing both the number of active units and correlations between said units), which is thought to be important for efficient coding in associative networks and facilitating navigation [29]. Simply put, an even number of both behaviorally correlated and anti-correlated cells would neither be sparse, nor orthogonal, and would therefore decrease the number of unique activation 'patterns' available to the network.
Despite the fact that it is known that the response of neurons to stimuli or behavior is diverse (even within a small region of interest), it is not well understood how the variety of neuronal responses related to behavior manifests itself in the overall information processing and specifically in the properties of neuronal avalanches. The movement induced suppression and enhancement of R and M cell activity (respectively), suggests that some cells participate in avalanches selectively, which is in line with the results reported in [54]. Furthermore, as previously mentioned, the unconditional avalanches derived from the sub-populations closely resemble the conditional avalanches of the total population. Therefore, just as some cells participate in avalanches selectively, avalanches within a sub-network of the total network participate in behavior selectively. This would support the notion that neuronal avalanches reflect transient cell assembly formation [51]. As it is known that neuronal firing patterns in response to sensory stimuli are preserved in spontaneous activity [55], future studies could also probe as to whether the M and R groups are similarly reflected in spontaneous network activity (e.g., during anesthesia).
Recent research regarding maintained criticality in the context of behavior has suggested that adaptive binning (or thresholding) is required to observed SF statistics in both the pre and post-task avalanche distributions (with the same exponent) [27]. It was argued that non-stationary levels of neural activity require time-dependent binning/thresholding. Without this adaptive aspect, the authors show the pre and post-task neuronal avalanches assume SF statistics and super-critical statistics, respectively. This differs from our analysis of the total population with a non-adaptive threshold in that we do not observe clearly super-critical distributions in the moving state (our analogue of the task state). One possible reason for this is that the resting and moving states are elongated in time, as opposed to the task in [27], which is relatively more brief and transient. As shown in [23], network response to the input of strong sustained stimuli (in our case self-initiated movement), is initially super-critical, but transitions rapidly toward criticality as time progresses.
Interestingly, while adaptive thresholding reduces the differences between the pre and post-transition avalanche distributions when compared against static thresholding (as in [27]), not all sub-populations generated SF statistics under adaptive thresholding. This suggests that adaptive methods are not in and of themselves sufficient to recover critical distributions. However it is currently not clear what delineates whether or not SF statistics can be observed through adaptive thresholding. In [38] it was found that SF statistics are lost when the temporal diversity, as measured by the coefficient of variation, is too low. One possibility might be that an analogous process occurs when cells are non-randomly sub-sampled, resulting in a relatively highly correlated population of neurons and subsequently a loss of neuronal response diversity. However, due to the significantly poorer temporal resolution and limited recording duration (20 min at 19 Hz versus 200 min at 30 kHz in [38]), our current data set does not allow us to meaningfully test this hypothesis.
Finally, we discuss the results presented here in the broader context of the critical brain hypothesis and behavior. As the M group shows in the running state, avalanche analysis of the total population can predominately reflect a single sub-population whose dynamics dominate the time series. However, it is not necessarily clear that would be true in general. In experiments involving large populations of neurons, thresholding is typically required as periods of quiescence are rare, unless the network topology is known such that causality between firings can be established [56,57]. It is possible that this is a manifestation of sub-populations of neurons which individually may display SF statistics, but are anti-correlated to other sub-populations. Such a pair-wise correlation analysis has recently been proposed as a potential coarse-graining procedure for neural populations [58], and future progress in that area may elucidate this idea more. Alternatively, it is possible that sub-populations are not necessarily anti-correlated with one another, but rather are correlated to behavior that is not being recorded or tracked. While we distinguish between two states (resting and moving), it should not be taken to mean that there is not other on-going behavior. Periods of 'resting' should not be interpreted as 'inactive' as these periods may contain whisking behavior, grooming, changes in pupil diameter, tail movement etc, which may be reflected in sub-populations that are different to (but not necessarily independent from) the R and M groups identified here. These behaviors may occur sequentially or simultaneously, but without additional information about behavior [44], or network structure [22], the neural activity that codes for these different behavior is mixed together. Given that spontaneous activity is not entirely random [59], it may be prudent to ask whether 'spontaneous' activity is truly spontaneous, or if it only appears to be so because the relevant underlying variables are not being measured. Indeed, even a sleeping brain displays state transitions (e.g., from non-REM to REM [60]). Viewing the problem from this perspective may provide insight to reconciling criticality during spontaneous and evoked (i.e., behavioral) brain activity [27].
In summary, the critical brain hypothesis has provided an attractive framework with which to understand large scale neural activity. However, tying criticality to behavior has thus far proven to be challenging. We show here that even in very local populations of neurons, there are cells that are both correlated and anticorrelated to behavior. Avalanches within these sub-groups may not resemble avalanches observed at the total population level, suggesting that the observed exponents are in fact effective exponents generated by the mixing of diverse neuronal dynamics. Based on these findings, we believe that in order to understand the differences between spontaneous and evoked neural avalanches in large populations of neurons further one needs to take the structural and functional network properties of these sub-groups into account.