Neutral theory and scale-free neural dynamics

TL;DR

This paper proposes that scale-free neural avalanches can arise from neutral drift dynamics, rather than criticality, challenging the common assumption that such behavior indicates the brain operates at a phase transition point.

Contribution

It introduces a model showing scale-invariant neural avalanches emerge from neutral dynamics, not criticality, emphasizing the role of causal information in analyzing neural activity.

Findings

Neural avalanches can be scale-free without criticality.

Neutral drift explains power-law distributed avalanches.

Causal information is crucial for accurate avalanche analysis.

Abstract

Avalanches of electrochemical activity in brain networks have been empirically reported to obey scale-invariant behavior --characterized by power-law distributions up to some upper cut-off-- both in vitro and in vivo. Elucidating whether such scaling laws stem from the underlying neural dynamics operating at the edge of a phase transition is a fascinating possibility, as systems poised at criticality have been argued to exhibit a number of important functional advantages. Here we employ a well-known model for neural dynamics with synaptic plasticity, to elucidate an alternative scenario in which neuronal avalanches can coexist, overlapping in time, but still remaining scale-free. Remarkably their scale-invariance does not stem from underlying criticality nor self-organization at the edge of a continuous phase transition. Instead, it emerges from the fact that perturbations to the system exhibit a neutral drift --guided by demographic fluctuations-- with respect to endogenous spontaneous activity. Such a neutral dynamics --similar to the one in neutral theories of population genetics-- implies marginal propagation of activity, characterized by power-law distributed causal avalanches. Importantly, our results underline the importance of considering causal information --on which neuron triggers the firing of which-- to properly estimate the statistics of avalanches of neural activity. We discuss the implications of these findings both in modeling and to elucidate experimental observations, as well as its possible consequences for actual neural dynamics and information processing in actual neural networks.

Equations

$$\begin{array}[]{ccc}I&\overset{\epsilon}{\longrightarrow}&A_{M(t)+1}\\ A_{k}+I&\overset{\lambda}{\longrightarrow}&A_{k}+A_{k}\\ A_{k}&\overset{\mu}{\longrightarrow}&I\end{array}$$

$${\dot{\rho }}_k=\left(\lambda (1-\rho )-\mu \right){\rho }_k+\sqrt{\frac{1}{N}\left(\lambda (1-\rho )+\mu \right){\rho }_k}{\xi }_k(t),$$

$${\dot{\rho }}_k=-\frac{\mu }{\lambda -\mu }ϵ{\rho }_k+\sqrt{\frac{\mu }{N}\left(2-\frac{ϵ}{\lambda -\mu }\right){\rho }_k}{\xi }_k(t).$$

$${\dot{\rho }}_k=\sqrt{{\rho }_k}\text{leavevmode}{\xi }_k(\hat{t}),$$

$${\dot{V}}_i=-\frac{V_i-V_r}{RC}+\sum _k\frac{I_{{\text{e}}_i}^k(t)}{C}+\frac{1}{C}\sum _{\begin{array}{c}{i}^{\prime }\in n.n.(i)\\ j,k\end{array}}\Theta (p_r{U}_{i^{\prime }j}(t_{s_{i^{\prime }}}^k)-{\zeta }_{i^{\prime }j}^k)I_{{\text{in}}_{i^{\prime }}}^k(t),$$

$${\dot{U}}_{ij}=\frac{1-U_{ij}}{{\tau }_R}-\sum _k{U}_{ij}\Theta (p_r-{\zeta }_{ij}^k)\delta (t-t_{s_i}^k).$$

Full text

\JournalInfo\Archive\PaperTitle

Neutral theory and scale-free neural dynamics \AuthorsMatteo Martinello1, Jorge Hidalgo1, Serena di Santo2,3, Amos Maritan1, Dietmar Plenz4 and Miguel A. Muñoz2,*

