Krotov and Hopfield (2021) proposed a biologically plausible two-layer associative memory network with memory storage capacity exponential in the number of visible neurons. However, the capacity was only linear in the number of hidden neurons. This limitation arose from the choice of nonlinearity between the visible and hidden units, which enforced winner-takes-all dynamics in the hidden layer, thereby restricting each hidden unit to encode only a single memory. We overcome this limitation by introducing a novel associative memory network with a threshold nonlinearity that enables distributed representations. In contrast to winner-takes-all dynamics, where each hidden neuron is tied to an entire memory, our network allows hidden neurons to encode basic components shared across many memories. Consequently, complex patterns are represented through combinations of hidden neurons. These representations reduce redundancy and allow many correlated memories to be stored compositionally. Thus, we achieve much higher capacity: exponential in the number of hidden units, provided the number of visible units is sufficiently larger than the number of hidden neurons. Exponential capacity arises because all binary states of the hidden units can become stable memory patterns with an appropriately chosen threshold. Moreover, the distributed hidden representation, which has much lower dimensionality than the visible layer, preserves class-discriminative structure, supporting efficient nonlinear decoding. These results establish a new regime for associative memory, enabling high-capacity, robust, and scalable architectures consistent with biological constraints.

Why do neurons encode information the way they do? Normative answers to this question model neural activity as the solution to an optimisation problem; for example, the celebrated efficient coding hypothesis frames neural activity as the optimal encoding of information under efficiency constraints. Successful normative theories have varied dramatically in complexity, from simple linear models (Atick & Redlich, 1990), to complex deep neural networks (Lindsay, 2021). What complex models gain in flexibility, they lose in tractability and often understandability. Here, we split the difference by constructing a set of tractable but flexible normative representational theories. Instead of optimising the neural activities directly, following (Sengupta et al. 2018), we instead optimise the representational similarity, a matrix formed from the dot products of each pair of neural responses. Using this, we show that a large family of interesting optimisation problems are convex. This includes problems corresponding to linear and some non-linear neural networks, and problems from the literature not previously recognised as convex such as modified versions of semi-nonnegative matrix factorisation or nonnegative sparse coding. We put these findings to work in two ways. First, we extend previous results on modularity and mixed selectivity in neural activity; in so doing we provide the first necessary and sufficient identifiability result for a form of semi-nonnegative matrix factorisations. Second, we seek to understand the meaningfulness of single neural tuning curves as compared to neural representations. In particular we derive an identifiability result stating that, for an optimal representational similarity matrix, if neural tunings are `different enough' then they are uniquely linked to the optimal representational similarity, partially justifying the use of single neuron tuning analysis in neuroscience. In sum, we identify an interesting space of convex problems, and use that to derive neural coding results.

Neural networks trained with gradient descent often learn solutions of increasing complexity over time, a phenomenon known as simplicity bias. Despite being widely observed across architectures, existing theoretical treatments lack a unifying framework. We present a theoretical framework that explains a simplicity bias arising from saddle-to-saddle learning dynamics for a general class of neural networks, incorporating fully-connected, convolutional, and attention-based architectures. Here, simple means expressible with few hidden units, i.e., hidden neurons, convolutional kernels, or attention heads. Specifically, we show that linear networks learn solutions of increasing rank, ReLU networks learn solutions with an increasing number of kinks, convolutional networks learn solutions with an increasing number of convolutional kernels, and self-attention models learn solutions with an increasing number of attention heads. By analyzing fixed points, invariant manifolds, and dynamics of gradient descent learning, we show that saddle-to-saddle dynamics operates by iteratively evolving near an invariant manifold, approaching a saddle, and switching to another invariant manifold. Our analysis also illuminates the effects of data distribution and initialization on the duration and number of plateaus in learning, dissociating previously confounding factors. Overall, our theory offers a framework for understanding when and why gradient descent progressively learns increasingly complex solutions.

Why do biological and artificial neurons sometimes modularise, each encoding a single meaningful variable, and sometimes entangle their representation of many variables? In this work, we develop a theory of when biologically inspired networks---those that are nonnegative and energy efficient---modularise their representation of source variables (sources). We derive necessary and sufficient conditions on a sample of sources that determine whether the neurons in an optimal biologically-inspired linear autoencoder modularise. Our theory applies to any dataset, extending far beyond the case of statistical independence studied in previous work. Rather we show that sources modularise if their support is ``sufficiently spread''. From this theory, we extract and validate predictions in a variety of empirical studies on how data distribution affects modularisation in nonlinear feedforward and recurrent neural networks trained on supervised and unsupervised tasks. Furthermore, we apply these ideas to neuroscience data, showing that range independence can be used to understand the mixing or modularising of spatial and reward information in entorhinal recordings in seemingly conflicting experiments. Further, we use these results to suggest alternate origins of mixed-selectivity, beyond the predominant theory of flexible nonlinear classification. In sum, our theory prescribes precise conditions on when neural activities modularise, providing tools for inducing and elucidating modular representations in brains and machines.

AbstractDuring many tasks the brain receives real-time feedback about performance. What should it do with that information, at the synaptic level, so that tasks can be performed as well as possible? The conventional answer is that it should learn by incrementally adjusting synaptic strengths. We show, however, that learning on its own is severely suboptimal. To maximize performance, synaptic plasticity should also operate on a much faster timescale – essentially, the synaptic weights should act as a control signal. We propose a normative plasticity rule that embodies this principle. In this, fast synaptic weight changes greedily suppress downstream errors, while slow synaptic weight changes implement statistically optimal learning. This enables near-perfect task performance immediately, efficient task execution on longer timescales, and confers robustness to noise and other perturbations. Applied in a cerebellar microcircuit model, the theory explains longstanding experimental observations and makes novel testable predictions.

While attention-based models have demonstrated the remarkable ability of in-context learning (ICL), the theoretical understanding of how these models acquired this ability through gradient descent training is still preliminary. Towards answering this question, we study the gradient descent dynamics of multi-head linear self-attention trained for in-context linear regression. We examine two parametrizations of linear self-attention: one with the key and query weights merged as a single matrix (common in theoretical studies), and one with separate key and query matrices (closer to practical settings). For the merged parametrization, we show that the training dynamics has two fixed points and the loss trajectory exhibits a single, abrupt drop. We derive an analytical time-course solution for a certain class of datasets and initialization. For the separate parametrization, we show that the training dynamics has exponentially many fixed points and the loss exhibits saddle-to-saddle dynamics, which we reduce to scalar ordinary differential equations. During training, the model implements principal component regression in context with the number of principal components increasing over training time. Overall, we provide a theoretical description of how ICL abilities evolve during gradient descent training of linear attention, revealing abrupt acquisition or progressive improvements depending on how the key and query are parametrized.

A remarkable demonstration of the flexibility of mammalian motor systems is primates’ ability to learn to control brain-computer interfaces (BCIs). This constitutes a completely novel motor behavior, yet primates are capable of learning to control BCIs under a wide range of conditions. BCIs with carefully calibrated decoders, for example, can be learned with only minutes to hours of practice. With a few weeks of practice, even BCIs with randomly constructed decoders can be learned. What are the biological substrates of this learning process? Here, we develop a theory based on a re-aiming strategy, whereby learning operates within a low-dimensional subspace of task-relevant inputs driving the local population of recorded neurons. Through comprehensive numerical and formal analysis, we demonstrate that this theory can provide a unifying explanation for disparate phenomena previously reported in three different BCI learning tasks, and we derive a novel experimental prediction that we verify with previously published data. By explicitly modeling the underlying neural circuitry, the theory reveals an interpretation of these phenomena in terms of biological constraints on neural activity.

The neural representations of prior information about the state of the world are poorly understood1. Here, to investigate them, we examined brain-wide Neuropixels recordings and widefield calcium imaging collected by the International Brain Laboratory. Mice were trained to indicate the location of a visual grating stimulus, which appeared on the left or right with a prior probability alternating between 0.2 and 0.8 in blocks of variable length. We found that mice estimate this prior probability and thereby improve their decision accuracy. Furthermore, we report that this subjective prior is encoded in at least 20% to 30% of brain regions that, notably, span all levels of processing, from early sensory areas (the lateral geniculate nucleus and primary visual cortex) to motor regions (secondary and primary motor cortex and gigantocellular reticular nucleus) and high-level cortical regions (the dorsal anterior cingulate area and ventrolateral orbitofrontal cortex). This widespread representation of the prior is consistent with a neural model of Bayesian inference involving loops between areas, as opposed to a model in which the prior is incorporated only in decision-making areas. This study offers a brain-wide perspective on prior encoding at cellular resolution, underscoring the importance of using large-scale recordings on a single standardized task.

Using multiple input streams simultaneously to train multimodal neural networks is intuitively advantageous but practically challenging. A key challenge is unimodal bias, where a network overly relies on one modality and ignores others during joint training. We develop a theory of unimodal bias with multimodal deep linear networks to understand how architecture and data statistics influence this bias. This is the first work to calculate the duration of the unimodal phase in learning as a function of the depth at which modalities are fused within the network, dataset statistics, and initialization. We show that the deeper the layer at which fusion occurs, the longer the unimodal phase. A long unimodal phase can lead to a generalization deficit and permanent unimodal bias in the overparametrized regime. Our results, derived for multimodal linear networks, extend to nonlinear networks in certain settings. Taken together, this work illuminates pathologies of multimodal learning under joint training, showing that late and intermediate fusion architectures can give rise to long unimodal phases and permanent unimodal bias. Our code is available at: https://yedizhang.github.io/unimodal-bias.html.

We investigate the implications of removing bias in ReLU networks regarding their expressivity and learning dynamics. We first show that two-layer bias-free ReLU networks have limited expressivity: the only odd function two-layer bias-free ReLU networks can express is a linear one. We then show that, under symmetry conditions on the data, these networks have the same learning dynamics as linear networks. This enables us to give analytical time-course solutions to certain two-layer bias-free (leaky) ReLU networks outside the lazy learning regime. While deep bias-free ReLU networks are more expressive than their two-layer counterparts, they still share a number of similarities with deep linear networks. These similarities enable us to leverage insights from linear networks to understand certain ReLU networks. Overall, our results show that some properties previously established for bias-free ReLU networks arise due to equivalence to linear networks.

Biological synaptic transmission is unreliable, and this unreliability likely degrades neural circuit performance. While there are biophysical mechanisms that can increase reliability, for instance by increasing vesicle release probability, these mechanisms cost energy. We examined four such mechanisms along with the associated scaling of the energetic costs. We then embedded these energetic costs for reliability in artificial neural networks (ANN) with trainable stochastic synapses, and trained these networks on standard image classification tasks. The resulting networks revealed a tradeoff between circuit performance and the energetic cost of synaptic reliability. Additionally, the optimised networks exhibited two testable predictions consistent with pre-existing experimental data. Specifically, synapses with lower variability tended to have 1) higher input firing rates and 2) lower learning rates. Surprisingly, these predictions also arise when synapse statistics are inferred through Bayesian inference. Indeed, we were able to find a formal, theoretical link between the performance-reliability cost tradeoff and Bayesian inference. This connection suggests two incompatible possibilities: evolution may have chanced upon a scheme for implementing Bayesian inference by optimising energy efficiency, or alternatively, energy efficient synapses may display signatures of Bayesian inference without actually using Bayes to reason about uncertainty.

Using multiple input streams simultaneously in training multimodal neural networks is intuitively advantageous, but practically challenging. A key challenge is unimodal bias, where a network overly relies on one modality and ignores others during joint training. While unimodal bias is well-documented empirically, our theoretical understanding of how architecture and data statistics influence this bias remains incomplete. Here we develop a theory of unimodal bias with deep multimodal linear networks. We calculate the duration of the unimodal phase in learning, as a function of the depth at which modalities are fused within the network, dataset statistics, and initialization. We find that the deeper the layer at which fusion occurs, the longer the unimodal phase. In addition, our theory reveals the modality learned first is not necessarily the modality that contributes more to the output. Our results, derived for multimodal linear networks, extend to ReLU networks in certain settings. Taken together, this work illuminates pathologies of multimodal learning under joint training, showing that late and intermediate fusion architectures can give rise to long unimodal phases and even prioritize learning a less helpful modality.

Full text Figures and data Side by side Abstract eLife assessment Introduction Results Discussion Materials and methods Data availability References Peer review Author response Article and author information Abstract Observations of power laws in neural activity data have raised the intriguing notion that brains may operate in a critical state. One example of this critical state is 'avalanche criticality', which has been observed in various systems, including cultured neurons, zebrafish, rodent cortex, and human EEG. More recently, power laws were also observed in neural populations in the mouse under an activity coarse-graining procedure, and they were explained as a consequence of the neural activity being coupled to multiple latent dynamical variables. An intriguing possibility is that avalanche criticality emerges due to a similar mechanism. Here, we determine the conditions under which latent dynamical variables give rise to avalanche criticality. We find that populations coupled to multiple latent variables produce critical behavior across a broader parameter range than those coupled to a single, quasi-static latent variable, but in both cases, avalanche criticality is observed without fine-tuning of model parameters. We identify two regimes of avalanches, both critical but differing in the amount of information carried about the latent variable. Our results suggest that avalanche criticality arises in neural systems in which activity is effectively modeled as a population driven by a few dynamical variables and these variables can be inferred from the population activity. eLife assessment This paper provides a simple example of a neural-like system that displays criticality, but not for any deep reason; it's just because a population of neurons are driven (independently!) by a slowly varying latent variable, something that is common in the brain. Moreover, criticality does not imply optimal information transmission (one of its proposed functions). The work is likely to have an important impact on the study of criticality in neural systems and is convincingly supported by the experiments presented. https://doi.org/10.7554/eLife.89337.3.sa0 About eLife assessments Introduction The neural criticality hypothesis – the idea that neural systems operate close to a phase transition, perhaps for optimal information processing – is both ambitious and banal. Measurements from biological systems are limited in the range of spatial and temporal scales that can be sampled, not only because of the limitations of recording techniques but also due to the fundamentally non-stationary behavior of most, if not all, biological systems. These limitations make proving that an observation indicates critical behavior difficult. At the same time, the idea that brain networks are critical echoes the anthropic principle: tuned another way, a network becomes quiescent or epileptic and in either state, seems unlikely to support perception, thought, or flexible behavior, yet these observations do not explain how such fine-tuning could be achieved. Further muddying the water, researchers have reported multiple kinds of criticality in neural networks, including through analysis of avalanches (Beggs and Plenz, 2003; Plenz et al., 2021; O'Byrne and Jerbi, 2022; Girardi-Schappo, 2021) and of coarse-grained activity (Meshulam et al., 2019), as well as of correlations (Dahmen et al., 2019). How these flavors of critical behavior relate to each other or any functional network mechanism is unknown. The phenomenon that we will refer to as 'avalanche criticality' appears remarkably widespread. It was first observed in cultured neurons (Beggs and Plenz, 2003) and later studied in zebrafish (Ponce-Alvarez et al., 2018), turtles (Shew et al., 2015), rodents (Ma et al., 2019), monkeys (Petermann et al., 2009), and even humans (Poil et al., 2008). The standard analysis, described later, requires extracting power-law exponents from fits to the distributions of avalanche size and of duration and assessing the relationship between exponents. There is debate over whether these observations reflect true power laws, but within the resolution achievable from experiments, neural avalanches exhibit power laws with exponent relationships predicted from theory developed in physical systems (Perkovic et al., 1995). Avalanche criticality is not the only form of criticality observed in neural systems. Zipf's law, in which the frequency of a network state is inversely proportional to its rank, appears in systems as diverse as fly motion estimation and the salamander retina (Mora and Bialek, 2011; Schwab et al., 2014; Aitchison et al., 2016). More recently, Meshulam et al., 2019 reported various statistics of population activity in the mouse hippocampus, including the eigenvalue spectrum of the covariance matrix and the activity variance. These were found to scale as populations were 'coarse-grained' through a procedure in which neural activities were iteratively combined based on similarity. Similar observations have been reported in spontaneous activity recorded across a wide range of brain areas in the mouse (Morales et al., 2023). Simple neural network models of such data explain neither Zipf's law nor coarse-grained criticality (Meshulam et al., 2019). Even though these three forms of criticality are observed through different analyses, they may originate from similar mechanisms. Numerous studies have reported relatively low-dimensional structure in the activity of large populations of neurons (Mazor and Laurent, 2005; Ahrens et al., 2012; Mante et al., 2013; Pandarinath et al., 2018; Stringer et al., 2019; Nieh et al., 2021), which can be modeled by a population of neurons that are broadly and heterogeneously coupled to multiple latent (i.e. unobserved) dynamical variables. Using such a model, we previously reproduced scaling under coarse-graining analysis within experimental uncertainty (Morrell et al., 2021). Zipf's law has been explained by a similar mechanism (Schwab et al., 2014; Aitchison et al., 2016; Humplik and Tkačik, 2017). A single quasi-static latent variable has been shown to produce avalanche power laws, but not the relationships expected between the critical exponents (Priesemann and Shriki, 2018), while a model including a global modulation of activity can generate avalanche criticality (Mariani et al., 2021), but has not demonstrated coarse-grained criticality (Morrell et al., 2021). It is not known under what conditions the more general latent dynamical variable model generates avalanche criticality. Here, we examine avalanche criticality in the latent dynamical variable model of neural population activity. We find that avalanche criticality is observed over a wide range of parameters, some of which may be optimal for information representation. These results demonstrate how criticality in neural recordings can arise from latent dynamics in neural activity, without need for fine-tuning of network parameters. Results Critical exponents values and crackling noise We begin by defining the metrics used to quantify avalanche statistics and briefly summarize experimental observations, which have been reviewed in detail elsewhere (Plenz et al., 2021; O'Byrne and Jerbi, 2022; Girardi-Schappo, 2021). Activity is recorded across a set of neurons and binned in time. Avalanches are then defined as contiguous time bins, in which at least one neuron in the population is active. The duration of an avalanche is the number of contiguous time bins and the size is the summed activity during the avalanche. The distributions of avalanche size and duration are fit to power laws (P(S)∼S−τ for size S, and P(D)∼D−α for duration D) using standard methods (Clauset et al., 2009). Power laws can be indicative of criticality, but they can also result from non-critical mechanisms (Touboul and Destexhe, 2017; Priesemann and Shriki, 2018). A more stringent test of criticality is the 'crackling' relationship (Perkovic et al., 1995; Touboul and Destexhe, 2017), which involves fitting a third power-law relationship, S¯(D)∼Dγfit, and comparing γfit to the predicted exponent γpred, derived from the size and duration exponents, τ and α: (1) γfit=?γpred≡α−1τ−1. Previous work demonstrating approximate power laws in size and duration distributions through the mechanism of a slowly changing latent variable did not generate crackling (Touboul and Destexhe, 2017; Priesemann and Shriki, 2018). Measuring power-laws in empirical data is challenging: it generally requires setting a lower cut-off in the size and duration, and the power-law behavior only has limited range due to the finite size and duration of the recording itself. Nonetheless, there is some consensus (Shew et al., 2015; Fontenele et al., 2019; Ma et al., 2019) that even if τ and α vary over a wide range (1.5 to about 3) across recordings, the values of γfit and γpred stay in a relatively narrow range, from about 1.1 to 1.3. Avalanche scaling in a latent dynamical variable model We study a model of a population of neurons that are not coupled to each other directly but are driven by a small number of latent dynamical variables – that is, slowly changing inputs that are not themselves measured (Figure 1A). We are agnostic as to the origin of these inputs: they may be externally driven from other brain areas, or they may arise from large fluctuations in local recurrent dynamics. The model was chosen for its simplicity, and because we have previously shown that this model with at least about five latent variables can produce power laws under the coarse-graining analysis (Morrell et al., 2021). In this paper, we examine avalanche criticality in the same model. Figure 1 Download asset Open asset Latent dynamical variable model produces avalanche criticality. Simulated network is N=1024 neurons. Other parameters in Table 1. (A) Model structure. Latent dynamical variables hμ(t) are broadly coupled to neurons si(t) in the recorded population. (B) Raster plot of a sample of activity binned at 3 ms resolution across 128 neurons with five latent variables, each with correlation timescale τF=15⁢s. (C) Projection of activity into a simulated field of view for illustration. (D-F) Avalanche analysis in a network (parameters NF=5, τF=104, η=4 and ϵ=12), showing size distribution (D), duration distribution (E), and size with duration scaling (F). Lower cutoffs used in fitting are shown with vertical lines and their values are indicated in the figures. There are Nobs=42725 avalanches of size S≥Smin in this simulated dataset. Estimated values of the critical exponents are shown in the titles of the panels. Specifically, we model the neurons as binary units (si) that are randomly (Jiμ∼N(0,1)) coupled to dynamical variables hμ(t). The probability of any pattern {si}, given the current state of the latent variables, is (2) P(si|{hμ(t)})=1Z({hμ(t)})exp⁡(−η∑μ=1NFsiJiμhμ(t)−ϵsi), where the parameter η controls the scaling of the variables and ϵ controls the overall activity level. We modeled each latent variable as an Ornstein-Uhlenbeck process with the time scale τF (see Materials and methods). Thus our model has four parameters: η (input scaling), ϵ (activity threshold), τF (dynamical timescale), and NF (number of latent variables). Distributions of avalanche size and avalanche duration within this model followed approximate power laws (Figure 1C; see Materials and methods). In the example shown (NF=5, τF=104, η=4 and ϵ=12), we found exponents τ=1.89±0.02 (size) and α=2.11±0.02 (duration). Further, the average size of avalanches with fixed duration scaled as S∼Dγ, with the fitted γfit=1.24±0.02, in agreement with the predicted value γpred=1.24±0.02. Thus, our model could generate avalanche scaling, at least for some parameter choices. In the following sections, we examine how avalanche scaling depends on model parameters (NF, τF, η and ϵ; see Table 2). We first focus on two sets of simulations: one set with NF=1 latent variable, which does not generate scaling under coarse-graining (Morrell et al., 2021), and one set with NF=5 latent variables, which can generate such scaling for some values of parameters τF, η, and ϵ (Morrell et al., 2021; Table 1). Table 1 Simulation parameters for Figure 1. ParameterDescriptionValueϵbias towards silenceϵ=12ηvariance multiplierη=4.0NFnumber of latent fieldsNF=5τFlatent field time constantτ=104Nnumber of cellsN=1024 Avalanche scaling depends on the number of latent variables We analyzed avalanches from one- and five-variable simulations, each with fixed latent dynamical timescale (τF=5×103 time steps; see Table 2 for parameters). In the following sections, time is measured in simulation time steps, see Materials and methods for converting time steps to seconds. We used established methods for measuring empirical power laws (Clauset et al., 2009). The minimum cutoffs for size (Smin) and duration (Dmin) are indicated by vertical lines in Figure 2. For the population coupled to a single latent variable, the avalanche size distribution was not well fit by a power law (Figure 2A). With a sufficiently high minimum cut-off (Dmin), the duration distribution was approximately power-law (Figure 2B). Table 2 Simulation parameters for Figure 2. ParameterDescriptionValueϵbias towards silenceϵ=8 (for NF=1) or ϵ=12 (for NF=5)ηvariance multiplierη=4.0NFnumber of latent fieldsNF=1 or 5τFlatent field time constantτF=103,...105Nnumber of cellsN=1024 Figure 2 with 4 supplements see all Download asset Open asset Multiple latent variables generate avalanche scaling at shorter timescales than a single latent variable. Simulated network is N=1024 neurons. Other parameters used for simulations for this figure are found in Table 2. (A-C) Scaling analysis for one variable models where the dynamic timescale is equal to 5×103 time steps. (A) Distribution of avalanche sizes. MLE value of exponent for best-fit power law is τ=1.98 (0.02 SE), with lower cutoff indicated by the vertical line. (B) Distribution of avalanche duration. MLE value of α is 1.81 (0.02 SE). (C) Average size plotted against avalanche duration (blue points), with power-law fit (black line) and predicted relationship (yellow line) from MLE values for exponents in A and B. Gray bar on the horizontal axis indicates range, over which a power law with γ=1.72 fits the data (see Materials and methods). (D-F) Analysis of avalanches from a simulation of a population coupled to five independent latent variables where the dynamic timescale is equal to 5×103 time steps. (G) Exponents τ for avalanche size distributions across timescales for one-variable (blue) and five-variable (red) simulations. Each circle is a simulation with independently drawn coupling parameters. Simulations had to show scaling over at least two decades to be included in panels (G–J). (H) Exponents α for avalanche duration distributions for simulations in G. (I) Fitted values of γ for simulations in G. (J) Difference between fitted and predicted γ values. Five-variable simulations produce crackling over a wider range of timescales than single-variable simulations. We next assessed whether the simulation produced crackling. If so, the value γfit obtained by fitting S¯(D)∼Dγfit would be similar to γpred=α−1τ−1. In many cases, such as the one-variable example shown in Figure 2C, the full range of avalanche durations were not fit by a single power law. Therefore, we determined the largest range, over which a power law was a good fit to the simulated observations. In this case, slightly over two decades of apparent scaling were observed starting from avalanches with minimum duration slightly less than 100 time steps (Figure 2C), with a best-fit value of γfit∈[1.69,1.74]. As we did not find a power-law in the size distribution, calculating γpred is meaningless. If we do it anyway, we obtain γpred=0.83±0.03 (yellow line in Figure 2C), which clearly deviates from the fitted value of γ. Thus, for the single latent dynamical variable model (τF=5000), power-law fits are poor, and there is no crackling. The five-variable model produces a different picture. We now find avalanches, for which size and duration distributions are much better fit by power-law models starting from very low minimum cutoffs (Figure 2D–E, Figure 2—figure supplement 2). Average size scaled with duration, again over more than two decades, with γfit=1.27±0.03, which was in close agreement with γpred=1.25±0.02 (Figure 2F). Holding other parameters constant, we thus found that scaling relationships and crackling arise in the multi-variable model but not the single-variable model. Avalanche scaling depends on the time scale of latent variables Based on simulations in the previous section, we surmised that the five-variable simulation generated scaling more readily due to creating an 'effective' latent variable that had slower dynamics than any individual latent variable. We reasoned that at any moment in time, the latent variable state hμ(t) is a vector in the latent space. Because coupling to the latent variables is random throughout the population, only the length (∼NF) and not the direction of this vector matters, and the timescale of changes in this length would be much slower than τF, the timescale of each of the components hμ(t). We therefore speculated that increasing the timescale of dynamics of the latent variables should eventually lead to scaling and crackling in the single-variable model as well as the five-variable one. To examine the dependence of avalanche scaling on this timescale, we simulated one-variable and five-variable networks at fixed η and ϵ coupled to latent variables with the correlation time of their Ornstein-Uhlenbeck dynamics of τF∈[103,105] time steps, spanning from a factor of 10 faster to a factor of 10 slower than the original τF in Figure 1. Simulations were replicated five times at each combination of parameters by drawing new latent variable coupling values (Jiμ), as well as new latent variable dynamics and instances of neural firing. For simulations that passed the criteria to be fitted by power laws, we plot the fitted values of τ , α, γfit, and γfit−γpred (Figure 2G–J). Missing points are those for which distributions did not pass the power law fit criteria. In the single-variable model, best-fit exponents changed abruptly for latent variable timescale around τF=104 (Figure 2G and H), while in the five-variable model, exponents tended to increase gradually (Figure 2G and H, red). The discontinuity in the single-variable case reflected a change in the lower cutoff values in the power-law fits: size and duration distributions generated with faster latent dynamics could be fit reasonably well to a power law by using a high value of the lower cutoff (Figure 2—figure supplement 3). For time scales greater than ∼104, the minimum cutoffs dropped, and the single-variable model generated power-law distributed avalanches and crackling (Figure 2J), similar to the five-variable model. In summary, in the latent dynamical variable model, introducing multiple variables generated scaling at faster timescales. However, by slowing the timescale of the latent dynamics, the model generated signatures of critical avalanche scaling for both multi- and single-variable simulations. Avalanche criticality, input scaling, and firing threshold In the previous section, we found that a very slow single latent dynamical variable generated avalanche criticality in the simulation population. Thus, from now on, we simplify the model in order to characterize avalanche statistics across values of input scaling η and firing threshold ϵ. Specifically, we modeled a population of N=128 neurons coupled to a single quasi-static latent variable. We simulated 103 segments of 104 steps each and drew a new value of the latent variable (h∼N(0,1)) for each segment. Ten replicates of the simulation were generated at each of the combinations of η and ϵ (see Materials and methods). Almost independent of η and ϵ, we found quality power law fits and crackling. Figure 3 shows the average (across n=10 network realizations) of the exponents extracted from size (τ, Figure 3A) and duration (α, Figure 3C) distributions. At small firing threshold (ϵ=2), we do not observe scaling because the system is always active, and all avalanches merge into one. At large firing threshold ϵ and low input scaling η, we do not observe scaling because activity is so sparse that all avalanches are small. At intermediate values of the parameters, the simulations generated plausible scaling relationships in size and duration. The difference between γfit and γpred was typically less than 0.1 (Figure 4D–F), which was consistent with previously reported differences between fit and predicted exponents (Ma et al., 2019). Thus, there appears to be no need for fine-tuning to generate apparent scaling in this model, at least in the limit of (near) infinite observation time. Wherever η and ϵ generate avalanches, there are approximate power-law distributions and crackling. Figure 3 Download asset Open asset Exponents across network simulations for networks of N=128 neurons. Each parameter combination η,ϵ was simulated for ten replicates, each time drawing randomly the couplings Ji, the latent variable values, and the neural activities. Other parameters in Table 3. (A) Average across replicates for the size exponent τ. (B) Scatter plot of α vs. τ for each network replicate for parameter combinations indicated in A. Linear relationships between τ and α, corresponding to the minimum and maximum values of γfit from panel E, are shown to guide the eye. (C) Same as A, for duration exponent α. (D) Predicted exponent, γpred, derived from A and C. (E) Value of γfit from fit to S¯∼Dγ. (F) Difference between γpred and γfit. Figure 4 with 1 supplement see all Download asset Open asset Avalanches in the latent dynamical variable model with a single quasistatic variable. Parameters in Table 3. (A) Number of avalanches in simulations from Figure 3 as a function of the calculated probability of avalanches at fixed η across values of ϵ and latent variable h. Line indicates equality. (B) Probability of avalanches with η=2 across values of ϵ and h. The latent variable h is normally distributed with mean 0 and variance 1. Where the distribution of h overlaps with regions of high probability (black), avalanches occur. (C) Probability of avalanches at ϵ=8 across values of η and h. Increasing η narrows the range of h that generates avalanches. (D) Probability of avalanches at h=0 for a populations of 128 neurons (black line) and for a varying ϵ. Size distributions corresponding to simulations marked by the green and orange crosses are in E, F. (E) Example of size distribution with ϵ<ϵ0 (orange marker in D). Size cutoff is close to 100. (F) Example of size distribution with ϵ>ϵ0 (green marker in D). Size cutoff is < 10. To determine where avalanches occur, we derive the avalanche rate across values of the latent variable h. In the quasi-static model, the probability of an avalanche initiation is the probability of a transition from the quiet to an active state. Because all neurons are conditionally independent, this is Pava=Psilence(1−Psilence). Then the expected number of avalanches N^ava is obtained by integrating Pava over h at each value of η and ϵ: (3) N^ava=∫Pava(ϵ,η,h;Ji,N)p(h)dh=∫∏i(11+e−ηJih−ϵ)(1−∏i(11+e−ηJih−ϵ))p(h)dh, where p(h) is the standard normal distribution. This probability tracks the observed number of avalanches across simulations, Figure 4A. To gain an intuition for the conditions under which avalanches occur, we show two slices of the avalanche probability, at fixed η (Figure 4B) and at fixed ϵ (Figure 4C). Black regions indicate where avalanches are likely to occur. If, for a given value of ϵ and η, there is no overlap between high avalanche probability regions and the distribution of h, then there will be no avalanches. For large ϵ, avalanches occur because neurons with large coupling to the latent variable (η|Ji|>>1, recall Ji∼N(0,1)) are occasionally activated by a value of the latent variable h that is sufficient to exceed ϵ (Figure 4B). Thus, the scaling parameter η controls the value of h for which avalanches occur most frequently (Figure 4C). As ϵ decreases, avalanches occur for smaller and smaller h until avalanches primarily occur when h=0. To calculate the probability of avalanches, we must integrate over all values of h, but we can gain a qualitative understanding of which avalanche regime the system is in by examining the probability of avalanches at h=0. At h=0, the avalanche probability (see Materials and methods) is (4) Pava(ϵ,η,h=0;Ji,N)=(11+e−ϵ)N(1−(11+e−ϵ)N), which is maximized at ϵ0=−log⁡(21/N−1), independent of Ji and η. After some algebra, we find that ϵ0∼log⁡N for large N. The dependence on N reflects that a larger threshold is required for larger networks: large networks (N→∞) are unlikely to achieve complete network silence, therefore preventing avalanches from occurring. Similarly, small networks (N∼1) are unlikely to fire consecutively and thus are unlikely to avalanche. We plot Pava(ϵ,η;Ji,N,h=0) as a function of ϵ in Figure 4D. The peak at ϵ0 divides the space into two regions. For ϵ<ϵ0, a power-law is only observed in the large-size avalanches, which are rare (Figure 4E, green). By contrast, when ϵ>ϵ0, minimum size cutoffs are low (Figure 4F, orange). Both regions, ϵ<ϵ0 and ϵ>ϵ0, exhibit crackling noise scaling. If observation times are not sufficiently long (estimated in Figure 4—figure supplement 1), then scaling will not be observed in the ϵ<ϵ0 region, whose scaling relations arise from rare events. Insufficient observation times may explain experiments and simulations where avalanche scaling was not found. Inferring the latent variable Our analysis of Pava(ϵ,η,h) at h=0 suggested that there are two types of avalanche regimes: one with high activity and high minimum cutoffs in the power law fit (Type 1), and the other with lower activity and size cutoffs (Type 2). Further, when Pava drops to zero, avalanches disappear because the activity is too high or too low. We now examine how information about the value of the latent variables represented in the network activity relates to the activity type. To delineate these types, we calculated numerically ϵ∗(η), the value of ϵ, for which the probability of avalanches is maximized, and the contours of Pava (Figure 5A). Curves for ϵ∗(η) and ϵ0 and Pava=10−3 are shown in Figure 5A and B. Figure 5 Download asset Open asset Information in the neural activity about the latent variable is higher in the low-ϵ avalanche region, compared to high-ϵ avalanche or high-rate avalanche-free activity. (A) Probability of avalanche per time step across values of η and ϵ. Solid magenta curve follows ϵ∗(η), the value of ϵ maximizing the probability of avalanches at fixed η. Dashed magenta line indicates ϵ0, calculated analytically, which matches ϵ∗ at η=0. (B) Information about latent variable, calculated from maximum likelihood estimate of h using population activity. MLE approximation is invalid in the dark-blue region bounded by gray curve. Magenta line marks the maximum values of Pava, reproduced from A. Dashed black curve indicates Pava=0.001. The highest information region falls between ϵ∗(η) and the contour for Pava=0.001. (C - E) Slices of B, showing IMLE(ϵ) for η={2,5,9}. Magenta and dashed black lines again indicate ϵ∗ and Pava=0.001, respectively, as in B. Black dashed line marks the approximate boundary between the high-activity/no avalanche and the high-cutoff avalanche, and magenta line marks boundary between high-cutoff and low-cutoff avalanche regions. We expect that the more cells fire, the more information they would convey, until the firing rate saturates, and inferring the value of the latent variable becomes impossible. Figure 5B supports the prediction: generally, information is higher in regions with more activity (lower ϵ, higher η), but only up to a limit: as ϵ→0, information decreases. This decrease begins approximately where the probability of avalanches drops to nearly zero (dashed black lines, Figure 5B–E) because all of the activity merges into a few very large avalanches. In other words, the Type-1 avalanche region coincides with the highest information about the latent variable. The critical brain hypothesis suggests that the brain operates in a critical state, and its functional role may be in optimizing information processing (Beggs, 2008; Chialvo, 2010). Under this hypothesis, we would expect the information conveyed by the network to be maximized in the regions we observe avalanche criticality. However, we see that critical regions do not always have optimal information transmission. In Figure 5, the region that displays crackling noise is that where avalanches exist (Pava>0.001), which corresponds to any η value and ϵ≳3. This avalanche region encompasses both networks with high information transmission and networks with low information transmission. In summary, observing avalanche criticality in a system does not imply a high-information processing network state. However, the scaling can be seen at smaller cutoffs, and hence with shorter recordings, in the high-information state. This parallels the discussion by Schwab et al., 2014, who noticed that the Zipf's law always emerges in neural populations driven by quasi-stationary latent fields, but it emerges at smaller system sizes when the information about the latent variable is high. Discussion Here, we studied systems with distributed, random coupling to latent dynamical variables and we found that avalanche crit