Scalefree behavioral dynamics directly linked with scalefree cortical dynamics
Abstract
Naturally occurring body movements and collective neural activity both exhibit complex dynamics, often with scalefree, fractal spatiotemporal structure. Scalefree dynamics of both brain and behavior are important because each is associated with functional benefits to the organism. Despite their similarities, scalefree brain activity and scalefree behavior have been studied separately, without a unified explanation. Here, we show that scalefree dynamics of mouse behavior and neurons in the visual cortex are strongly related. Surprisingly, the scalefree neural activity is limited to specific subsets of neurons, and these scalefree subsets exhibit stochastic winnertakeall competition with other neural subsets. This observation is inconsistent with prevailing theories of scalefree dynamics in neural systems, which stem from the criticality hypothesis. We develop a computational model which incorporates known celltypespecific circuit structure, explaining our findings with a new type of critical dynamics. Our results establish neural underpinnings of scalefree behavior and clear behavioral relevance of scalefree neural activity.
Editor's evaluation
This paper is an important study that is of interest to neuroscientists studying the organization of neural activity and behavior. The authors present compelling evidence to link the apparently scalefree distributions of behavioral metrics with scalefree distributions of neural activity. They then explore computationally mechanistic models that could account for these observations. The simulations of mechanistic models are provocative and suggest interesting networkconnectivity hypotheses to test in future experiments.
https://doi.org/10.7554/eLife.79950.sa0eLife digest
As we go about our days, how often do we fidget, compared to how frequently we make larger movements, like walking down the hall? And how rare is a trek across town compared to that same walk down the hall? Animals tend to follow a mathematical law that relates the size of our movements to how often we do them.
This law posits that smalltomedium movements and largetohuge movements are related in the same way, that is, the law is ‘scalefree’, it holds the same for different scales of movement. Surprisingly, measurements of brain activity also follow this scalefree law: the level of activation of a group of neurons relates to how often they are activated in the same way for different levels of activation.
Although body movements and brain activity behave in a mathematically similar way, these two facts had not previously been linked. Jones et al. studied body movements and brain activity in mice, and found that scalefree body movements were linked to scalefree brain activity, but only in certain subsets of neurons. This observation had been hidden because other subsets of neurons compete with scalefree neurons. When the scalefree neurons turn on, the competing groups turn off. When averaged together, these fluctuations cancel out.
The findings of Jones et al. provide a new understanding of how brain and body dynamics are orchestrated in healthy organisms. In particular, their results suggest that the complex, multiscale nature of behavior and body movements may emerge from brain activity operating at a critical tipping point between order and disorder, at the edge of chaos.
Introduction
From fidgeting to expeditions, natural body movements manifest over a very broad range of spatiotemporal scales. The complexity of such movements is often organized with fractal structure, scalefree fluctuations spanning multiple spatiotemporal orders of magnitude (Anteneodo and Chialvo, 2009; Proekt et al., 2012; Sims et al., 2008; Viswanathan et al., 1999). Such behavioral complexity may be beneficial for foraging (Garg and Kello, 2021; Sims et al., 2008; Viswanathan et al., 1999; Wosniack et al., 2017), visual search (Viswanathan et al., 1999), decisionmaking based on priority (Barabási, 2005; Sorribes et al., 2011), flexible switching of behavior (Abe, 2020), and perhaps more. Similarly, fluctuations of ongoing neural activity in the cerebral cortex can exhibit fractal, scalefree fluctuations like the spatiotemporal cascades sometimes referred to as ‘neuronal avalanches’ (Beggs and Plenz, 2003; Bellay et al., 2015; Ma et al., 2019; Priesemann et al., 2014; Scott et al., 2014; Shew et al., 2015; Shriki et al., 2013; Tagliazucchi et al., 2012; Yu et al., 2017) and longrange temporal correlations (Hardstone et al., 2012; Kello et al., 2010; Palva et al., 2013; Smit et al., 2013). Considering these observations of behavior and brain activity together, a simple, yet unanswered question arises. Are scalefree dynamics in cerebral cortex related to scalefree body movements? Or are these statistical similarities merely superficial, without any direct link?
Two previous studies in humans have shown that temporally scalefree ongoing brain dynamics are related to some subtle aspects of behavior – fluctuations in success at rhythmic finger tapping (Smit et al., 2013) and sensory detection tasks (Palva et al., 2013). But, here, we are concerned with less subtle, largerscale body movements and invasive measurements of brain activity with single neuron resolution.
According to prevailing theories of scalefree ongoing neural activity, these dynamics are often interpreted as ‘background’ activity, not directly linked to behavior. This view arises in part because computational models of scalefree neural activity generate this type of dynamics autonomously, due to internal interactions, without explicit modeling of a body that behaves (Dalla Porta and Copelli, 2019; di Santo et al., 2018; GirardiSchappo et al., 2020; Li and Shew, 2020; Wilting et al., 2018). Similarly, in experiments, scalefree neural activity has typically been observed in animals that are at rest, not engaged by any particular motor task or sensory input, or under anesthesia (Fontenele et al., 2019; Gautam et al., 2015; Gireesh and Plenz, 2008), or in vitro (Beggs and Plenz, 2003; Bellay et al., 2020; Shew et al., 2009; Tetzlaff et al., 2010) where there is no behavior to be observed. Based on these previous experimental observations and modeling efforts, one might naturally conclude that scalefree neural activity is internally generated, emerging spontaneously, without any link to behavior.
However, recent experiments in awake mice have revealed that behaviors typically ignored by experimenters, such as whisking, eye and face movement, fidgeting, walking, and running, are strongly related to ongoing neural activity in the cortex (Clancy et al., 2019; Musall et al., 2019; Salkoff et al., 2020; Stringer et al., 2019). These studies did not address our question here about scalefreeness of neural activity and behavior, but they suggest that if ongoing activity is scalefree, then it may be strongly related to behavior. More specifically, these studies point out the possibility that previous observations of scalefree neural activity in awake animals didn’t notice a link to behavior because those studies did not measure ongoing, spontaneous movements of the face, whiskers, and body. If this is the case, it would dramatically revise the typical interpretation of scalefree ongoing activity and provide direct evidence that scalefree behavior and scalefree brain activity are inseparably related. However, in the work by Stringer et al., 2019, the measured behaviors could explain only 10–30% of the total variance in neural activity fluctuations. Thus, the question remains: is this behaviorallyrelated 10–30% of neural activity scalefree? Here, our aim was to answer this question.
Results
Here, we analyzed simultaneous recordings of more than 10,000 pyramidal neurons during 1–2 hr in visual cortex of awake mice ( seven mice, nine recordings, recording duration 1.52 ± 0.51 hr), performed and first reported by Stringer et al., 2019. The neurons were distributed throughout a volume (300–400 μm indepth, 0.9 × 0.935 mm in the area) and sampled at 3 Hz. In addition to the neural activity, several aspects of behavior and body movements were also recorded (Figure 1). Here, we focused on four aspects of behavior. First, we studied run speed, which was assessed using an optical sensor and a floating spherical treadmill. Second, we examined pupil diameter, which was obtained from a camera targeting the faces of the mice. Pupil diameter dynamics are associated with changes in arousal and other aspects of body movement and sensory input (Reimer et al., 2016; Reimer et al., 2014; Vinck et al., 2015). Third, we studied changes in direction of gaze by tracking the speed of the center of the pupil as the mice looked around. Finally, we used previously developed methods to study whisker motion (Stringer et al., 2019), which was also obtained from the face camera.
Scalefree behavior
The behavior tended to occur in bursts; rest periods were punctuated with welldefined bouts of body movement (Figure 1C and Figure 2A). All four of the behavioral variables we studied tended to start and stop together, approximately, but differed in their detailed, fast fluctuations (Figure 1C). We first sought to determine whether the behaviors exhibited scalefree dynamics. To this end, we first defined each bout of elevated body activity, hereafter called a behavioral ‘event,’ based on a threshold (the median of the entire time series). Each behavioral event began when the behavioral time series, run speed for example (Figure 2A), exceeded the threshold and ended when it returned below the threshold. The ‘size’ of each behavioral event was defined as the area between the threshold and the variable during the event (Figure 2A). This definition of events and their sizes is motivated by previous studies of scalefree neural activity (Gautam et al., 2015; Larremore et al., 2014; Li and Shew, 2020; Poil et al., 2012), typically referred to as neuronal avalanches. In the case of run speed, the event size corresponds to an effective distance traveled during the event. If the animal behavior is scalefree, one expects the distribution Pr(s) of behavioral event sizes (s) to have a powerlaw form, Pr(s)~s^{τ}. We found that these distributions often were in the form of a powerlaw over a wide range of event sizes (Figure 2B). For example, the run speed event size distribution for mouse 1 was scalefree for nearly four orders of magnitude of event sizes (powerlaw range, r=3.9 decades). To rigorously assess the range over which the data is powerlaw distributed, we used a maximum likelihood fitting algorithm that accounts for the number of events observed and the possibility of a cutoff at the head and tail of the distribution (Methods). Our approach builds on that used in previous studies (Shew et al., 2015). We report the powerlaw range r in decades, i.e., the number of orders of magnitude. If no range meets our fit criteria for statistical significance, we report r=0.
Although not all behaviors were scalefree over such a large range for all recordings, we found that most of the mice exhibited scalefree behaviors, particularly for running and pupil fluctuations (Figure 2C). In addition to the powerlaw range for larger behavioral events, some behaviors exhibited a range of smallscale behavioral events with a flat distribution consistent with low amplitude noise in the measurement methods (Figure 2—figure supplement 1). These conclusions were robust to changes in the two most important parameters for the powerlaw fitting algorithm: the goodnessoffit criterion and the outlier exclusion threshold (more on outlier exclusion in Figure 2—figure supplement 2).
Scalefree neural activity
Next, we turned to the neural activity to test whether it was scalefree like the behavior. First, we averaged over all neurons to obtain a single population activity time series. Then, we treated the neural data in a similar way to the behavioral data. We defined neural events with a median threshold and defined event size as the area between the threshold and the data during the event (Figure 3A). This definition of events and event sizes is wellestablished in previous studies of scalefree neural activity (Gautam et al., 2015; Larremore et al., 2014; Li and Shew, 2020; Poil et al., 2012). (For more discussion of why this definition was chosen over other established alternatives, see Methods). The resulting distributions of neural event sizes were not scalefree; they were poorly fit by a powerlaw distribution (Figure 3B). Considering all 9 recordings, we found that r=1.5 ± 0.9 decades (mean ± SD). Such a short powerlaw range should be interpreted as evidence against any powerlaw at all, because nearly any function can be fit on a short range of data (more on why a large powerlaw range is important below).
One note of warning: at first glance, the distribution in Figure 3B may appear to have a decent powerlaw fit for the small event sizes. This is a misleading artifact resulting from using logarithmic bins and very few samples to estimate probability density. The inset in Figure 3B shows the number of events rather than probability density. See also Figure 2—figure supplement 1 for more clarification of this point. We emphasize that our powerlaw fitting algorithm does not use any bins and accounts for sample size, thus avoiding such misleading artifacts.
The lack of scalefree dynamics for the neural population was surprising considering two facts together: first, the behavior exhibited scalefree dynamics (Figure 2), and second, many neurons are strongly correlated with behavior (Figure 1 and Stringer et al., 2019). One possible explanation for this puzzling observation is that summing the entire population together obscures the scalefree dynamics of certain subsets of neurons. Next, we set out to test this possibility. Since we did not know, a priori, which neurons might be in these subsets, we adopted a bruteforce, shotgun search. First, we picked a ‘seed’ neuron at random from the entire population. Then we identified the 50 neurons that were most correlated with the seed neuron (we also tried 200 neurons without major impacts on the following results, Figure 3E). We averaged the activity of these 50 neurons to obtain a single time series and proceeded to define events and examine their size distribution. Two examples of time series obtained from subsets of neurons in this way are shown in Figure 3A. We repeated this process for 1000 different seed neurons. We found that some subsets of neurons were indeed scalefree, with powerlaw scaling over more than four orders of magnitude, while other subsets were not scalefree (Figure 3C). With so many neurons to choose from and a shotgun approach like this, it is important to avoid chancelevel false positive conclusions. As a conservative control for this possibility, we repeated the analysis, but with surrogate data, generated by applying a random time shift to each neuron’s activity time series relative to other neurons. This timeshifted control data, thus, has identical statistics to the original data at the single neuron level, but correlations among neurons occur only by chance. For timeshifted controls, we found that large powerlaw ranges (greater than 3.5 decades) were rare and not found at all in the majority of recordings (Figure 3D and E). This conclusion was robust for different choices of group size and goodnessoffit criteria used for the powerlaw fitting (Figure 3E).
Two concerns with the data analyzed here are that it was sampled at a slow time scale (3 Hz frame rate) and that the deconvolution methods used to obtain the data here from the raw GCAMP6s Ca imaging signals are likely to miss some of the underlying spike activity (Huang et al., 2021). Since our analysis of neural events hinges on summing up activity across neurons, could it be that the missed activity creates systematic biases in our observed event size statistics? To address this question, we analyzed some timeresolved spike data (Neuropixel recording from Stringer et al., 2019). Starting from the spike data, we created a slow signal, similar to that we analyzed here by convolving with a Catransient, down sampling, deconvolving, and zscoring (Figure 3—figure supplements 1–2). We compared neural event size distributions to ‘ground truth’ based on the original spike data (with no loss of spikes) and found that the neural event size distributions were very similar, with the same exponent and the same powerlaw range (Figure 3—figure supplement 1). Thus, we conclude that our reported neural event size distributions are reliable.
Linking scalefree brain activity and behavior
Having established that certain subsets of neurons have scalefree dynamics and that behavior has scalefree dynamics, we next sought to assess how the neural and behavioral dynamics are related. We first computed the correlation coefficient between each behavior time series and each neural subset time series. The correlation was very weak between the total neural population time series and behavior (red x’s in Figure 4). For the 50neuron subsets, we found a very wide range of behavioral correlation values, from 0.9 to –0.75 (Figure 4). Interestingly, the subsets with the greatest correlation with behavior also tended to have the largest ranges of powerlaw scaling. Moreover, the neural subsets with the smallest powerlaw range were most often near zero correlation with behavior. These results strongly suggest that the scalefree subsets of neurons are directly related to scalefree behavior. It is also interesting to note that many of the subsets that were strongly anticorrelated with behavior also had a large powerlaw range (the Ushaped relationships in Figure 4). These findings were very robust in 6 out of 9 recordings (mouse 1–4 c, Figure 4), but less prominent for some behaviors in the remaining three recordings (mouse 5–7). We note that these three ‘outlier’ animals also had tendencies to run much more than usual (mouse 6 and 7 ran 50 times more than the next most active animal) or much less than usual (mouse 5 ran 10 times less than the next least active animal) as shown in Figure 4—figure supplement 1. This observation suggests that our main findings hold best in vigilance states that are moderate, not too hyperactive, or hypoactive.
Powerlaw range, exponents, and scaling relations
Previous studies of scalefree neural activity have typically not reported powerlaw range. Readers may wonder why we focused our analyses here on the powerlaw range. So far, our results offer two answers to this question. First, and most importantly, only neural subsets with large powerlaw ranges were strongly correlated with behavior (Figure 4). Second, any claim of powerlaw scaling is stronger if that scaling extends over a greater range. Previous studies often rely on a ‘rule of thumb’ that one or two decades of scaling is not bad. For the data analyzed here, we found that a powerlaw range of 1–2 decades is insignificant compared to chance for event size distributions (Figure 3).
One reason previous studies of scalefree brain activity have not reported powerlaw range stems from their motivating hypothesis that scalefree activity results from the neural circuits operating in a special dynamical regime near the tipping point of a phase transition, i.e., near criticality. Rather than powerlaw range, these studies focused on scaling relationships between the powerlaw exponents of the event size and duration distributions as predicted by the theory of critical phenomena. Based on our results, we contend that these previous studies might be strengthened by also accounting for powerlaw range. Conversely, previous studies also raise a natural question about whether the data we analyzed here supports predictions from the theory of criticality. Our next goal is to answer this question.
According to prevailing theories of criticality in neural systems, both the sizes and durations of events should be powerlaw distributed. Moreover, the theory of crackling noise, which explains a broad class of stickslip critical phenomena (Sethna et al., 2001), predicts that the size powerlaw exponent τ and the duration powerlaw exponent α should be related to how size scales with duration. More specifically, S(T)~T^{β} where β_{crackling} = (α  1)/(τ –1). Many previous observations of scalefree neural activity support this prediction with values of β_{crackling} ≈ 1.1–1.3 (Fontenele et al., 2019; Fosque et al., 2021; Friedman et al., 2012; Ma et al., 2019; Shew et al., 2015). To test this theory, we first found the bestfit powerlaw and its exponent τ for the size distribution (Figure 5A). Then, we fit a powerlaw to the event duration distributions for the same set of events that were in the range of powerlaw scaling for event sizes (Figure 5B), thus obtaining α. And finally, we computed a bestfit β by doing a linear fit to log(S) versus log(T) as shown in Figure 5C. We found that τ, α, and β varied dramatically across different neural subsets (Figure 5D). However, if we examined only the behaviorallycorrelated subsets, i.e., those with a large powerlaw range (r>3.5), we found that the powerlaw exponents were relatively consistent: τ=1.0 ± 0.1, α=1.2 ± 0.2, and fit β=2.2 ± 0.2 (mean ± SD, Figure 5E). Thus, we conclude that powerlaw scaling exists for both sizes and durations of neural events and that size scales with duration according to another powerlaw. The existence of these scaling laws is consistent with criticality, but we found that our measured values of β were not consistent with crackling noise theory predictions; best fit β was not close to (α  1)/(τ –1). This was typically because τ was near one, resulting in small, sometimes negative values of τ–1 in the denominator of the predicted β formula. Moreover, our observations of β were greater than two, which is not close to previously reported values from experiments – around 1.1–1.3. This suggests that if the scalefree neural activity we observe has its origins in criticality, it may be a type of criticality that is different than previous observations. Indeed, this discrepancy with the existing theory is part of the motivation for the new model we present below. Additional scaling laws from the experiments are compared to the model results below.
Finally, we examined the exponents for distributions of behavioral events. Considering the behaviors with a large powerlaw range (r>3), we found τ=1.2 ± 0.2, α=1.4 ± 0.3, and fit β=1.8 ± 0.4. These exponents were related to the exponents and powerlaw ranges for neural event distributions, but only if we considered the neural subsets with large powerlaw range (r>3.5 decades) and strong correlation to behavior (R>0.6). We found that τ, α, r, and β were significantly correlated for neural and behavioral events (Pearson correlation R_{τ}=0.5, p<10^{–17}; R_{α}=0.2, p<0.002; R_{r} = 0.4, p<10^{–66}; R_{β}=0.3, p<10^{–4}). For these comparisons, we paired each neural subset with the type of behavior that was most correlated to that subset.
Correlations of individual behavioral and neural events
The calculation of correlations between neural subsets and behavior, like those in Figure 4, represents the entire recording duration, leaving open several questions. Are the high correlations in Figure 4 simply due to coarse, onoff dynamics of the neural and behavioral activity, or are there strong correlations at the more detailed level of the faster fluctuations during a single behavioral event? To answer these questions, we computed eventspecific correlations – computed between the behavior time series and the neural subset time series during each behavioral event, for each neural subset (Figure 6A). For instance, for 1000 neural subsets and 1000 behavioral events (which are typical numbers for this data), we would calculate 1 million correlation values. Each correlation is based on the relatively short time series between the start and stop of the corresponding behavioral event. These eventspecific correlations were often quite high (>0.9). However, the chancelevel occurrence of such high correlations can also be high for such short time series. We account for such chancelevel correlations by repeating the calculations for timeshifted control data (Methods), which defines the eventspecific, subsetspecific chancelevel occurrence rate.
First, we asked how many neural subsets were strongly correlated with each behavioral event and whether this count depended on behavioral event size. We found that many behavioral events had a significant number of such strongly correlated subsets and that the larger behavioral events tended to have a larger number of strongly correlated subsets (Figure 6B). Next, we asked how many behavioral events were strongly correlated with each neural subset and whether this count depended on the powerlaw range of the neural subset. We found that many subsets were strongly correlated with a significant number of behavioral events and that the subsets with greater powerlaw range tended to be correlated with more behavioral events (Figure 6C). The example results in Figure 6 are shown for all mice and all types of behaviors in Figure 6—figure supplements 1–2. Finally, we considered the time during which the mouse was active and asked what fraction of that time had a significant number of strongly correlated neural subsets. This fraction of time was often as high as 0.5 (Figure 6D). Thus, we conclude that the detailed, fast fluctuations of behavior are significantly and strongly correlated with scalefree neural activity.
Anticorrelated neural subsets
How is it possible for the total population to not exhibit scalefree fluctuations, while certain subsets are scalefree? A clue to this mystery is apparent upon close inspection of the results in Figure 4. Some neural subsets are strongly correlated with behavior, while other subsets are strongly anticorrelated, similar to previous reports from the motor cortex (Zagha et al., 2015) and prefrontal cortex (GarciaJuncoClemente et al., 2017). Since the scalefree subsets are often strongly correlated with behavior, it stands to reason that some neural subsets are strongly anticorrelated with the scalefree subsets. We show in Figure 7 that this fact explains why the total population does not exhibit scalefree dynamics. When two anticorrelated subsets are averaged together, they cancel out, resulting in relatively small fluctuations at the level of the total population (Figure 7A). The broad range of correlations and anticorrelations in each recording are more apparent using ‘correlation spectra’ (Figure 7B). For each previously defined seed neuron, we ranked all the other neurons according to their correlation with the seed neuron. These correlation spectra often reveal large numbers of strongly anticorrelated neurons. To directly show how cancelation of anticorrelated subsets abolishes scalefree dynamics, we recomputed the powerlaw range, but this time defined each subset to include the 25 most correlated and the 25 most anticorrelated neurons, relative to each seed neuron. The powerlaw range for these canceling subsets was greatly reduced compared to the original correlated subsets for all recordings (Figure 7C–E).
Stochastic winnertakeall competition and criticality
What mechanisms might be responsible for scalefree dynamics among a certain subset of neurons that are obscured due to cancelation with anticorrelated subsets of neurons? Theoretical models for studying scalefree dynamics in neural networks have often employed randomly connected (ErdosRenyi), networks of excitatory and inhibitory neurons, tuned to operate at criticality (Gautam et al., 2015; Li and Shew, 2020; Yang et al., 2012). However, such models do not manifest qualitatively different slow fluctuations in different subsets of neurons. More specifically, large groups of anticorrelated neurons do not occur in previous models. How can traditional models be changed to account for our observations? Here, we propose that a nonrandom network structure is needed. Indeed, experiments have shown that connectivity among cortical neurons is far from random (Cossell et al., 2015; Jiang et al., 2015; Song et al., 2005); this nonrandom structure is particularly pronounced among different types of inhibitory neurons (Jiang et al., 2015; Pfeffer et al., 2013; Pi et al., 2013; Wall et al., 2016). To our knowledge, this celltypespecific connectivity has been ignored in previous models of scalefree neural activity.
Which aspects of this network structure might explain our findings? We initially sought clues from wellknown circuit motifs that generate the anticorrelated activity underlying many types of animal locomotion (Kiehn, 2016). A common motif in these circuits, sometimes termed ‘crossing inhibition’ or ‘winnertakeall,’ entails two excitatory circuit nodes, say e+ and e−, whose activity is anticorrelated because they interact via at least one inhibitory node; e+ excites the inhibitory node, which suppresses e, and vice versa, e excites the inhibitory node which suppresses e+ (Figure 8A–D). Such circuit motifs certainly are common in the cortex, but they are mixed and interconnected with numerous other motifs. Even ErdosRenyi networks contain many motifs of this type. Is it plausible that ‘winnertakeall’ motifs are responsible for our findings?
As shown in Figure 8A, we tested this possibility using a computational model of 1000 binary neurons divided into two excitatory groups (e+ and e− with 400 neurons each) and two groups of inhibitory neurons (i+ and i− with 100 neurons each). We considered two inhibitory groups, instead of just one, to account for previous reports of anticorrelations between VIP and SOM inhibitory neurons in addition to anticorrelations between groups of excitatory neurons (GarciaJuncoClemente et al., 2017). The excitatory groups were densely connected within each group (50% connectivity) and sparsely connected across groups (5%), consistent with experimental observations (Cossell et al., 2015; Song et al., 2005). We also included dense crossing inhibition according to known connectivity among some inhibitory cell types (Jiang et al., 2015; Pfeffer et al., 2013; Pi et al., 2013; Wall et al., 2016). Excitatory and inhibitory synaptic weights were of the same absolute magnitude, but opposite in sign, and the connectivity matrix was normalized such that its largest eigenvalue was one (Figure 8B), following previous studies of scalefree dynamics and criticality (Gautam et al., 2015; Larremore et al., 2011; Li and Shew, 2020). Each neuron fired probabilistically, in proportion to the sum of its inputs (Methods).
We found that the e+ and e− neural subsets in this model produced large, slow, scalefree fluctuations (Figure 8C) like those observed in the scalefree subsets of neurons from our experimental observations (Figure 7A). The dynamics of e+ were strongly anticorrelated with those of e−, and i+ was anticorrelated with i− (Figure 8C), resulting in cancelation and a lack of scalefree dynamics when the entire population was considered together (Figure 8D). All neurons in the model were driven by unstructured random input; there was no statistical difference between the input to different subsets. Thus, the anticorrelated switching behavior in our model emerges stochastically. These results were robust to a wide range of different input rates (Figure 8E), provided that the largest eigenvalue of the connectivity matrix was near 1. If the largest eigenvalue was increased beyond one, the dynamics stopped switching between e+ and e−, instead getting locked into either e+ active and e− inactive or vice versa (this occurred to the right of the red dashed lines in Figure 8E and F). Thus, the stochastic switching between e+ and e− activity is a type of critical dynamics occurring at a boundary (red dashed line) in the model parameter space. Such critical winnertakeall dynamics never actually ‘choose a winner,’ instead switching randomly between the two ‘equally deserving winners.’ Scalefree anticorrelated fluctuations were also abolished when the density of connections between e+ and e− was increased beyond about 10% (Figure 8F). Above this density of connectivity between e+ and e−, the network effectively behaves like a single network, like previous models of criticality and scalefree dynamics.
In addition to testing dependence on the largest eigenvalue, the input rate, and the density of e+toe− connectivity, we did a systematic exploration of different network structures. We sought to understand if the network structure considered in Figure 8A–F is the only structure that reproduces our experimental observations, or if some other network structures could also work. Figure 8G summarizes the results of testing several other network structures. We found that crossing inhibition mediated by a single inhibitory population could suffice. The disinhibitory motif highlighted by other recent studies (GarciaJuncoClemente et al., 2017) is, thus, not necessary to reproduce our observations. Importantly, we showed that our observations can also occur in a more realistic network including a third population of inhibitory neurons with connectivity like that of PV neurons, according to known connectivity (Jiang et al., 2015; Pfeffer et al., 2013). We employed a brute force search considering 873,000 different combinations of dense, sparse, or zero connectivity for the 16 connections among the four network nodes in Figure 8A (more detail in Figure 8—figure supplement 1). We found only 31 configurations that matched the four criteria in Figure 8G. All matching configurations included one of two dense crossing inhibition motifs as well as segregated excitatory groups e+ and e− (Figure 8H and Figure 8—figure supplement 1). In fact, these two motifs were sufficient for predicting which configurations match our experimental results with 80% accuracy (Methods).
Finally, we considered the exponents and scaling relationships for the event size and duration distributions from the (example shown in Figure 9A, B and C). We considered the parts of model parameter space with nonzero powerlaw range for e+ (Figure 9D inset). We compared to experimental results from all neural subsets and all recordings (Figure 9E). For smaller powerlaw ranges, τ and α were more scattered and inconsistent. For the model, we found τ=1.3 ± 0.12, α=1.4 ± 0.15, and β=1.4 ± 0.04, for powerlaw range above three decades. For both the experiment and the model, the bestfit β was typically not well predicted by crackling noise theory; a confirmation of crackling noise theory would have the points falling on the unity line in the right panels of Figure 9D and E.
While our model offers a simple explanation of anticorrelated scalefree dynamics, its simplicity comes with limitations. Perhaps the most obvious limitation of our model is that it does not include neurons with weak correlations to both e+ and e− (those neurons in the middle of the correlation spectrum shown in Figure 7B). In Figure 8—figure supplement 2, we show that our model can be modified in a simple way to include asynchronous neurons. Another limitation is that we assumed that all nonzero synaptic connections were equal in weight. We loosen this assumption allowing for variable weights in Figure 8—figure supplement 2, without changing the basic features of anticorrelated scalefree fluctuations. Future work might improve our model further by accounting for neurons with intermediate correlations, i.e., less sharp boundaries of the e+ and e− groups.
Discussion
We have shown that ongoing, untrained locomotion, whisking, and pupil diameter changes in mice often exhibit scalefree dynamics. This scalefree behavior is directly related to concurrent scalefree cortical neural activity among certain subsets of neurons, with a strong onetoone correspondence between many behavioral and neural events. Moreover, the subsets of neurons with the greatest range of scalefreeness (i.e. the greatest powerlaw range) were the most correlated with behavior.
Why might a large powerlaw range among neurons be important for behavior? The behavioral variables we analyzed are likely to involve both local interactions among neighboring neurons and longrange coordination across many brain regions. For example, volitional movements like running and whisking presumably require coordinated interactions within and among motor, sensory, and other brain regions associated with decisions to move. One way to manifest such multiscale coordination from local to brainwide circuits is to have neural events with sizes that occur with diverse scales, i.e., neural events with large powerlaw range.
Previous studies suggest additional functional implications of this close relationship between scalefree behavior and brain activity. First, scalefree neural activity has been associated with multiple functional benefits for sensory information processing and memory (Clawson et al., 2017; Gautam et al., 2015; Shew and Plenz, 2013; Shriki and Yellin, 2016; Suárez et al., 2021), while scalefree behavior has been associated with benefits for foraging, search, and decision making (Abe, 2020; Barabási, 2005; Garg and Kello, 2021; Sims et al., 2008; Sorribes et al., 2011; Viswanathan et al., 1999; Wosniack et al., 2017). Our results suggest that these two lists of functional benefits, previously considered separately, may in fact be a single unified list that occurs together. This long list of benefits may explain why several types of neural plasticity have been shown to maintain scalefree neural activity (Hoseini et al., 2017; Ma et al., 2019; Shew et al., 2015).
What mechanisms are responsible for scalefree neural activity? The prevailing view is the criticality hypothesis (see for example Beggs and Plenz, 2003; Muñoz, 2018; Shew and Plenz, 2013; Wilting and Priesemann, 2019). In this view, scalefree dynamics occur because the system operates near criticality, i.e., near the tipping point between two distinct dynamical regimes  one synchronous, the other asynchronous. However, previously studied models of criticality fail to explain several crucial components of our observations. First, we observed that scalefree activity is limited to subsets of neurons, which were strongly anticorrelated with other subsets of neurons. This apparent stochastic winnertakeall competition between subsets results in the cancelation of fluctuations at the level of the total population, thus obscuring the scalefree dynamics of the subsets. Previous models also do not agree with the powerlaw scaling relationships we found – especially our measurements of β in the relation S(T)~T^{β}. Our results suggest that a new type of criticality may be required to explain these observations. One possibility suggested by our model is that the scalefree dynamics we observe occur at the boundary between winnerless switching and singlewinner lockedin dynamics (the red dashed line in Figure 8E and F). Additional theoretical efforts are necessary to more fully explore how the traditional criticality hypothesis relates to the competitive criticality suggested by our model.
At first glance, our results appear to not only contradict existing theory of criticality in neural systems, but also prior experimental observations based on different measurement methods. Most prominently, many previous studies have observed scalefree neural activity based on fast electrophysiological measurements using the total measured population, without accounting for anticorrelated subsets of neurons. Can fast ephys dynamics be reconciled with the slow calcium imaging signals we analyzed here? A detailed reconciliation of these different timescales of measurement will require more extensive investigation, but we provide a partial answer here (Figure 3—figure supplements 1–2). We analyzed an electrophysiological recording of spike times from 294 single units, also reported in the Stringer et al., 2019 paper. We found that when the activity is analyzed at a fast time scale (5 ms resolution), scalefree activity is clear using the full population and anticorrelations were very weak. The exponent of the power law was near –2, consistent with some previous reports based on timeresolved spike recordings in awake animals (Fontenele et al., 2019; Ma et al., 2019). In contrast, when the activity was analyzed at a slow time scale (330 ms resolution) like the calcium imaging data, strong anticorrelations emerged and scalefree dynamics were much more apparent (i.e. power law range was greater) when analyzed in subsets of neurons. This result (Figure 3—figure supplement 2) suggests that our findings are not inconsistent with prior studies of scalefree activity based on timeresolved electrophysiological data and calls for a more thorough study of how correlation structure and scalefree dynamics differ and relate across spatiotemporal scales. At the coarse time resolution of the calcium imaging data, we note that the neural subsets with large powerlaw range tended to be composed of neurons with high activity variance and high pairwise covariance (Figure 8—figure supplement 3). We note that a similar point about the emergence of anticorrelations at slow timescales was made in the original Stringer et al., 2019 paper.
What causes the scalefree fluctuations in behavior we observed from experiments? From the point of view of our model and the criticality hypothesis more generally, scalefreeness originates in the dynamics of neural circuits. Indeed, the model generates scalefree fluctuations without any interactions with an ‘outside world.’ These internally generated scalefree neural dynamics then, might manifest as scalefree behavior. However, recalling that the experimental measurements reported here were from primary visual cortex, one might ask whether the scalefree fluctuations we observed in V1 are driven by scalefree fluctuations in visual input. In our view, this is unlikely given several facts. First, the gaze motion often had the least convincing scalefreeness in our results (i.e. smallest powerlaw range, Figure 2C) and gaze motion would be the most obvious reason for changes in visual input (there were no dynamic visual stimuli on a viewed screen for these recordings). Second, the Stringer et al., paper (2019), from which our data originated, showed that spatially distributed electrophysiological recordings from multiple brain regions (not just sensory) were correlated to behavior in a similar way. This suggests that our results here may be a brainwide phenomenon, not just a V1 phenomenon.
Setting aside speculations, our results firmly establish that complex fluctuations and events that make up scalefree neural activity are not ‘background noise’ nor ‘internal’ cognitive processes as is often supposed. Our findings call for a revision of these traditional views; scalefree neural activity directly causes (or perhaps is caused by) scalefree behavior and body motion.
Materials and methods
Animals and data acquisition
Request a detailed protocolThe data analyzed here were first published in Stringer et al., 2019 and are publicly available at doi:10.25378/janelia.6163622.v4.
All animal protocols and data acquisition methods are described in the original publication. All experimental procedures were conducted according to the UK Animals Scientific Procedures Act (1986). Experiments were performed at University College London under personal and project licenses released by the Home Office following appropriate ethics review. In brief, multiplane calcium imaging during periods with no visual stimulation was performed in seven adult mice (P35 to P125) bred to express GCaMP6s in excitatory neurons. Our labeling of the different recordings corresponds to the original authors’ labeling system of the data as follows:
mouse 1  spont_M150824_MP019_20160405,
mouse 2  spont_M160825_MP027_20161212,
mouse 3  spont_M160907_MP028_20160926,
mouse 4 a  spont_M161025_MP030_20161120,
mouse 4b  spont_M161025_MP030_20170616,
mouse 4 c  spont_M161025_MP030_20170623,
mouse 5  spont_M170714_MP032_20170804,
mouse 6  spont_M170717_MP033_20170818,
mouse 7  spont_M170717_MP034_20170825.
Neural data preprocessing
Request a detailed protocolBeginning from the deconvolved fluorescence traces (variable called Fsp in the original data set), first, the data were zscored. That is, for each neuron, we subtracted its timeaveraged activity and divided by its standard deviation. Second, we applied a lowpass filter (cutoff frequency 0.2 Hz, 2nd order butterworth filter, using Matlab filtfilt function).
Definition of neural events
Request a detailed protocolFrom previous studies of scalefree neural activity, there are three established strategies for defining neural events. The original approach was to first bin time (with time bin durations usually defined by the average interspikeinterval, or average interval between LFP peaks) and define an event or avalanche as a sequence of time bins with at least one spike or LFP peak, preceded and terminated by an empty time bin (Beggs and Plenz, 2003). This approach has been very useful when the data come from a relatively small number of neurons or electrodes but runs into difficulty with very large numbers of single neurons. The simple reason for this difficulty is that for very large numbers of neurons, there is simply no empty time bins (down to the time resolution of measurements), resulting in a single event that never ends. A second approach that better accounts for the possibility of more than one event occurring at the same time came from analyzing fMRI brain activity (Tagliazucchi et al., 2012). This approach assumes that an event must evolve in a spatially contiguous way as it spreads through the brain. This approach, while good for spatially coarse measurements like fMRI, does not make sense for a small local circuit with a single neuron resolution like that analyzed here. The third approach is the one we adopted here, described in our Results section, which is the only good option, to our knowledge, for single neuron resolution and very large numbers of neurons.
Powerlaw fitting and range
Request a detailed protocolTo assess the scalefreeness of behavioral events and neural population events, we developed an algorithm for finding the range of event sizes that are well fit by a power law. The algorithm is based on a maximum likelihood fitting procedure as established in previous work (e.g. Clauset et al., 2009; Langlois et al., 2014; Shew et al., 2015). We fit the measured event size distribution with a truncated powerlaw with minimum event size s_{min} and maximum event size s_{max}, excluding data that fell outside these bounds during the fit process. Similar to previous studies, our algorithm has two fitting parameters  s_{min} and the powerlaw exponent. s_{max} was not a fitting parameter; it was chosen to be the largest observed event size (for subsets with powerlaw range >3, the mean ± SD for s_{max} was 59 ± 26 s for event durations and 102 ± 48 for event sizes). We tested exponents between 0.7 and 2 in steps of 0.02. We tested s_{min} values between the smallest observed event size and the largest observed event size, increasing in 10 logarithmically spaced increments per decade. Thus, the powerlaw range reported in the manuscript has a resolution of 0.1 decades. Importantly, our fitting algorithm is entirely independent of any choice of binning that might be used to visualize the event size distribution.
For our purposes of measuring power law range, previously reported algorithms required improvements to account for (exclude) confounding outlier event sizes. These ‘outliers’ were rare events that caused noise in the extremes of the distribution tail or head and, in some cases, caused spuriously large powerlaw range estimates. Accounting for these outliers provides a more conservative estimate of powerlaw range. We defined outliers by first ranking all event sizes from smallest to largest. Next, we computed differences in sizes (in decades) for consecutive sizes in this ranked list. Outliers have a large difference in size compared to the following size (or preceding size). We define outliers as events with a size difference greater than 3% of the total range in decades. We also tried 6% as an outlier threshold and found that our results were not very sensitive to this choice (Figure 2C).
The fitting algorithm executed the following steps. First, outliers were excluded. Second, events with size less than s_{min} were excluded. Third, the maximum likelihood powerlaw exponent was calculated. Fourth, we assessed the goodnessoffit. We repeated these four steps for all the possible s_{min} values, in the end, identifying the largest range (smallest s_{min}) that passed our goodnessoffit criterion.
The goodnessoffit criterion we adopted here was also new, to our knowledge, and better suited to our goals of assessing powerlaw range, compared to previous methods. The steps for quantifying goodnessoffit were as follows. First, we created a cumulative distribution function (CDF) of the real data. Second, we created 500 surrogate data sets drawn from the bestfit truncated powerlaw. Third, we created 500 CDFs, one for each of the surrogate data sets. Fourth, we resampled all 501 CDFs with 10 logarithmically spaced points per decade, linearly interpolated. Fifth, we calculated the fraction F of points in the resampled CDF of the real data that fell within the bounds of the 500 resampled surrogate CDFs. F is our goodnessoffit measure. F=1 means that the entire range of the real data falls within the expected variation for a perfect powerlaw, given the number of samples in the data. We explored various goodnessoffit criteria between F=0.75 and F=0.9 for the behavioral data (Figure 2C) and the neural data (Figure 3E). This is a rather conservative goodnessoffit test, much more conservative than the more typically used KolmogorovSmirnov statistic for example.
The 500 surrogate data sets were also used for the purpose of creating the grayshaded regions in each plotted distribution of avalanche sizes and durations (Figure 2C, Figure 5A and B, Figure 9A and B). First, each surrogate data set was used to create a probability density function (PDF) with the same bins as the real data. Each of these PDFs was slightly different due to statistical noise (i.e. finite sample size). The grayshaded region represents the range between which 90% of the surrogated PDFs fell. Thus, the gray region indicates how much variability one should expect for a perfect power law with the best fit exponent and the same range and sample size as the real data. We also use the upper and lower bounds of the gray region to estimate error bars for our exponents (as reported in Figure 2C, Figure 5A and B, Figure 9A and B) We calculate τ_{up} from the slope of the upper boundary of the gray region and τ_{lo} from the slope of the lower boundary of the gray region. Then we report τ = τ_{MLE} ± τ_{VAR}, where τ_{MLE} is the maximum likelihood exponent and τ_{VAR} = (τ_{up}  τ_{lo})/2. Considering the event size exponents for all neural subsets with powerlaw range greater than 3, we found that τ_{VAR}=0.07 ± 0.05 (mean, SD); corresponding event duration exponents had τ_{VAR}=0.11 ± 0.07. For all behavioral event size distributions, τ_{VAR}=0.19 ± 0.17.
In Figure 2—figure supplement 2, we report the results of a benchmarking analysis in which we ran our powerlaw fitting algorithm on synthetic data sets drawn from known powerlaw distributions. We found that the exponents and powerlaw ranges found by our fitting algorithm are robust for a reasonable range of outlier exclusion thresholds (including 3% and 6%), different exponents, and sample sizes.
Statistical significance of correlations
Request a detailed protocolFigure 6 is based on computing correlation coefficients of short time series of behavior and neural activity. The number of such correlation coefficients computed for each behavioral event was 1000 because there were 1000 neural subsets. In the manuscript, we concluded that a significant number of these correlations were strong. This conclusion was arrived at by comparing to timeshifted control neural data (same as described for the timeshifted controls in the manuscript for Figure 3). The time shifts were randomly chosen from the interval 0 s up to the entire duration of the recording. A cyclic permutation was performed such that a time shift of + N samples would move the ith sample to become the (i+N)^{th} and the last N samples become the first N. For each behavioral event, we obtained 1000 timeshifted control correlation coefficients. We then counted how many of the real correlation coefficients were greater than 999 of the timeshifted control correlation coefficients. By chance, we should expect this count to be one. We concluded that there was a significant number of strong correlations if this count was greater than four. A similar reasoning was used to conclude that there were a significant number of strongly correlated behavioral events for each neural subset.
Computational model
Request a detailed protocolThe model consisted of N=1000 binary neurons; the state of the ith neuron at time t is ${s}_{i}\left(t\right)=0$ (quiescent) or ${s}_{i}\left(t\right)=1$ (firing). The population was divided into four groups called e+, e−, i+, and i−, as diagramed in Figure 8A. The e+ and e− groups include 400 excitatory neurons each; i+ and i− include 100 inhibitory neurons each. The dynamics of the neurons are updated synchronously according to
where the activation probability is
Here, $\eta $ is a constant representing input from outside the model network and/or a tendency to fire spontaneously without input. We characterized how changes in $\eta $ affect model dynamics in Figure 8E (vertical axis, labeled ‘external drive’). The connection from neuron j to neuron i is given by ${C}_{ij}$ , one element of the connection matrix depicted in Figure 8A. Certain pairs of groups were sparsely connected, while other pairs of groups were densely connected, or not connected at all. 50% of connections (randomly chosen) were zero for two densely connected groups. 95% of connections (randomly chosen) were zero for a sparsely connected pair of groups. All nonzero excitatory connections were set to the same constant c; all nonzero inhibitory connections were set to c. The value of c was set by first normalizing the entire connection matrix by its largest eigenvalue and then multiplying by another constant $\Lambda $ to obtain the desired largest eigenvalue. We studied how the model dynamics depend on $\Lambda $ in Figure 8F and Figure 8G (horizontal axis). For the shotgun search of many different circuit configurations, we set $\mathrm{\Lambda}=1$. We also studied how the density of connectivity between e+ and e− impacted model dynamics in Figure 8F (vertical axis, called ‘crossing excitation’).
All model data analysis was done on the population activity and averaged over neurons for each group.
We did a shotgun search for model configurations that result in dynamics consistent with our experimental results. To do this, we first generated a list of all possible configurations of the 16 pairwise group connections. We considered three possible densities for each of these connections: disconnected (0%), sparse (5%), or dense (50%). Thus, the list of all possible connection configurations includes 3^{16} (more than 43 million) possibilities. Clearly testing all these possibilities would take a long time, so we pared down the list with the following constraints, which exclude several very unrealistic possibilities:
Both e+ and e− groups must have nonzero withingroup connectivity.
No disconnected components were allowed. Every group had to be reachable by at least one connection.
There must be at least one excitatory to inhibitory connection.
There must be at least one inhibitory to excitatory connection.
After applying these constraints, the list of possible configurations still included 18,576,000 possible circuits. We tested 8,73,000 of these possibilities, chosen at random. We considered a circuit to be consistent with experimental results if it met the following conditions:
Correlation coefficient of e+ and e− must be less than –0.5
Correlation coefficient of i+ and i− must be less than –0.5
Power law range of e+ or e− must be greater than 3.5 decades
Power law range of the total population must be less than 2 decades
We found that only 31 circuit configurations met these conditions. These matching networks are illustrated in Figure 8—figure supplement 1. As noted in the main text, all matching configurations featured two structural motifs. First, they included dense selfconnections within e+ and within e−, but weak or nonexistent connections between e+ and e−. Second, all matching configurations included one of two dense crossing inhibition configurations (Figure 8H). Moreover, we found that these two conditions were sufficient to predict which networks would match the experimental results with 80% accuracy if we loosened the criteria for matching experiments slightly. These looser criteria were:
Correlation coefficient of e+ and e− must be less than –0.5
Correlation coefficient of i+ and i− must be less than –0.2
Power law range of e+ or e− must be at least 1.4 times that of the total population
Data availability
The data analyzed here were first published in Stringer et al., 2019 and are publicly available on figshare at https://doi.org/10.25378/janelia.6163622.v6. Analysis code is publicly available on figshare at https://doi.org/10.1101/2021.05.12.443799.