\KeywordsNeural avalanches — Criticality — Neutral theory – Generic scale invariance \AbstractAvalanches of electrochemical activity in brain networks have been empirically reported to obey scale-invariant behavior –characterized by power-law distributions up to some upper cut-off– both in vitro and in vivo. Elucidating whether such scaling laws stem from the underlying neural dynamics operating at the edge of a phase transition is a fascinating possibility, as systems poised at criticality have been argued to exhibit a number of important functional advantages. Here we employ a well-known model for neural dynamics with synaptic plasticity, to elucidate an alternative scenario in which neuronal avalanches can coexist, overlapping in time, but still remaining scale-free. Remarkably their scale-invariance does not stem from underlying criticality nor self-organization at the edge of a continuous phase transition. Instead, it emerges from the fact that perturbations to the system exhibit a neutral drift –guided by demographic fluctuations– with respect to endogenous spontaneous activity. Such a neutral dynamics –similar to the one in neutral theories of population genetics– implies marginal propagation of activity, characterized by power-law distributed causal avalanches. Importantly, our results underline the importance of considering causal information –on which neuron triggers the firing of which– to properly estimate the statistics of avalanches of neural activity. We discuss the implications of these findings both in modeling and to elucidate experimental observations, as well as its possible consequences for actual neural dynamics and information processing in actual neural networks.

Introduction

The introduction by Kimura in 1968 of the neutral theory –hypothesizing that most evolutionary change is the result of genetic drift acting on neutral alleles [1]– caused much debate and a revolution in the way population genetics and molecular evolution were understood. In a similar endeavor, Hubbell proposed that most of the variability in some ecological communities could be ascribed to neutral dynamics of similar species which expand or decline as a result of stochasticity [2, 3]. Neutral theories have in common that they neglect any a priori intrinsic difference between coexisting individuals, regardless of their “species” (allele, tree,…) type, implying that the dynamics is purely driven by random demographic effects. For instance, the introduction of a novel species within an established population triggers a random cascade of changes, or “avalanche”, which –as a result of the implicit neutrality– does not have an inherent tendency to neither shrink nor to expand at the expenses of others. This marginal-propagation process generates scale-free avalanches, which resembles critical ones even if the system is not necessarily posed at the edge of a phase transition [3, 4] (a brief mathematical summary of neutral theory can be found in SI Appendix S1). Neutral models have been successfully employed to explain the emergence of scale-free distributions in (i) epidemic outbreaks with neutral microbial strains [4], (ii) viral-like propagation of neutral memes [5], (iii) the evolution of the microbiome [6], and (iv) the renewal of the intestinal epithelium from neutral stem cells [7]. Could neutral theory be applied to neural dynamics of the brain? And, in particular, could it explain the emergence of neuronal avalanches reported for spontaneous activity?

The human brain has a special feature that is common to all mammalians: it is endogenously active; i.e. cascades of electrochemical activity at multiple timescales spontaneously pervade its dynamical state even in the absence of any apparent stimuli or task. Mounting evidence suggests that such an endogenous activity is not random, but structured, and it contributes significantly to stimulus-related responses, being essential to brain functioning. Specifically, spontaneous, spatiotemporal bursts of neural activity were reported to appear in the form of avalanches [8], whose sizes $s$ and durations $t$ are distributed as $P_{s}(s)\sim s^{-\tau}\leavevmode\nobreak\ \leavevmode\nobreak\ {\cal{F}}(s/s_{c})$, and $P_{t}(t)\sim t^{-\alpha}\leavevmode\nobreak\ \leavevmode\nobreak\ {\cal{G}}(t/t_{c})$, respectively, where $\tau\approx 3/2$ and $\alpha\approx 2$ are critical exponents similar to those of an unbiased branching process [9, 10, 11], $\cal{F}$ and $\cal{G}$ are scaling functions and $s_{c}$ and $t_{c}$ are system-size dependent cut-offs obeying finite-size scaling [12]. Similar results have been obtained both in vitro and in vivo, as well as for different tissues, preparation types, experimental techniques, and animal species (see e.g. [13, 14, 15, 16, 17, 18, 19]). Remarkably, signs of scale-invariance have been reported to vanish under abnormal circumstances such as under the influence of modified pharmacological conditions, under anesthesia, or in pathological conditions [20]. We refer to [21, 22, 23, 24, 25] for overviews and discussions on the state of the art. Taken together, these observations suggest that scale-free avalanches are a generic feature of spontaneous activity in cortical tissues, suggesting that they stem from an underlying critical phenomenon (see however [26]), and this conclusion seems to back the hypothesis that biological computing systems might operate at the edge of phase transitions [27, 28, 29], providing them with optimal transmission and storage of information, exquisite sensitivity to signals, and a number of other important functional advantages [30, 31].

Scale-free distributed events or bursts of “activity” such as earthquakes, vortex avalanches in superconductors, and Barkhaussen noise are common place in Nature (see e.g. [32, 33]) and are often ascribed to their underlying dynamics being poised at a critical point. The paradigm of “self-organized criticality” was developed to explain how and why natural systems could self-tune to the vicinity of critical points [32, 34, 35]; scale-free distributed avalanches turn out to be the fingerprint of critical points of a phase transitions into quiescent (or “absorbing”) states [36, 37]. Despite the success and conceptual beauty of this framework, not all scale-invariant episodes of activity can be ascribed to underlying criticality [38, 39]; for instance, power-law distributed excursion sizes and times are generated by unbiased random walks [40], self-organization to the edge of a discontinuous phase transition [41] and, as discussed above, neutral dynamics [4, 7].

In this paper we explore the possibility that empirically observed neural avalanches could be scale-free as a result of an underlying neutral dynamics –i.e. that each single event of activity is indistinguishable from others and can marginally propagate through the network– alternatively to being self-organized to the edge of a phase transition. This is, we explore whether scale-free avalanches could stem from the neutral competition for available space of activity generated from different sources or stimuli. We put forward a subtle but important difference between such causal avalanches and existing empirical evidence, and discuss how neutral patterns of activity –i.e. coexisting causal avalanches of many different shapes, sizes and durations– could be exploited by real neural systems for efficient coding, optimal transmission of information and, thus, for memory and learning [42],

Results

Computational model and its phenomenology

Using a model of leaky integrate-and-fire neurons regulated by synaptic plasticity, Millman, et al. [43] were able to capture the empirical observation of bistability in cortical networks, i.e. two well differentiated stable patterns of cortical activity, called Up and Down states (see e.g. [44, 45] and Refs. therein). Briefly, the model consists of $N$ leaky integrate-and-fire excitatory neurons forming a directed random Erdős-Rényi network with average connectivity $K$. Neurons integrate synaptic inputs from other neurons and fire action potentials, which rapidly deplete associated synaptic resources. These latest recover at a slow time scale, thereby regulating the overall level of activity in the network (see MM). The model can be tuned by controlling e.g. its average synaptic strength. For weak synaptic strengths, a quiescent phase with very low levels of activity, the Down state, exists, whereas a second, stable state with high firing rates, the Up state, emerges for large synaptic strengths (see Fig. 1A). For intermediate strengths, spontaneous fluctuations allow for rapid Up and Down states alternations (see Fig. 1B). This phenomenology – which can also be reproduced by keeping synaptic strength fixed and varying the synaptic recovery time or some other parameter of the model– corresponds to a discontinuous phase transition (see Fig. 1 and Fig. 1 in [43]) and therefore lacks the critical point characteristic of continuous transitions. Remarkably, when tracking cascades of neuronal firing based on participating neurons, i.e. causal avalanches (see below), the model was shown to exhibit scale-free distributions of sizes and durations during Up-states, with associated exponents $\tau\approx 3/2$ and $\alpha\approx 2$, i.e. the hallmark of actual neuronal avalanches. Accordingly, the authors considered the Up state as “self-organized critical”, in contrast to the Down state which was “subcritical” with causal cascades that were not scale-free [43]. Given that critical dynamics emerge at continuous phase-transitions, the presence of scale-invariant Up state avalanches in the absence of any such transition in this model is puzzling, and prompted us to identify possible alternative mechanisms for the emergence of scale-free avalanches.