figshareCode for computational model, power law fitting, power law range [Shew Lab].https://doi.org/10.6084/m9.figshare.21954389.v1

FigshareRecordings of ten thousand neurons in visual cortex during spontaneous behaviors.https://doi.org/10.25378/janelia.6163622.v6
References

Neuronal avalanches in neocortical circuitsThe Journal of Neuroscience 23:11167–11177.https://doi.org/10.1523/JNEUROSCI.233511167.2003

Selective participation of single cortical neurons in neuronal avalanchesFrontiers in Neural Circuits 14:620052.https://doi.org/10.3389/fncir.2020.620052

Locomotiondependent remapping of distributed cortical networksNature Neuroscience 22:778–786.https://doi.org/10.1038/s4159301903578

Criticality between cortical statesPhysical Review Letters 122:208101.https://doi.org/10.1103/PhysRevLett.122.208101

Evidence for quasicritical brain dynamicsPhysical Review Letters 126:98101.https://doi.org/10.1103/PhysRevLett.126.098101

Universal critical dynamics in high resolution neuronal avalanche dataPhysical Review Letters 108:1–5.https://doi.org/10.1103/PhysRevLett.108.208102

Fast online deconvolution of calcium imaging dataPLOS Computational Biology 13:e1005423.https://doi.org/10.1371/journal.pcbi.1005423