“Causal” avalanches

Following [43], we tracked causal cascades/avalanches.111Here, we use indistinctly the terms “avalanche” and “cascade”. Each one is initiated when an external input depolarizes a neuron’s membrane potential above its threshold to fire an action potential, it unfolds as the membrane potential of a subsequent neuron surpasses its threshold as a result of a synaptic input from an existing cascade member, and stops when this does not happen. The size of a cascade is the total number of action potentials it triggered, while the cascade duration is the timespan between its initiation and the time of its last action potential [43]. Avalanches were analyzed separately for Up and Down states in a network with $N=3000$ neurons, using different values of the external firing rate, $f_{e}$; in particular we analyzed the slow driving limit $f_{e}\rightarrow 0$. Our results are in perfect agreement with the phenomenology found in [43]: cascades in the Down state do not exhibit scale invariance but instead have a characteristic scale (see SI Appendix S2). In contrast, cascades during Up-states distribute in size and duration according to power-laws with exponent close to $\tau=3/2$ and $\alpha=2$, respectively (see Fig. 2A). As already observed in [43], these results are quite robust, do not depend on how deep into the Up state –i.e. how far from the transition point– simulations are run, nor do they change upon introducing inhibitory neurons (see SI Appendix S3).

“Time-correlated” avalanches from time binning

A key point of the previous analysis is that causal information between activation events (i.e. “who triggers who”) is essential to define avalanches. However, in empirical analyses it is not clear whether events occurring nearby in time –usually ascribed to the same avalanche in statistical analyses– are actually causally connected or not. The standard approach, that has been successfully used in the analysis of experimental data, where causal information of event propagation is typically not accessible [8, 13], consists in defining cascades from a series of discrete supra-threshold events, by choosing a discrete time bin $\Delta t$. An avalanche is defined as a sequence of successive time windows with at least one event in each that is preceded and ended by an empty bin. In principle, one could expect different scaling relations when varying the time window $\Delta t$, as was demonstrated in analyses of empirical data [8]. For comparison, and following [8, 13], we take $\Delta t$ to be equal to the average inter-event interval (IEI), defined as the average time interval between successive events.222Even though the IEIs can vary for different experimental situations, size and duration distributions have been claimed to exhibit universal behavior with exponent $\tau\approx 3/2$ and $\alpha\approx 2$, respectively, provided that data are binned using the IEI. Using this binning procedure in timeseries from the computational model, we find that cascade duration and size distributions obtained from Up states are exponentially distributed with a characteristic scale, showing no signs of scale-invariant behavior (see Fig. 2B). Distributions did not change qualitatively for different values of $\Delta t$. Thus, in the model of Millman et al., cascades based on temporal proximity differ significantly from cascades based on causal information. This finding is in contrast to the established scale-free avalanche distributions that emerge from experimental data based on temporal proximity.

In the model of Millman et al., causal avalanches can (and do) coexist in time (see Fig. 1C and D); thus, the temporal proximity approach does necessarily fail to uncover true causal avalanches. Thus, our observations, together with the lack of a continuous phase transition, question the origin of true scale invariance within Up states and its actual relationship with experimental findings. To shed light on this problem, in the next section we analyze a minimal model which captures the main ingredients for activity propagation, showing that the observed scale-free causal avalanches in the model of Millman et al. stems from an underlying neutral dynamics [4, 7].

Neutral (causal) avalanches in a minimal model for

activity propagation

In archetypical models of activity propagation such as the contact process, directed percolation and the susceptible-infected-susceptible model [10, 46], “active” sites propagate activity to their nearest neighbors or become de-activated at some rates. As a result, depending on rate values, there exist a quiescent and an active phase, as well as a critical point separating them [10, 46]; avalanches started from a single initial event exhibit scale invariance only at criticality (see SI Appendix S4), and if they are triggered at a sufficiently slow rate, they do not overlap.

In contrast, within the framework of neutral dynamics, multiple avalanches can propagate simultaneously. The difference between critical and neutral avalanches can be vividly illustrated by considering a variant of the contact process, consisting of many different but equivalent “species”. It can be studied arbitrarily far from the phase transition to explore the statistics of causal avalanches. More specifically, we consider a fully-connected network with $N$ nodes that can be either active ($A$) or inactive ($I$). At every time, each single active site is assigned to a unique individual avalanche/species $k$ (the one from which it derives) and labeled by $A_{k}$. The dynamics is as follows: i) a new avalanche, with a new label, is initiated by the spontaneous activation of an inactive site at small driving rate $\epsilon$; ii) active sites propagate the activity to neighboring inactive places at rate $\lambda$, and iii) active sites become inactive at rate $\mu$. This is equivalent to the following set of reactions for $k=1,...,M(t)$:

[TABLE]

where $M(t)$ is the total number of avalanches triggered up to time $t$. This dynamical process is neutral (or symmetrical) among avalanches/species as parameter rates do not depend on label $k$ (see SI Appendix S1 for an extended discussion on neutral theories). The duration (resp. size) of an avalanche $k$ is the time elapsed (resp. total number of activations) between its spontaneous generation and the extinction of its label. Observe that different avalanches can coexist (all the most in the active phase) and that the total number of coexisting avalanches can vary in time. The state of the system is determined by $M(t)$ and the number of $k-$type active sites, $n_{k}(t)$, or, equivalently, the corresponding density $\rho_{k}(t)=n_{k}(t)/N$. The total density of active sites is defined as $\rho(t)=\sum_{k=1}^{M(t)}\rho_{k}(t)$. Importantly, the system of Eq. (1) is nothing but the standard contact process (with a non-vanishing rate for spontaneous activation $\epsilon$) if avalanche labels are ignored. Therefore, in the limit $\epsilon\rightarrow 0$, it exhibits a continuous phase transition for the total activity density at the critical point given by $\lambda_{c}=\mu$ [10, 46].

We performed computer simulations of the dynamics described by Eq. (1) by means of the Gillespie algorithm [47] in a fully-connected network of size $N=10^{4}$. Parameter values are chosen for the system to lie well inside the active phase, $\lambda=2$, $\mu=1$ (i.e. $\lambda=2\lambda_{c}$), and $\epsilon$ taking small values such as $10^{-1},10^{-2},10^{-3}$ and $10^{-4}$. Typical timeseries for individual avalanches, $\rho_{k}$, as well as for the total activity, $\rho$, are depicted in Fig. 3A. Observe that the steady-state overall density (gray color) coincides, on average, with that of the contact process in the infinite size limit, $\rho^{*}\simeq(1-\mu/\lambda)+\epsilon\mu/(\lambda(\lambda-\mu))$ (see MM). On the other hand, individual avalanches (colored curves in Fig. 3A) experience wild fluctuations as a function of time. The statistics of avalanches is illustrated in Fig. 3B revealing that avalanche sizes and durations are power-law distributed with exponents $\tau=3/2$ and $\alpha=2$ in the limit of small spontaneous activation rate $\epsilon\rightarrow 0$. Remarkably, scale-free avalanches appear all across the active phase, $\lambda>\lambda_{c}$ (see SI Appendix S4).

Analytical approach

To shed light on this result, we study analytically this simplified model in the large network-size limit. Starting from the master equation associated to Eq. (1), performing a system-size expansion for large but finite system sizes [48], the dynamics of a just-created avalanche is well-described by the following equation:

[TABLE]