An inhibitory pullpush circuit in frontal cortexNature Neuroscience 20:389–392.https://doi.org/10.1038/nn.4483

Efficient Lévy walks in virtual human foragingScientific Reports 11:5242.https://doi.org/10.1038/s4159802184542w

Maximizing sensory dynamic range by tuning the cortical state to criticalityPLOS Computational Biology 11:e1004576.https://doi.org/10.1371/journal.pcbi.1004576

Synaptic balance due to homeostatically selforganized quasicritical dynamicsPhysical Review Research 2:12042.https://doi.org/10.1103/PhysRevResearch.2.012042

Detrended fluctuation analysis: a scalefree view on neuronal oscillationsFrontiers in Physiology 3:450.https://doi.org/10.3389/fphys.2012.00450

Scaling laws in cognitive sciencesTrends in Cognitive Sciences 14:223–232.https://doi.org/10.1016/j.tics.2010.02.005

Decoding the organization of spinal circuits that control locomotionNature Reviews. Neuroscience 17:224–238.https://doi.org/10.1038/nrn.2016.9

Maximum likelihood estimators for truncated and censored powerlaw distributions show how neuronal avalanches may be misevaluatedPhysical Review. E, Statistical, Nonlinear, and Soft Matter Physics 89:012709.https://doi.org/10.1103/PhysRevE.89.012709