with the initial condition $\rho_{k}=1/N$, and where $\xi_{k}(t)$ represents a zero-mean Gaussian white noise of unit variance (to be interpreted in the Itô sense [48]). If the system is very large, and when the spreading rate lies within the active phase ($\lambda>\mu$), the total activity density exhibits very small fluctuations, remaining quite stable around the steady-state value, as illustrated by the gray-colored timeseries in Fig. 3A.

To understand this variability, let us assume $\rho(t)\simeq\rho^{*}$ in Eq. (2) and $\rho_{k}\ll\rho$ (which is good approximation for large system sizes where many avalanches coexist); thus,

[TABLE]

In the limit of small driving, $\epsilon\rightarrow 0$, the deterministic/drift term in Eq. (3) vanishes, and the dynamics of avalanche $k$ can be simply written as:

[TABLE]

where for simplicity in the notation, a factor $2\mu/N$ has been reabsorbed into the time scale $\hat{t}$. Eq. (4) represents a freely-moving random-walk with demographic fluctuations, and –as further discussed in SI Appendix S1– it describes the evolution of a species density in any neutral-type of dynamics. In other words, once an avalanche starts, its statistics is entirely driven by neutral demographic fluctuations, regardless of the distance to the critical point333Not surprisingly, Eq. (4) corresponds also to the mean-field description of an unbiased branching process [9].. Furthermore, the avalanche exponents associated with this neutral, noise-driven, dynamics are $\alpha=2$ and $\tau=3/2$. Actually, the previous reasoning holds all across the active (Up) phase; on the other hand, in the quiescent (Down) state, the steady state activity $\rho^{*}$ goes to [math], as the deterministic driving force in Eq. (2) is negative, leading to subcritical avalanches, as indeed reported in [43].

Thus, a simple approach allowed us to explicitly show that neutral dynamics among coexisting dynamically-indistinguishable avalanches leads to scale-free distributions all across the active (Up) phase, i.e. arbitrarily far away from the edge of the phase transition. The same conclusion extends to the model of Millman et al. even if detailed analytical calculations for such a case are more difficult to perform.

Discussion

A remarkable observation –that has elicited a great deal of interest– is that neural activity in the brain of mammals, including humans, occurs in the form of neuronal avalanches consisting of outbursts of neural activity intervened by periods of relative quiescence, across many resolution scales in a robust way [8, 21]. For in vitro studies of relatively small networks it seems plausible to assume that events occurring during one of such outbursts are causally connected, so that activity emerges at some location and transiently propagates through the network, causing a cascade of co-activations. However, there is no clear empirical validation that this is actually the case; diverse causally-connected cascades could, in principle, occur simultaneously, hindering their experimental discrimination as individual avalanches. Obviously, the situation is much more involved in large neural networks as analyzed in vivo as diverse scales of resolution, e.g. from local field potential measurements, magneto-encephalography, functional magnetic resonance imaging, etc. There is no known empirical procedure to actually disentangle causal influences, nor to discern whether different causal cascades of activations overlap (as they probably do in functional brains). Thus, in the absence of a better indicator, events of activity are customarily clustered together as individual avalanches, relying on a criterion of temporal proximity.

It thus remains to be fully elucidated what is the true nature of scale-free avalanches in actual neural systems. To shed light on this, here we scrutinized the most commonly referred model –introduced by Millman and coauthors [43]– justifying the emergence of power-law distributed avalanches in networks of integrate-and-fire neurons with synaptic plasticity. First of all, we reproduced the findings in [43], and confirmed that the model exhibits two different phases in parameter space, an Up-state characterized by large average firing rates and a Down-one with small firing, separated by a discontinuous phase transition. We carefully analyzed the dynamics within the active phase, and corroborated that diverse avalanches can coexist, and that their sizes and durations are scale-free (with exponents, $3/2$ and $2$, respectively) if and only if precise information on which neuron triggers the firing of which –which is accessible in computational models– is used to identify (causal) avalanches. On the other hand, a different analysis –which is the one customarily applied to empirical data– based on defining avalanches through a time-binning procedure, blind to detailed causal information between activation events, does not reveal any trace of scale-freedom in avalanche distributions.

These observations naturally pose two important questions. First, if this model is not self-organized to the edge of a phase transition, where do the computationally-reported scale-free (causal) avalanches within this model stem from? And second, does this model constitute a faithful representation of actual neural dynamics, including the experimentally observed scale-invariant avalanches?

To answer the first question we designed a simplified dynamical model with an overall phenomenology very similar to that of the model of [43]: i.e. it exhibits scale-invariant causal avalanches all along its active phase, regardless of the distance to a phase-transition point (which actually can be either a continuous or a discontinuous one depending on model details). This simplified model –a variant of the contact process with many different types of active particles— allowed us to uncover that scale-invariant avalanches within the active phase stem from the neutral dynamics among diverse coexisting (causal) avalanches. In particular, if new seeds of activity are injected at a very slow rate in a system with recurrent background activity (i.e. in its active phase) each one does not have a net drift toward contracting or expanding in the background of recurrent activity in which it unfolds; its dynamics just follows demographic fluctuations, much as in neutral theories of population genetics. Moreover, the branching ratio is equal to unity, and causal avalanches are power-law distributed (as in the unbiased branching processes) without the model being posed at the edge of a phase transition. In summary, the observed scale-invariance in a well-accepted computational model for neuronal dynamics as well as in a simplified model stems from the neutrality or symmetry between diverse competing cascades of causally-related events which coexists in a background of recurrent activity.

In what respects the second question above, it might occur that the discussed computational model does not reproduce all the phenomenology of actual neural dynamics in real networks. For instance, activity exhibits clear temporal clustering (so that measured power-laws disappear when times are reshuffled [8]) and, as we have shown, this fact is mostly lacking in the model of Millman et al.. This drawback was overcome in a more recent and detailed computational model including many additional neuro-physiologically realistic ingredients (such as e.g. inhibitory plasticity) which exhibits temporal clustering of activity together with scale-free avalanches [49]. In this case, avalanches are also measured employing causal information so that scale-invariance is likely to stem from underlying neutrality, rather than criticality. It would be highly desirable to have a study of purely time-binned avalanches in this type of approach, allowing to compare them with causal ones. From a broader perspective, more complete computational models and/or analyses allowing to scrutinize the possible emergence and interplay between neutrality and criticality are highly needed.

Finally, the main question that remains to be answered is: given that various types of functional advantages are ascribed to criticality, do these same advantages still exist if neuronal scale-free avalanches turn out to be the consequence of underlying neutrality rather than of the tuning to the edge of a phase transition? While we do not have a definite answer to this, we speculate that this type of power-law distributed coexisting causal avalanches could play a fundamental role in neural functioning. In particular, there are known biological mechanisms, such as learning rules, that take into account causal information (i.e. which neuron triggers the firing of which); a well-documented example is “synaptic timing dependent plasticity” (STDP) [42] by which synaptic weights are either reinforced or weakened depending on the relative spike timing between the pre- and post-synaptic neuron. STDP has been found to stabilize the dynamics of neural networks and to maintain reproducible patterns of causal neuronal avalanches [49]. Thus, patterns of activity, generated by neutral dynamics, and consisting on mostly-non-interacting scale-free avalanches could be stored and stabilized or “ingrained” by such a mechanism, allowing the network to spontaneously generate a large set of attractors and a broad dynamical repertoire, in a similar way in which cellular diversity –stemming from underlying neutral dynamic of stem cells– entails functional advantages in epithelial tissues [7]. These speculative ideas need to be much more carefully scrutinized, and we plan to do so in forthcoming work.

Summing up, some of the existing confusion surrounding different types of scale-invariance in neural activity can be rationalized in the framework of neutral theories, posing new and fascinating questions that may contribute to clarify the criticality hypothesis in the cortex and its implications for function and learning.

Materials and Methods

Model for neural dynamics