Inhibition causes ceaseless dynamics in networks of excitable nodesPhysical Review Letters 112:138103.https://doi.org/10.1103/PhysRevLett.112.138103

Tuning network dynamics from criticality to an asynchronous statePLOS Computational Biology 16:e1008268.https://doi.org/10.1371/journal.pcbi.1008268

Colloquium: criticality and dynamical scaling in living systemsReviews of Modern Physics 90:31001.https://doi.org/10.1093/ejcts/ezx068

Singletrial neural dynamics are dominated by richly varied movementsNature Neuroscience 22:1677–1686.https://doi.org/10.1038/s4159301905024

Spike avalanches in vivo suggest a driven, slightly subcritical brain stateFrontiers in Systems Neuroscience 8:108.https://doi.org/10.3389/fnsys.2014.00108

Tuning to odor solubility and sorption pattern in olfactory epithelial responsesThe Journal of Neuroscience 34:2025–2036.https://doi.org/10.1523/JNEUROSCI.373613.2014

Neuronal avalanches imply maximum dynamic range in cortical networks at criticalityThe Journal of Neuroscience 29:15595–15600.https://doi.org/10.1523/JNEUROSCI.386409.2009

The functional benefits of criticality in the cortexThe Neuroscientist 19:88–100.https://doi.org/10.1177/1073858412445487

Adaptation to sensory input tunes visual cortex to criticalityNature Physics 11:659–663.https://doi.org/10.1038/nphys3370

Neuronal avalanches in the resting MEG of the human brainThe Journal of Neuroscience 33:7079–7090.https://doi.org/10.1523/JNEUROSCI.428612.2013

Optimal information representation and criticality in an adaptive sensory recurrent neuronal networkPLOS Computational Biology 12:e1004698.https://doi.org/10.1371/journal.pcbi.1004698

Longrange temporal correlations in restingstate α oscillations predict human timingerror dynamicsThe Journal of Neuroscience 33:11212–11220.https://doi.org/10.1523/JNEUROSCI.281612.2013

The origin of behavioral bursts in decisionmaking circuitryPLOS Computational Biology 7:e1002075.https://doi.org/10.1371/journal.pcbi.1002075

Learning function from structure in neuromorphic networksNature Machine Intelligence 3:771–786.https://doi.org/10.1038/s42256021003761

Selforganized criticality in developing neuronal networksPLOS Computational Biology 6:e1001013.https://doi.org/10.1371/journal.pcbi.1001013