The model of Millman et al. [43] consists of a population of $N$ leaky integrate-and-fire excitatory neurons which are randomly connected in a directed graph to, on average, other $K$ neurons in the population (i.e. forming a Erdős-Rényi network [50]). External inputs, $I_{\text{e}}^{k}(t)$, are Poisson-distributed with rate $f_{e}$ and internal inputs, $I_{\text{in}}^{k}(t)$, are generated from spiking neurons in the network ($k$ accounts for the input number). Both internal and external currents are modeled by exponentials functions of amplitude $w_{\text{e/in}}$ and characteristic time $\tau_{s}$, $I_{\text{e/in}_{i}}^{k}(t)=w_{\text{e/in}}\exp(-(t-t_{s_{i}}^{k})/\tau_{s})$, where $t_{s_{i}}^{k}$ represents the corresponding spiking time of neuron $i$. Each individual neuron $i$ is described by a dynamical variable $V_{i}$ representing its membrane potential. When this value reaches a threshold value $\theta$, the neuron spikes and it may open –with probability $p_{r}$– each of its $n_{r}$ associated release sites per synapse, inducing a postsynaptic current. After spiking, the membrane potential is reset to the resting potential value, $V_{r}$, for a refractory period $\tau_{rp}$, during which its dynamics is switched-off. Synaptic depression is implemented by means of a dynamical “utility” variable $U_{ij}(t)\in[0,1]$ (for neuron $i$ and release site $j$), which modulates the release probability $p_{r}\rightarrow U_{ij}p_{r}$. The membrane potential obeys the following equation:

[TABLE]

where $R$ is the membrane resistance $C$ its capacitance, $k$ is the spike number, $i^{\prime}$ runs over presynaptic neurons linking to $i$, and $j^{\prime}$ over its release sites; $\zeta_{i^{\prime}j^{\prime}}^{k}$ is a uniform random number in $[0,1]$ and $\Theta(x)$ the Heaviside step function. On the other hand, the synaptic utility $U_{ij}$ is set to [math] immediately after a release and recovers exponentially to $1$ at constant rate, $\tau_{R}$:

[TABLE]

As equations (5) and (6) are linear during successive events, they can be integrated exactly, which allowed us to implement both synchronous (or clock-driven) and asynchronous (or event driven) methods [51], leading to essentially indistinguishable results. When not specified, model parameters were taken as in [43]: $K=7.5$, $n_{r}=6$, $R=2/3$ G$\Omega$, $C=30$ pF, $V_{r}=-70$ mV, $\theta=-50$ mV, $w_{\mathrm{e}}=95$ pA, $w_{\mathrm{in}}=50$ pA, $p_{r}=0.25$, $\tau_{rp}=1$ ms, $\tau_{s}=5$ ms and $\tau_{R}=0.1$ s. We also studied versions of the model including inhibitory couplings, but this did not alter the main conclusions (see SI Appendix S3).

Steady state density in the minimal model for activity propagation

Neglecting fluctuations from finite size effects, the dynamics of the total density of activity for the process described by Eq. (1) becomes deterministic in the limit $N\rightarrow\infty$: $\dot{\rho}=(\lambda(1-\rho)-\mu)\rho+\epsilon(1-\rho)$. Its stationary solution, $\dot{\rho}=0$, is $\rho^{*}=(\lambda-\mu-\epsilon-\sqrt{4\epsilon\lambda+(\lambda-\mu-\epsilon)^{2}})/2\lambda$, that, up to first order in $\epsilon$, can be written as $\rho^{*}\simeq 1-\mu/\lambda+\epsilon\mu/(\lambda(\lambda-\mu))$ in the active phase, $\lambda>\mu$, and $\rho^{*}\simeq\epsilon/(\mu-\lambda)$ in the quiescent one, $\lambda<\mu$.

Acknowledgments

We are grateful to the Spanish-MINECO for financial support (under grant FIS2013-43201-P; FEDER funds) and to P. Villegas, P. Moretti, S. Suweis and S. Vassanelli for extremely useful discussions and comments.