BrainWide maps of synaptic input to cortical interneuronsThe Journal of Neuroscience 36:4000–4009.https://doi.org/10.1523/JNEUROSCI.396715.2016

Operating in a reverberating regime enables rapid tuning of network states to task requirementsFrontiers in Systems Neuroscience 12:55.https://doi.org/10.3389/fnsys.2018.00055

25 years of criticality in neuroscienceestablished results, open controversies, novel conceptsCurrent Opinion in Neurobiology 58:105–111.https://doi.org/10.1016/j.conb.2019.08.002

The evolutionary origins of Lévy walk foragingPLOS Computational Biology 13:e1005774.https://doi.org/10.1371/journal.pcbi.1005774

Maximal variability of phase synchrony in cortical networks with neuronal avalanchesThe Journal of Neuroscience 32:1061–1072.https://doi.org/10.1523/JNEUROSCI.277111.2012
Decision letter

Gordon J BermanReviewing Editor; Emory University, United States

Timothy E BehrensSenior Editor; University of Oxford, United Kingdom

Keith B HengenReviewer; Washington University, United States
Our editorial process produces two outputs: (i) public reviews designed to be posted alongside the preprint for the benefit of readers; (ii) feedback on the manuscript for the authors, including requests for revisions, shown below. We also include an acceptance summary that explains what the editors found interesting or important about the work.
Decision letter after peer review:
Thank you for submitting your article "Scalefree behavioral dynamics directly linked with scalefree cortical dynamics" for consideration by eLife. Your article has been reviewed by 3 peer reviewers, and the evaluation has been overseen by a Reviewing Editor and Ronald Calabrese as the Senior Editor. The following individual involved in review of your submission has agreed to reveal their identity: Keith B Hengen (Reviewer #1).
The reviewers have discussed their reviews with one another, and although there were some concerns, there was general enthusiasm for the approach and for the ideas presented in the work. Accordingly, the Reviewing Editor has drafted this to help you prepare a revised submission.
Essential revisions:
1. The reviewers would like to see a detailed and thoughtful reflection on the role that 3 Hz Ca imaging might play in the conclusions that the authors derive. While the dataset in question offers many neurons, this approach is, from other perspectives, impoverished – calcium intrinsically misses spikes, a 3 Hz sampling rate is two orders of magnitude slower than an action potential, and the recordings are relatively short for amassing substantial observations of low probability (large) avalanches. The potential concern is that some of this disconnect may reflect optophysiological constraints. One argument against this is that a truly scale free system should be observable at any temporal or spatial scale and still give rise to the same sets of power laws. This quickly falls apart when applied to biological systems which are neither infinite in time nor space. As a result, the severe mismatch between the spatial resolution (single cell) and the temporal resolution (3 Hz) of the dataset, combined with filtering intrinsic to calcium imaging, raises the possibility that the conclusions are influenced by the methods.
2. Another reservation expressed by the referees has to do with the generality of the conclusions drawn from the mechanistic model. One of the connectivity motifs identified appears to be i+ to e and i to e+, where potentially i+/i are SOM and VIP (or really any specific inhibitory type) cells. The specific connections to subsets of excitatory cells appear to be important (based on the solid lines in Figure 8). This seems surprising: is there any experimental support for excitatory cells to preferentially receive inhibition from either SOM or VIP, but not both? More broadly, there was concern that the neat diagrams drawn here are misleading. The sample raster, showing what appears to be the full simulation, certainly captures the correlated/anticorrelated pattern of the 100 cells most correlated with a seed cell and 100 cells most anticorrelated with it, but it does not contain the 11,000 cells in between with zero to moderate levels of correlation. We probably expect that the full covariance matrix has similar structure from any seed (see Meshulam et al. 2019, PRL, for an analysis of scaling of coarsegrained activity covariance), and this suggests multiple crossover inhibition constraints, which seem like they could be hard to satisfy. The motifs identified in Figure 8 likely exist, but I am left with many questions of what we learned about connectivity rules that would account for the full distribution of correlations. Would starting with an ErdosRenyi network with slight overrepresentation of these motifs be sufficient? How important is the homogeneous connection weights from each pool assumption – would allowing connection weights with some dispersion change the results?
3. Putting 2) another way, it's unclear why the averaging is required in the first place. This operation projects the entire population down in an incredibly lossy way and removes much of the complexity of the population activity. Second, the authors state that it is highly curious that subsets of the population exhibit power laws while the entire population does not. While the discussion and hypothesizing about different ei interactions is interesting, it is possible that there's a discussion to be had on a much more basic level of whether there are topology independent explanations, such as basic distributions of correlations between neurons that can explain the subnetwork averaging. Specifically, if the correlation to any given neuron falls off, e.g., with an exponential falloff (i.e., a Gaussian Process type covariance between neurons), it seems that similar effects should hold. This type of effect can be easily tested by generating null distributions existing code bases. This is an important point since local (broadly defined) correlations of neurons implying the observed subnetwork behavior means that many mechanisms that have local correlations but don't cluster in any meaningful way could also be responsible for the local averaging effect.
4. In general, the discussion of "two networks" seems like it relies on the correlation plot of Figure~7B. The decay away from the peak correlation is sharp, but there does not seem to be significant clustering in the anticorrelation population, instead a very slow decay away from zero. The authors do not show evidence of clustering in the neurons, nor any biophysical reason why e and i neurons are present in the imaging data. The alternative explanation (as mentioned in (b)) is that the there is a more continuous set of correlations among the neurons with the same result. In fact one of the reviewers tested this themself using code to generate some data with the desired statistics, and the distribution of events seems to also describe this same observation. Obviously, the full test would need to use the same event identification code, and so it is quite important that the authors consider the much more generic explanation for the subnetwork averaging effect. We recommend assessing the possibility that broader explanations, e.g., in the form of the distributions of correlations accounts for the observed phenomenon. Even with 10K neurons, there are many other forces at play influencing the observed network and while it is nice that ei networks are one explanation, much less constraining explanations that are still biophysically feasible should be discussed and compared against. I have provided one possible approach (see PDF for code and example figure: https://submit.elifesciences.org/eLife_files/2022/05/19/00107349/00/107349_0_attach_6_474_convrt.pdf).
5. Another important aspect here is how single neurons behave. It was not clear if single neurons were stated to exhibit a power law. If they do, then that would help in that there are different limiting behaviors to the averaging that pass through the observed stated numbers. If not, then there is an additional oddity that one must average neurons at all to obtain a power law. We recommend the authors show the full curve so that the readers get a more detailed sense of how averaging effects the powerlaw interpretation of the data.
6. There is something that seems off about the range of \β values inferred with the ranges of \tau and $\α$. With \tau in [0.9,1.1], then the denominator 1\tau is in [0.1, 0.1], which the authors state means that \β (found to be in [2,2.4]) is not near \β_{crackling} = (\α1)/(1\tau). It seems as this is the opposite, as the possible values of the \β_{crackling} is huge due to the denominator, and so \β is in the range of possible \β_{crackling} almost vacuously. Was this statement just poorly worded?
7. It is not clear if there is more to what the authors are trying to say with the specifics of the scale free fits for behavior. Apparently, these results are used to motivate the neural studies, but aside from that, the details of those ranges don't seem to come up again. Given that the primary connection between neuronal and behavioral activity seems to be Figure 4. The distribution of points in these plots seem to be very lopsided, in that some plots have large ranges of fewtono data points. It would be very helpful to get a sense of the distribution of points that are a bit hard to see given the overlapping points and superimposed lines. We recommend that the authors add distribution information to the plots in Figure 4B to give a sense of how points are spread through the [correlation with behavior]by[power law range] space. Potential plots might be a colocated histogram, or perhaps an uncertainty estimate as a function of correlation based on the number of points and variance. This would help show significance of the curves in a way that accounts for the uneven spread of datapoints.
8. Neural activity correlated with some behavior variables can sometimes be the most active subset of neurons. This could potentially skew the maximum sizes of events and give behaviorally correlated subsets an unfair advantage in terms of the scalefree range. In a similar vain to 8), what are the typical dynamic ranges for subsets correlated and uncorrelated with behavior? We recommend showing a number of these to see if those dynamic ranges are impacting the possible ranges in the [correlation with behavior]by[power law range] plots. Perhaps something like curve in each plot showing the minimum maximum value of the power law range per correlation range. In general, the reviewers struggled with the interpretation of Figure 4b in the sense that there seems to be such variability between mice. How much do the authors feel that this is a difference in neural populations imaged, vs changes in imaging conditions (illumination, window clarity, optical alignment) or differences in mouse activity levels?
Reviewer #1 (Recommendations for the authors):
This paper feels highly polished and thorough in its presentation. I truly enjoyed reading it and believe it will be of value to the community.
A nuanced question: mouse #5 and mouse #6 consistently seems to break the rules implied by the rest of the dataset. The behavioral descriptors look normal in figure 2, but neural structure is notably different from the rest of the group in 3E, 4B, 6C, 7E, and S4 (that one is behavior). Is there anything meaningfully different about these recordings (n cells, mean event rate, anatomical location etc)?
Reviewer #2 (Recommendations for the authors):
Comments, questions, and technical issues for the authors.
1. The values of tau (size exponent) from neural avalanches are surprisingly low. For instance, in the Ma et al. (Hengen) 2019 Neuron paper, tau ranged from 1.5 to 1.9, sometimes more than 2, but never less than 1. I am wondering how the imaging preprocessing affects the estimate of tau. For instance, if you started with simulated spiking activity that has avalanches with size pdf that go as a power law with exponent tau, and then use a forward model to generate "calcium imaging," and then apply deconvolution, zscoring, and lowpass filtering, and then measure the avalanches again: what is the new tau?
2. This may be an ignorant question (apologies). The power law range is quantified in decibels (dB) throughout the paper; do you actually mean decades?
3. Related to the setup of the model, was there a reason that there are no adaptation mechanisms in this network model, as there often are in mechanistic models for avalanche criticality (including past work by the authors of this paper, e.g. Shew's 2015 Nature Physics paper)? Also, there appears to have been an error with the reference manager, as this reference shows up twice in the reference list.
4. It would be helpful the authors could elaborate on predictions that their results make for future studies. Maybe this is rather technical, but can you tell us when you expect to find a powerlaw distribution as a function of how much of the population is sampled and for how long? What if you were analyzing Neuropixels data, where you lose the extensive spatial sampling (and the restriction to pyramidal cells only) but you gain 3 orders of magnitude in temporal resolution?
5. On page 8, you ask "are all behavioral events equally correlated to their concurrent neural events, or are certain neural events from certain subsets of neurons more strongly related to behavioral events?" I don't understand the question. What does it mean for all behavioral events to be equally correlated to concurrent neural events? Aren't "concurrent neural events" in specific subsets of neurons?
Reviewer #3 (Recommendations for the authors):
1. Limits of calcium imaging: My recommendation is to assess mathematically the potential impact of missing data on the range and powerlaw slope estimates, which are the primary values used throughout the paper.
2. Correlations and powerlaws in subsets.
2ac. My recommendation is to assess the possibility that broader explanations, e.g., in the form of the distributions of correlations accounts for the observed phenomenon. Even with 10K neurons, there are many other forces at play influencing the observed network and while it is nice that ei networks are one explanation, much less constraining explanations that are still biophysically feasible should be discussed and compared against. I have provided one possible approach (see PDF for code and example figure: https://submit.elifesciences.org/eLife_files/2022/05/19/00107349/00/107349_0_attach_6_474_convrt.pdf) that I hope will be useful to the authors.
2d. I recommend the authors show the full curve so that the readers get a more detailed sense of how averaging effects the powerlaw interpretation of the data.
3. Please check that this calculation and interpretation is correct.
4. Connection between brain and behavior:
4b. I recommend that the authors add distribution information to the plots in Figure~4B to give a sense of how points are spread through the [correlation with behavior]by[power law range] space. Potential plots might be a colocated histogram, or perhaps an uncertainty estimate as a function of correlation based on the number of points and variance. This would help show significance of the curves in a way that accounts for the uneven spread of datapoints.
4c. In a similar vein, what are the typical dynamic ranges for subsets correlated and uncorrelated with behavior? I recommend showing a number of these to see if those dynamic ranges are impacting the possible ranges in the [correlation with behavior]by[power law range] plots. Perhaps something like curve in each plot showing the minimum maximum value of the power law range per correlation range.
4d. In general I'm struggling with the interpretation of Figure~4b in the sense that there seems to be such variability between mice. How much do the authors feel that this is a difference in neural populations imaged, vs changes in imaging conditions (illumination, window clarity, optical alignment) or differences in mouse activity levels?
[Editors’ note: further revisions were suggested prior to acceptance, as described below.]
Thank you for resubmitting your work entitled "Scalefree behavioral dynamics directly linked with scalefree cortical dynamics" for further consideration by eLife. Your revised article has been evaluated by Timothy Behrens (Senior Editor) and a Reviewing Editor.
The manuscript has been improved, and the reviewers concurred that an eventual acceptance is likely, but there are some remaining issues that need to be addressed. Specifically, the reviewers think that it's important that you address the points raised by Reviewer #3 (see below) regarding the "Mechanisms vs. Statistics" questions.
Reviewer #1 (Recommendations for the authors):
The authors have responded effectively to previous reviews both in the updated writing (discussion and intro) as well as the models and results.
Reviewer #2 (Recommendations for the authors):
Overall the authors addressed my concerns adequately. The paper has important results and makes a valuable contribution.
Reviewer #3 (Recommendations for the authors):
I appreciate your time and effort to respond to my review points. I believe the majority of my points have been addressed. The main weakness I still see is the lack of discussion about the broader mechanisms beyond the basic structures described that could account for the observations (see below).
Mechanisms vs. Statistics
I would like to clarify my reason for bringing up this "statistical" viewpoint that I believe may have been lost in translation. The paper as it stands makes the following logical steps (in terms of the mechanistic model) 1) The data exhibits power law scaling under certain binning of neurons (effect E) and 2) One way to account for this effect E is to consider certain ei models. This is reasonable, but the overall search for possible mechanisms seems to be a combination of intuition and trial and error. Different architectures were tried and either "passed" or "failed".
The purpose of bringing up the statistical properties needed was to hopefully raise a conversation around what core properties are needed to replicate this effect. As activity statistics and connectivity are intimately related in mechanistic models, such a characterization would point to how broad or narrow the family of mechanisms that are possible is: i.e. how significant is it that the given sampling of models in the papers points to certain configurations. The closest I can see is the discussion point here that I think fails to completely address this point:
"One possibility suggested by our model is that the scalefree dynamics we observe occur at the boundary between winnerless switching and singlewinner lockedin dynamics (the red dashed line in Figures8E and F). Additional theoretical efforts are necessary to more fully explore how the traditional criticality hypothesis relates to the competitive criticality suggested by our model."
On the author's distinction between mechanisms and statistics, I do not believe the two are independent paths to choose between. Mechanisms have statistical signatures (socalled "statistical model") and data statistics inform which mechanisms are possible given data. At the end of the day, these are mathematical models that need to connect core concepts to data. My point for this particular paper, which I will try to say more clearly and succinctly now is that there are many possible mechanisms that could explain observed effect E. I was hoping in my prior review to spark a slightly longer discussion on what overall properties would such a family of mechanisms share. I fail to see how identifying core statistics that would pare down the possible mechanisms is at odds with looking for a mechanism that explains an effect.
https://doi.org/10.7554/eLife.79950.sa1Author response
Essential revisions:
1. The reviewers would like to see a detailed and thoughtful reflection on the role that 3 Hz Ca imaging might play in the conclusions that the authors derive. While the dataset in question offers many neurons, this approach is, from other perspectives, impoverished – calcium intrinsically misses spikes, a 3 Hz sampling rate is two orders of magnitude slower than an action potential, and the recordings are relatively short for amassing substantial observations of low probability (large) avalanches. The potential concern is that some of this disconnect may reflect optophysiological constraints. One argument against this is that a truly scale free system should be observable at any temporal or spatial scale and still give rise to the same sets of power laws. This quickly falls apart when applied to biological systems which are neither infinite in time nor space. As a result, the severe mismatch between the spatial resolution (single cell) and the temporal resolution (3 Hz) of the dataset, combined with filtering intrinsic to calcium imaging, raises the possibility that the conclusions are influenced by the methods.
We quite agree with the reviewer that reconciling different scales of measurement is an important and interesting question. One clue comes from Stringer et al’s original paper (2019 Science). They analyzed timeresolved spike data (from Neuropixel recordings) alongside the Ca imaging data we analyzed here. They showed that if the ephys spike data was analyzed with coarse time resolution (300 ms time bins, analogous to the Ca imaging data), then the anticorrelated activity became apparent (50/50 positive/negative loadings of PC1). When analyzed at faster time scales, anticorrelations were not apparent (mostly positive loadings of PC1). This interesting point was shown in their Supplementary Fig 12.
This finding suggests that our findings about anticorrelated neural groups may be relevant only at coarse time scales. Moreover, this point suggests that avalanche statistics may differ when analyzed at very different time scales, because the cancelation of anticorrelated groups may not be an important factor at faster timescales.
In our revised manuscript, we explored this point further by analyzing spike data from Stringer et al 2019. We focused on the spikes recorded from one local population (one Neuropixel probe). We first took the spike times of ~300 neurons and convolved them with a fast rise/slow fall, like typical Ca transient. Then we downsampled to 3 Hz sample rate. Next, we deconvolved using the same methods as those used by Stringer et al (OASIS nonnegative deconvolution). And finally, we zscored the resulting activity, as we did with the Ca imaging data. With this Calike signal in hand, we analyzed avalanches in four ways and compared the results. The four ways were: (1) the original timeresolved spikes (5 ms resolution), (2) the original spikes binned at 330 ms time res, (3) the full population of slow Calike signal, and (4) a correlated subset of neurons from the slow Calike signal. Based on the results of this new analysis (now in Figs S3 and S4), we found several interesting points that help reconcile potential differences between fast ephys and slow Ca signals:
1. In agreement with Sup Fig 12 from Stringer et al, anticorrelations are minimal in the fast, timeresolved spike data, but can be dominant in the slow, Calike signal.
2. Avalanche size distributions of spikes at fast timescales can exhibit a nice power law, consistent with previous results with exponents near 2 (e.g. Ma et al Neuron 2019, Fontenele et al PRL 2019). But, the same data at slow time scales exhibited poor powerlaws when the entire population was considered together.
3. The slow time scale data could exhibit a better power law if subsets of neurons were considered, just like our main findings based on Ca imaging. This point was the same using coarse timebinned spike data and the slow Calike signals, which gives us some confidence that deconvolution does not miss too many spikes.
In our opinion, a more thorough understanding of how scalefree dynamics differs across timescales will require a whole other paper, but we think these new results in our Figs S3 and S4 provide some reassurance that our results can be reconciled with previous work on scale free neural activity at faster timescales.
2. Another reservation expressed by the referees has to do with the generality of the conclusions drawn from the mechanistic model. One of the connectivity motifs identified appears to be i+ to e and i to e+, where potentially i+/i are SOM and VIP (or really any specific inhibitory type) cells. The specific connections to subsets of excitatory cells appear to be important (based on the solid lines in Figure 8). This seems surprising: is there any experimental support for excitatory cells to preferentially receive inhibition from either SOM or VIP, but not both?
There is indeed direct experimental support for the competitive relationship between SOM, VIP, and functionally distinct groups of excitatory neurons. This was shown in the paper by Josh Trachtenberg’s group: GarciaJuncoClemente et al 2017. An inhibitory pullpush circuit in frontal cortex. Nat Neurosci 20:389–392. However, we emphasize that we also showed (lower left motif in Fig 8G) that a simpler model with only one inhibitory group is sufficient to explain the anticorrelations and scalefree dynamics we observe. We opted to highlight the model with two inhibitory groups since it can also account for the GarciaJuncoClemente et al results.
In the section where we describe the model, we state, “We considered two inhibitory groups, instead of just one, to account for previous reports of anticorrelations between VIP and SOM inhibitory neurons in addition to anticorrelations between groups of excitatory neurons (GarciaJuncoClemente et al., 2017).”
More broadly, there was concern that the neat diagrams drawn here are misleading. The sample raster, showing what appears to be the full simulation, certainly captures the correlated/anticorrelated pattern of the 100 cells most correlated with a seed cell and 100 cells most anticorrelated with it, but it does not contain the 11,000 cells in between with zero to moderate levels of correlation.
We agree that our original model has several limitations and that one of the most obvious features lacking in our model is asynchronous neurons (The limitations are now discussed more openly in the last paragraph of the model subsection). In the data from the GarciaJuncoClemente et al paper above there are many asynchronous neurons as well. To ameliorate this limitation, we have now created a modified model that now accounts for asynchronous neurons together with the competing anticorrelated neurons (now shown and described in Fig 8 – figure supplement 2). We put this modified model in supplementary material and kept the simpler, original model in the main findings of our work, because the original model provides a simpler account of the features of the data we focused on in our work – i.e. anticorrelated scalefree fluctuations. The addition of the asynchronous population does not substantially change the behavior of the two anticorrelated groups in the original model.
We probably expect that the full covariance matrix has similar structure from any seed (see Meshulam et al. 2019, PRL, for an analysis of scaling of coarsegrained activity covariance), and this suggests multiple crossover inhibition constraints, which seem like they could be hard to satisfy.
We agree that it remains an outstanding challenge to create a model that reproduces the full complexity of the covariance matrix. We feel that this challenge is beyond the scope of this paper, which is already arguably squeezing quite a lot into one manuscript (one reviewer already suggested removing figures!).
We added a paragraph at the end of the subsection about the model to emphasize this limitation of the model as well as other limitations. This new paragraph says:
While our model offers a simple explanation of anticorrelated scalefree dynamics, its simplicity comes with limitations. Perhaps the most obvious limitation of our model is that it does not include neurons with weak correlations to both e+ and e (those neurons in the middle of the correlation spectra shown in Fig 7B). In Fig 8  figure supplement 2, we show that our model can be modified in a simple way to include asynchronous neurons. Another limitation is that we assumed that all nonzero synaptic connections were equal in weight. We loosen this assumption allowing for variable weights in Fig 8  figure supplement 2, without changing the basic features of anticorrelated scalefree fluctuations. Future work might improve our model further by accounting for neurons with intermediate correlations.
The motifs identified in Figure 8 likely exist, but I am left with many questions of what we learned about connectivity rules that would account for the full distribution of correlations. Would starting with an ErdosRenyi network with slight overrepresentation of these motifs be sufficient? How important is the homogeneous connection weights from each pool assumption – would allowing connection weights with some dispersion change the results?
First, we emphasize that our specific goal with our model was to identify a possible mechanism for the anticorrelated scalefree fluctuations that played the key role in our analyses. We agree that this is not a complete account of all correlations, but this was not the goal of our work. Nonetheless, our new modified model in Fig 8  figure supplement 2 now accounts for additional neurons with weak correlations. However, we think that future theoretical/modeling work will be required to better account for the intermediate correlations that are also present in the experimental data.
We confirmed that an ErdoRenyi network of E and I neurons can produce scalefree dynamics, but cannot produce substantial anticorrelated dynamics (Fig 8G, top right motif). Additionally, the parameter space study we performed with our model in Fig 8 showed that if the interactions between the two excitatory groups exceed a certain tipping point density, then the model behavior switches to behavior expected from an ErdosRenyi network (Fig 8F). Finally, we have now confirmed that some nonuniformity of synaptic weights does not change the main results (Fig 8  figure supplement 2). In the model presented in Fig 8  figure supplement 2, the value of each nonzero connection weight was drawn from a uniform distribution [0,0.01] or [0.01,0] for excitatory and inhibitory connections, respectively. All of these facts are described in the model subsection of the paper results.
3. Putting 2) another way, it's unclear why the averaging is required in the first place. This operation projects the entire population down in an incredibly lossy way and removes much of the complexity of the population activity.
Our population averaging approach is motivated by theoretical predictions and previous work. According to established theoretical accounts of scalefree population events (i.e. nonequilibrium critical phenomena in neural systems) such populationsummed event sizes should have power law statistics if the system is near a critical point. This approach has been used in many previous studies of scalefree neural activity (e.g. all of those cited in the introduction in relation to scalefree neuronal avalanches). One of the main results of our study is that the existing theories and models of critical dynamics in neural systems fail to account for small subsets of neurons with scalefree activity amid a larger population that does not conform to these statistics. We could not make this conclusion if we did not test the predictions of those existing theories and models.
Second, the authors state that it is highly curious that subsets of the population exhibit power laws while the entire population does not. While the discussion and hypothesizing about different ei interactions is interesting, it is possible that there's a discussion to be had on a much more basic level of whether there are topology independent explanations, such as basic distributions of correlations between neurons that can explain the subnetwork averaging. Specifically, if the correlation to any given neuron falls off, e.g., with an exponential falloff (i.e., a Gaussian Process type covariance between neurons), it seems that similar effects should hold. This type of effect can be easily tested by generating null distributions existing code bases. This is an important point since local (broadly defined) correlations of neurons implying the observed subnetwork behavior means that many mechanisms that have local correlations but don't cluster in any meaningful way could also be responsible for the local averaging effect.
We appreciate the reviewer’s effort, trying out some code to generate a statistical model. We agree that we could create such a *statistical model* that describes the observed distribution of pairwise correlations among neurons. For instance, it would be trivial to directly measure the covariance matrix, mean activities, and autocorrelations of the experimental data, which would, of course, provide a very good statistical description of the data. It would also be simple to generate more approximate statistical descriptions of the data, using multivariate gaussians, similar to the code suggested by the reviewer. However, we emphasize, this would not meet the goal of our modeling effort, which is mechanistic, not statistical. The aim of our model was to identify a possible biophysical mechanism from which emerge certain observed statistical features of the data. We feel that a statistical model is not a suitable strategy to meet this aim. Nonetheless, we agree with the reviewer that clusters with sharp boundaries (like the distinction between e+ an e in our model) are not necessary to reproduce the cancelation of anticorrelated neurons. In other words, we agree that sharp boundaries of the e+ and e groups of our model are not crucial ingredients to match our observations.
4. In general, the discussion of "two networks" seems like it relies on the correlation plot of Figure~7B. The decay away from the peak correlation is sharp, but there does not seem to be significant clustering in the anticorrelation population, instead a very slow decay away from zero. The authors do not show evidence of clustering in the neurons, nor any biophysical reason why e and i neurons are present in the imaging data.
First a small reminder: As stated in the paper, the data here is only showing activity of excitatory neurons. Inhibitory neurons are certainly present in V1, but they are not recorded in this data set. Thus we interpret our e+ and e groups as two subsets of anticorrelated excitatory neurons, like those we observed in the experimental data. We agree that our simplified model treats the anticorrelated subsets as if they are clustered, but this clustering is certainly not required for any of the data analyses of experimental data. We expect that our model could be improved to allow for a less sharp boundary between e+ and e groups, but we leave that for future work, because it is not essential to most of the results in the paper. This limitation of the model is now stated clearly in the last paragraph of the model subsection.
The alternative explanation (as mentioned in (b)) is that the there is a more continuous set of correlations among the neurons with the same result. In fact one of the reviewers tested this themself using code to generate some data with the desired statistics, and the distribution of events seems to also describe this same observation. Obviously, the full test would need to use the same event identification code, and so it is quite important that the authors consider the much more generic explanation for the subnetwork averaging effect. We recommend assessing the possibility that broader explanations, e.g., in the form of the distributions of correlations accounts for the observed phenomenon. Even with 10K neurons, there are many other forces at play influencing the observed network and while it is nice that ei networks are one explanation, much less constraining explanations that are still biophysically feasible should be discussed and compared against. I have provided one possible approach (see PDF for code and example figure: https://submit.elifesciences.org/eLife_files/2022/05/19/00107349/00/107349_0_attach_6_474_convrt.pdf).
As discussed above, we respectfully disagree that a statistical model is an acceptable replacement for a mechanistic model, since we are seeking to understand possible biophysical mechanisms. A statistical model is agnostic about mechanisms. We have nothing against statistical models, but in this case, they would not serve our goals.
To emphasize our point about the inadequacy of a statistical model for our goals, consider the following argument. Imagine we directly computed the mean activities, covariance matrix, and autocorrelations of all 10000 neurons from the real data. Then, we would have in hand an excellent statistical model of the data. We could then create a surrogate data set by drawing random numbers from a multivariate gaussian with same statistical description (e.g. using code like that offered by reviewer 3). This would, by construction, result in the same numbers of correlated and anticorrelated surrogate neurons. But what would this tell us about the biophysical mechanisms that might underlie these observations? Nothing, in our opinion.
5. Another important aspect here is how single neurons behave. It was not clear if single neurons were stated to exhibit a power law. If they do, then that would help in that there are different limiting behaviors to the averaging that pass through the observed stated numbers. If not, then there is an additional oddity that one must average neurons at all to obtain a power law. We recommend the authors show the full curve so that the readers get a more detailed sense of how averaging effects the powerlaw interpretation of the data.
We understand that our approach may seem odd from the point of view of centrallimittheoremtype argument. However, as mentioned above (reply R3b) and in our paper, there is a wellestablished history of theory and corresponding experimental tests for powerlaw distributed population events in neural systems near criticality. The prediction from theory is that the population summed activity will have powerlaw distributed events or fluctuations. That is the prediction that motivates our approach. In these theories, it is certainly not necessary that individual neurons have powerlaw fluctuations on their own. In most previous theories, it is necessary to consider the collective activity of many neurons before the powerlaw statistics become apparent, because each individual neurons contributes only a small part to the emergent, collective fluctuations. This phenomenon does not require that each individual neuron have powerlaw fluctuations.
At the risk of being pedantic, we feel obliged to point out that one cannot understand the peculiar scalefree statistics that occur at criticality by considering the behavior of individual elements of the system; hence the notion that critical phenomena are “emergent”. This important fact is not trivial and is, for example, why there was a Nobel prize awarded in physics for developing theoretical understanding of critical phenomena.
6. There is something that seems off about the range of \β values inferred with the ranges of \tau and $\α$. With \tau in [0.9,1.1], then the denominator 1\tau is in [0.1, 0.1], which the authors state means that \β (found to be in [2,2.4]) is not near \β_{crackling} = (\α1)/(1\tau). It seems as this is the opposite, as the possible values of the \β_{crackling} is huge due to the denominator, and so \β is in the range of possible \β_{crackling} almost vacuously. Was this statement just poorly worded?
The point here is that theory of crackling noise predicts that the fit value of beta should be equal to (1alpha)/(1tau). In other words, a confirmation of the theory would have all the points on the unity line in the rightmost panels of Fig9D and 9E, not scattered by more than an order of magnitude around the unity line. (We now state this explicitly in the text where Fig 9 is discussed.) Broad scatter around the unity line means the theory prediction did not hold. This is well established in previous studies of scalefree brain dynamics and crackling noise theory (see for example Ma et al Neuron 2019, Shew et al Nature Physics 2015, Friedman et al PRL 2012). A clearer single example of the failure of the theory to predict beta is shown in Fig 5A,B, and C.
7) It is not clear if there is more to what the authors are trying to say with the specifics of the scale free fits for behavior. Apparently, these results are used to motivate the neural studies, but aside from that, the details of those ranges don't seem to come up again.
The reviewer is correct, the primary point in Fig 2 is that scalefree behavioral statistics often exist. Beyond this point about existence, reporting of the specific exponents and ranges is just standard practice for this kind of analysis; a natural question to ask after claiming that we find scale behavior is “what are the exponents and ranges”. We would be remiss not to report those numbers.
Given that the primary connection between neuronal and behavioral activity seems to be Figure 4. The distribution of points in these plots seem to be very lopsided, in that some plots have large ranges of fewtono data points.
We agree that this whitespace in the figure panels is a somewhat awkward, but we chose to keep the horizontal axis the same for all panels of Fig 4B, because this shows that not all behaviors, and not all animals had the same range of behavioral correlations. We felt that hiding this was a bit misleading, so we kept the white space.
It would be very helpful to get a sense of the distribution of points that are a bit hard to see given the overlapping points and superimposed lines. We recommend that the authors add distribution information to the plots in Figure 4B to give a sense of how points are spread through the [correlation with behavior]by[power law range] space. Potential plots might be a colocated histogram, or perhaps an uncertainty estimate as a function of correlation based on the number of points and variance. This would help show significance of the curves in a way that accounts for the uneven spread of datapoints.
We also agree that characterizing uncertainty in power law range as a function of correlation and providing distribution information are good ideas. If we have not misunderstood the reviewer’s suggestion, we think we have already done such characterization. The three lines in each panel of Fig 4B meet the goal of characterizing variability as a function of correlation. The middle line is median, the top and bottom lines span the quartile range, which is a good way to characterize variability around the median for nonnormally distributed variability. We also provide information on how points are distributed by plotting the points with partially transparent markers. In this way, higher density of overlapping points creates darker regions in the clouds of points. We feel that this approach avoids hiding outlier points and accounts for distributions of overlapping points. Adding more information (actual distributions, e.g.) to these plots would be largely redundant and make them too cluttered. We added the following sentence to Fig 4 caption: “Points have partially transparent markers, thus darker areas reveal higher density of points.”
8. Neural activity correlated with some behavior variables can sometimes be the most active subset of neurons. This could potentially skew the maximum sizes of events and give behaviorally correlated subsets an unfair advantage in terms of the scalefree range. In a similar vain to 8), what are the typical dynamic ranges for subsets correlated and uncorrelated with behavior? We recommend showing a number of these to see if those dynamic ranges are impacting the possible ranges in the [correlation with behavior]by[power law range] plots. Perhaps something like curve in each plot showing the minimum maximum value of the power law range per correlation range.
In general, the reviewers struggled with the interpretation of Figure 4b in the sense that there seems to be such variability between mice. How much do the authors feel that this is a difference in neural populations imaged, vs changes in imaging conditions (illumination, window clarity, optical alignment) or differences in mouse activity levels?
Here we follow a standard convention for Ca imaging data (e.g. original Stringer 2019 Science paper). Each neuron’s activity was zscored before defining groups and events. This means that the “dynamic ranges” (ranges of activity amplitudes) of each single neuron were similar. This should “normalize” the maximum possible event size in a fair way. We note, however, that the primary reason this zscoring is done is that the measurement signaltonoise ratio can vary across neurons due to variable expression of fluorescent indicator molecules across neurons.
As mentioned in above, the scatter of points in Fig 4 panels already directly shows the mintomax span of power law range for each correlation.
Reviewer #1 (Recommendations for the authors):
This paper feels highly polished and thorough in its presentation. I truly enjoyed reading it and believe it will be of value to the community.
We truly appreciate the positive feedback. Thanks
A nuanced question: mouse #5 and mouse #6 consistently seems to break the rules implied by the rest of the dataset. The behavioral descriptors look normal in figure 2, but neural structure is notably different from the rest of the group in 3E, 4B, 6C, 7E, and S4 (that one is behavior). Is there anything meaningfully different about these recordings (n cells, mean event rate, anatomical location etc)?
We agree that 3 of the 9 recordings (mice 57) are somewhat “outlying” relative to the results from the other 6 recordings. We took a closer look at those recordings and found that at least one possible reason is that the animals were much more active than usual (mice 6 and 7 were 50 times more active the next most active mouse, in terms of median run speed) or much less activity than usual (mouse 5 was 10 time less active than the next least active mouse). We have pointed this out in the main text near the description of Figure 4 and showed it directly in Figure 4 —figure supplement 1.
Reviewer #2 (Recommendations for the authors):Comments, questions, and technical issues for the authors.
1. The values of tau (size exponent) from neural avalanches are surprisingly low. For instance, in the Ma et al. (Hengen) 2019 Neuron paper, tau ranged from 1.5 to 1.9, sometimes more than 2, but never less than 1. I am wondering how the imaging preprocessing affects the estimate of tau. For instance, if you started with simulated spiking activity that has avalanches with size pdf that go as a power law with exponent tau, and then use a forward model to generate "calcium imaging," and then apply deconvolution, zscoring, and lowpass filtering, and then measure the avalanches again: what is the new tau?
We appreciate the suggestion to do some forward modeling of Calike signals from spikes. This is the approach we have now taken for FiguresS3 and S4.
2. This may be an ignorant question (apologies). The power law range is quantified in decibels (dB) throughout the paper; do you actually mean decades?
The typical definition of dB is the ratio of two log quantities and then multiplied by a factor of 10 or 20 (depending on the context). In our usage, we did the ratio of two log quantities, but we never multiplied by 10 or 20. To avoid confusion, we have replaced “dB” with “decades” throughout the paper.
3. Related to the setup of the model, was there a reason that there are no adaptation mechanisms in this network model, as there often are in mechanistic models for avalanche criticality (including past work by the authors of this paper, e.g. Shew's 2015 Nature Physics paper)? Also, there appears to have been an error with the reference manager, as this reference shows up twice in the reference list.
There was no reason for the omission of adaptation mechanisms other than to increase simplicity. We aimed to make the model as simple as possible and still explain the features of the data we were most concerned with (i.e. scalefreeness and anticorrelations).
Thanks for catching the double reference. We fixed that now.
4. It would be helpful the authors could elaborate on predictions that their results make for future studies. Maybe this is rather technical, but can you tell us when you expect to find a powerlaw distribution as a function of how much of the population is sampled and for how long? What if you were analyzing Neuropixels data, where you lose the extensive spatial sampling (and the restriction to pyramidal cells only) but you gain 3 orders of magnitude in temporal resolution?
Please see above and Figure 3 —figure supplements 12.
5. On page 8, you ask "are all behavioral events equally correlated to their concurrent neural events, or are certain neural events from certain subsets of neurons more strongly related to behavioral events?" I don't understand the question. What does it mean for all behavioral events to be equally correlated to concurrent neural events? Aren't "concurrent neural events" in specific subsets of neurons?
We agree that the wording of this sentence was confusing, so we deleted it. We think the goals of the analysis in Figure 6 are clearer with this deletion.
Reviewer #3 (Recommendations for the authors):
1. Limits of calcium imaging: My recommendation is to assess mathematically the potential impact of missing data on the range and powerlaw slope estimates, which are the primary values used throughout the paper.
We appreciate the concern here and agree that event size statistics could in principle be biased in some systematic way due to missed spikes due to deconvolution of Ca signals. To directly test this possibility, we performed a new analysis of spike data recorded with high time resolution electrophysiology. We began with forwardmodeling process to create a lowtimeresolution, Calike signal, using the same deconvolution algorithm (OASIS) that was used to generate the data we analyzed in our work here. In agreement with the reviewer’s concern, we found that spikes were sometimes missed, but the loss was not extreme and did not impact the neural event size statistics in a significant way compared to the ground truth we obtained directly from the original spike data (with no loss of spikes). This new work is now described in a new paragraph at the end of the subsection of results related to Fig 3 and in a new Fig 3  figure supplement 1. The new paragraph says…
“Two concerns with the data analyzed here are that it was sampled at a slow time scale (3 Hz frame rate) and that the deconvolution methods used to obtain the data here from the raw GCAMP6s Ca imaging signals are likely to miss some activity (Huang et al., 2021). Since our analysis of neural events hinges on summing up activity across neurons, could it be that the missed activity creates systematic biases in our observed event size statistics? To address this question, we analyzed some timeresolved spike data (Neuropixel recording from Stringer et al 2019). Starting from the spike data, we created a slow signal, similar to that we analyzed here by convolving with a Catransient, down sampling, deconvolving, and zscoring (Fig 3  figure supplement 1). We compared neural event size distributions to “ground truth” based on the original spike data (with no loss of spikes) and found that the neural event size distributions were very similar, with the same exponent and same powerlaw range (Fig 3  figure supplement 1). Thus, we conclude that our reported neural event size distributions are reliable.”
However, although loss of spikes did not impact the event size distributions much, the timescale of measurement did matter. As discussed above and shown in Fig 3  figure supplement 2, changing from 5 ms time resolution to 330 ms time resolution does change the exponent and the range of the power law. However, in the test data set we worked with, the existence of a power law was robust across time scales.
2. Correlations and powerlaws in subsets.
2ac. My recommendation is to assess the possibility that broader explanations, e.g., in the form of the distributions of correlations accounts for the observed phenomenon. Even with 10K neurons, there are many other forces at play influencing the observed network and while it is nice that ei networks are one explanation, much less constraining explanations that are still biophysically feasible should be discussed and compared against. I have provided one possible approach (see PDF for code and example figure: https://submit.elifesciences.org/eLife_files/2022/05/19/00107349/00/107349_0_attach_6_474_convrt.pdf) that I hope will be useful to the authors.
discussed above, we respectfully disagree that a statistical model is an acceptable replacement for a mechanistic model, since we are seeking to understand possible biophysical mechanisms. A statistical model is agnostic about mechanisms. We have nothing against statistical models, but in this case, they would not serve our goals.
To emphasize our point about the inadequacy of a statistical model for our goals, consider the following argument. Imagine we directly computed the mean activities, covariance matrix, and autocorrelations of all 10000 neurons from the real data. Then, we would have in hand an excellent statistical model of the data. We could then create a surrogate data set by drawing random numbers from a multivariate gaussian with same statistical description (e.g. using code like that offered by reviewer 3). This would, by construction, result in the same numbers of correlated and anticorrelated surrogate neurons. But what would this tell us about the biophysical mechanisms that might underlie these observations? Nothing, in our opinion.
2d. I recommend the authors show the full curve so that the readers get a more detailed sense of how averaging effects the powerlaw interpretation of the data.
It is not clear to us what is meant by the “full curve”.
3. Please check that this calculation and interpretation is correct.
The point here is that theory of crackling noise predicts that the fit value of beta should be equal to (1alpha)/(1tau). In other words, a confirmation of the theory would have all the points on the unity line in the rightmost panels of Fig9D and 9E, not scattered by more than an order of magnitude around the unity line. (We now state this explicitly in the text where Fig 9 is discussed.) Broad scatter around the unity line means the theory prediction did not hold. This is well established in previous studies of scalefree brain dynamics and crackling noise theory (see for example Ma et al Neuron 2019, Shew et al Nature Physics 2015, Friedman et al PRL 2012). A clearer single example of the failure of the theory to predict beta is shown in Fig 5A,B, and C.
4. Connection between brain and behavior:
4b. I recommend that the authors add distribution information to the plots in Figure~4B to give a sense of how points are spread through the [correlation with behavior]by[power law range] space. Potential plots might be a colocated histogram, or perhaps an uncertainty estimate as a function of correlation based on the number of points and variance. This would help show significance of the curves in a way that accounts for the uneven spread of datapoints.
We also agree that characterizing uncertainty in power law range as a function of correlation and providing distribution information are good ideas. If we have not misunderstood the reviewer’s suggestion, we think we have already done such characterization. The three lines in each panel of Figure 4B meet the goal of characterizing variability as a function of correlation. The middle line is median, the top and bottom lines span the quartile range, which is a good way to characterize variability around the median for nonnormally distributed variability. We also provide information on how points are distributed by plotting the points with partially transparent markers. In this way, higher density of overlapping points creates darker regions in the clouds of points. We feel that this approach avoids hiding outlier points and accounts for distributions of overlapping points. Adding more information (actual distributions, e.g.) to these plots would be largely redundant and make them too cluttered. We added the following sentence to Figure 4 caption:
“Points have partially transparent markers, thus darker areas reveal higher density of points.”
4c. In a similar vein, what are the typical dynamic ranges for subsets correlated and uncorrelated with behavior? I recommend showing a number of these to see if those dynamic ranges are impacting the possible ranges in the [correlation with behavior]by[power law range] plots. Perhaps something like curve in each plot showing the minimum maximum value of the power law range per correlation range.
Here we follow a standard convention for Ca imaging data (e.g. original Stringer 2019 Science paper). Each neuron’s activity was zscored before defining groups and events. This means that the “dynamic ranges” (ranges of activity amplitudes) of each single neuron were similar. This should “normalize” the maximum possible event size in a fair way. We note, however, that the primary reason this zscoring is done is that the measurement signaltonoise ratio can vary across neurons due to variable expression of fluorescent indicator molecules across neurons.
As mentioned in R3j, the scatter of points in Figure 4 panels already directly shows the mintomax span of power law range for each correlation.
4d. In general I'm struggling with the interpretation of Figure~4b in the sense that there seems to be such variability between mice. How much do the authors feel that this is a difference in neural populations imaged, vs changes in imaging conditions (illumination, window clarity, optical alignment) or differences in mouse activity levels?
We agree that 3 of the 9 recordings (mice 57) are somewhat “outlying” relative to the results from the other 6 recordings. We took a closer look at those recordings and found that at least one possible reason is that the animals were much more active than usual (mice 6 and 7 were 50 times more active the next most active mouse, in terms of median run speed) or much less activity than usual (mouse 5 was 10 time less active than the next least active mouse). We have pointed this out in the main text near the description of Fig 4 and showed it directly in Fig 4 – figure supplement 1.
[Editors’ note: further revisions were suggested prior to acceptance, as described below.]
Reviewer #3 (Recommendations for the authors):
I appreciate your time and effort to respond to my review points. I believe the majority of my points have been addressed. The main weakness I still see is the lack of discussion about the broader mechanisms beyond the basic structures described that could account for the observations (see below).
Mechanisms vs. Statistics
I would like to clarify my reason for bringing up this "statistical" viewpoint that I believe may have been lost in translation. The paper as it stands makes the following logical steps (in terms of the mechanistic model) (1) The data exhibits power law scaling under certain binning of neurons (effect E) and (2) One way to account for this effect E is to consider certain ei models. This is reasonable, but the overall search for possible mechanisms seems to be a combination of intuition and trial and error. Different architectures were tried and either "passed" or "failed".
The purpose of bringing up the statistical properties needed was to hopefully raise a conversation around what core properties are needed to replicate this effect. As activity statistics and connectivity are intimately related in mechanistic models, such a characterization would point to how broad or narrow the family of mechanisms that are possible is: i.e. how significant is it that the given sampling of models in the papers points to certain configurations. The closest I can see is the discussion point here that I think fails to completely address this point:
"One possibility suggested by our model is that the scalefree dynamics we observe occur at the boundary between winnerless switching and singlewinner lockedin dynamics (the red dashed line in Figures8E and F). Additional theoretical efforts are necessary to more fully explore how the traditional criticality hypothesis relates to the competitive criticality suggested by our model."
On the author's distinction between mechanisms and statistics, I do not believe the two are independent paths to choose between. Mechanisms have statistical signatures (socalled "statistical model") and data statistics inform which mechanisms are possible given data. At the end of the day, these are mathematical models that need to connect core concepts to data. My point for this particular paper, which I will try to say more clearly and succinctly now is that there are many possible mechanisms that could explain observed effect E. I was hoping in my prior review to spark a slightly longer discussion on what overall properties would such a family of mechanisms share. I fail to see how identifying core statistics that would pare down the possible mechanisms is at odds with looking for a mechanism that explains an effect.
We appreciate the clarifications and agree that mechanisms and statistics are intimately related. Indeed, in this context, much of our main results (event size and duration distributions, correlation spectra, etc.) can be viewed as statistical properties. Admittedly, these are rather “high level” statistics. To expand on this, we have now quantified some more basic statistics. We have now added a supplementary figure with single neuron variance and pairwise covariance values, comparing between neural subsets that exhibited good power laws and those that did not (Figure 8 —figure supplement 3). We show in this new figure that neural subsets with large powerlaw range tended to be composed of neurons with relatively high variance and high pairwise covariance. Perhaps these statistical properties will be shared by the family of mechanistic models suggested by the reviewer. However, in our view, it is beyond the scope of our manuscript to fully identify the necessary and sufficient statistical properties a model must generate to agree with our main results.
https://doi.org/10.7554/eLife.79950.sa2Article and author information
Author details
Funding
National Institutes of Health (NIHR15NS116742)
 Woodrow L Shew
National Science Foundation (NSF1912352)
 Woodrow L Shew
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
Ethics
All experimental procedures were conducted according to the UK Animals Scientific Procedures Act (1986). Experiments were performed at University College London under personal and project licenses released by the Home Office following appropriate ethics review.
Senior Editor
 Timothy E Behrens, University of Oxford, United Kingdom
Reviewing Editor
 Gordon J Berman, Emory University, United States
Reviewer
 Keith B Hengen, Washington University, United States
Version history
 Preprint posted: May 13, 2021 (view preprint)
 Received: May 3, 2022
 Accepted: January 6, 2023
 Accepted Manuscript published: January 27, 2023 (version 1)
 Version of Record published: February 15, 2023 (version 2)
Copyright
© 2023, Jones et al.
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.
Metrics

 1,831
 Page views

 314
 Downloads

 6
 Citations
Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.
Download links
Downloads (link to download the article as PDF)
Open citations (links to open the citations from this article in various online reference manager services)
Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)
Further reading

 Neuroscience
Blindness affects millions of people around the world. A promising solution to restoring a form of vision for some individuals are cortical visual prostheses, which bypass part of the impaired visual pathway by converting camera input to electrical stimulation of the visual system. The artificially induced visual percept (a pattern of localized light flashes, or ‘phosphenes’) has limited resolution, and a great portion of the field’s research is devoted to optimizing the efficacy, efficiency, and practical usefulness of the encoding of visual information. A commonly exploited method is noninvasive functional evaluation in sighted subjects or with computational models by using simulated prosthetic vision (SPV) pipelines. An important challenge in this approach is to balance enhanced perceptual realism, biologically plausibility, and realtime performance in the simulation of cortical prosthetic vision. We present a biologically plausible, PyTorchbased phosphene simulator that can run in realtime and uses differentiable operations to allow for gradientbased computational optimization of phosphene encoding models. The simulator integrates a wide range of clinical results with neurophysiological evidence in humans and nonhuman primates. The pipeline includes a model of the retinotopic organization and cortical magnification of the visual cortex. Moreover, the quantitative effects of stimulation parameters and temporal dynamics on phosphene characteristics are incorporated. Our results demonstrate the simulator’s suitability for both computational applications such as endtoend deep learningbased prosthetic vision optimization as well as behavioral experiments. The modular and opensource software provides a flexible simulation framework for computational, clinical, and behavioral neuroscientists working on visual neuroprosthetics.

 Neuroscience
Extinction is a specific example of learning where a previously reinforced stimulus or response is no longer reinforced, and the previously learned behaviour is no longer necessary and must be modified. Current theories suggest extinction is not the erasure of the original learning but involves new learning that acts to suppress the original behaviour. Evidence for this can be found when the original behaviour recovers following the passage of time (spontaneous recovery) or reintroduction of the reinforcement (i.e. reinstatement). Recent studies have shown that pharmacological manipulation of noradrenaline (NA) or its receptors can influence appetitive extinction; however, the role and source of endogenous NA in these effects are unknown. Here, we examined the role of the locus coeruleus (LC) in appetitive extinction. Specifically, we tested whether optogenetic stimulation of LC neurons during extinction of a foodseeking behaviour would enhance extinction evidenced by reduced spontaneous recovery in future tests. LC stimulation during extinction trials did not change the rate of extinction but did serve to reduce subsequent spontaneous recovery, suggesting that stimulation of the LC can augment rewardrelated extinction. Optogenetic inhibition of the LC during extinction trials reduced responding during the trials where it was applied, but no longlasting changes in the retention of extinction were observed. Since not all LC cells expressed halorhodopsin, it is possible that more complete LC inhibition or pathwayspecific targeting would be more effective at suppressing extinction learning. These results provide further insight into the neural basis of appetitive extinction, and in particular the role of the LC. A deeper understanding of the physiological bases of extinction can aid development of more effective extinctionbased therapies.