Stochastic quorum percolation and noise focusing in neuronal networks

The spontaneous activity of developing dissociated neuronal cultures is characterized by bursts of high-frequency collective activity followed by silent periods, with little activity [1,2]. Burst emergence can be explained by a process of nucleation and front propagation, where activity first increases in a small region of the culture, and then quickly propagates throughout the remaining of the system [3]; this rapid propagation depletes neurons of available neurotransmitters, and the system remains almost silent until it recovers and another burst can develop. In any typical culture, multiple nucleation sites can coexist, i.e., zones of high nucleation probability, where a burst can develop and propagate. The presence of specific nucleation sites and its spatial heterogeneity arise due to noise focusing [3], a symmetry-breaking mechanism that amplifies the quenched disorder in the network wiring, modifying the transport properties of the spontaneous activity. The resulting coarse-grained physical picture is that of noise flowing anisotropically through the two-dimensional space where the network is embedded. The combination of small inhomogeneities in network structure, together with the non-linear dynamics of the neurons, result in the noise being amplified while it concentrates in specific locations, i.e., the nucleation sites.

While this scenario has been recently understood at a mesoscopic level [4], the mechanisms of dynamical and topological amplification [3] that are responsible for this phenomenon are difficult to grasp at a quantitative level. To gain insight into this central point and in particular to elucidate the role of metric connectivity correlations (which play a major role in the behavior of neuronal cultures [5]) we present here a minimal statistical model, called stochastic quorum percolation (SQP), that captures the generic features of the noise focusing mechanism by extending the previous theory of quorum percolation [6] into a fully dynamical model. SQP is also closely related to other theoretical frameworks such as directed, compact percolation, branching processes and cellular automata [7].

Previously the concept of quorum percolation (QP) was introduced [6] to describe the collective dynamics of neuronal networks under external stimulation [8]. In these experiments, the network connectivity is weakened with different chemical synaptic blockers and an external current is applied to the system to study its response. QP is an extension of percolation theory in which a node requires a minimum quorum of m simultaneous inputs to become active. For a given initial fraction f of active nodes, the quorum condition is checked iteratively until no more nodes can be activated. QP is characterized by a discontinuous phase transition, where a small increase in f generates a large response in the size of the final active fraction Φ. QP is similar to bootstrap [9] and k-core percolation [10], although several differences exist, specially in the presence of disorder [11,12].

In the absence of any external stimuli, however, neuronal systems are spontaneously active, and the propagation of neuronal activity is often described in terms of avalanches [13], where any neuron firing can trigger a firing on its neighbours with probability p₁. By extension, a neuron receiving k inputs will fire with probability $p_k=1-(1-p_1)^k$. In QP we can define a similar expression, where a neuron fires if, and only if, it receives m inputs, $p_k=\Theta(m-k)$, where m is called the quorum threshold. However, real neurons are driven by noise, and always have a finite spontaneous firing rate. In the stochastic quorum percolation (SQP) model we assume that the internal noise of the neuron is characterized by a Poisson process with a given rate λ (a shot noise). This noise model is inspired by minis (spontaneous miniature post synaptic currents), which in sparse cultures have the same strength as evoked currents [14,15]. This description is also equivalent to considering an unobserved population projecting onto the observed network and firing with Poisson statistics. Hence, for simplicity, we assume that each discharge of our shot noise has the same effect as receiving one input from another neuron. Within this description, the probability that a neuron fires spontaneously within a time window $\Delta t$ is

$$\begin{eqnarray} &&p_0 = e^{-\lambda\Delta t}\sum_{i=m}^\infty\frac{(\lambda\Delta t)^i}{i!}={\cal P}(m,\lambda\Delta t), \end{eqnarray} \tag{ 1 }$$

where ${\cal P}$ is the regularized gamma function, and we choose $\Delta t$ as the characteristic integration time of synaptic currents $\Delta t\approx 20\ \text{ms}$. Hence, the spontaneous firing frequency of a neuron is $\omega_0\approx p_0/\Delta t$. Accordingly, an isolated neuron spontaneously fires when it accumulates the required quorum of m shots from a Poisson process with rate λ within $\Delta t$. The case where a neuron is receiving k external inputs within $\Delta t$ is equivalent to lowering the firing quorum from m to m –k, i.e.,

$$\begin{eqnarray} &&p_k = {\cal P}(m-k, x), \end{eqnarray} \tag{ 2 }$$

where we define $x\equiv\lambda\Delta t$ as the noise strength. In the limit of large m and weak noise $x\to0$, we recover the expression of the original QP model $p_k=\Theta(m-k)$.

At its essence, SQP is a dynamical process where spontaneous activations can induce other activations throughout the network, and, given enough inputs (quorum), generate a deterministic response. This minimal set of features can be related, one by one, to classical models of non-equilibrium phase transitions: its dynamical nature is that of directed percolation [16]; the spontaneous activations are related to the presence of external fields in non-equilibrium systems, e.g., the spontaneous generation of new particles in the pair-contact process [17]; and the quorum condition is akin to the presence of an upper absorbing state in compact directed percolation and in the Domany-Kinzel model [18,19].

The temporal evolution of SQP can be described as a Markov process within a network. The network is characterized by its size N and its adjacency matrix A, whose entries $A_{ij}=1$ denote the existence of a connection from node i to node j. And the state of each node in the system S_i is either active $S_i=1$ or inactive $S_i=0$. The evolution of the system is characterized by its transition probabilities

$$\begin{eqnarray} &&\textrm{Pr} \left({S_i(t+1)=1}\,|\,{\vec S(t)}\right)={\cal P}(m-k_i(t),x), \end{eqnarray} \tag{ 3 }$$

where $\vec S(t)=\left(S_1(t),\ldots,S_N(t)\right)$ and $k_i(t)=(A^\mathsf{T}\vec S(t))_i=\sum_{j}A_{ji}S_j(t)$, i.e., the number of active inputs of i at time t.

The model can be analyzed at the mean-field level in the limit of a maximally entropic infinite random graph with fixed in and out degree distributions, $p(K_i)=p(K_o)=p(K)$, i.e., N → ∞, $\langle CC\rangle\to0$ (average clustering coefficient). We characterize the system by the fraction of active nodes at time t, $\Phi(t)=N^{-1}\sum_iS_i(t)$ and look for the steady state solution $\Phi(t+1)=\Phi(t)$. Note from eq. (3) that taking into account all possible combinations of active inputs we can write

$$\begin{eqnarray} &&\Phi = \sum_K p(K)\sum_{n=0}^K\binom{K}{n}{\cal P}(m-n,x)\Phi^n(1-\Phi)^{K-n},\qquad \end{eqnarray} \tag{ 4 }$$

which can be solved numerically ¹ . Note that ${\cal P}$ has a sigmoid shape, and for all relevant cases, eq. (4) has either one or three solutions (see fig. 1). Without loss of generalization, we will assume a connectivity distribution $p(K)=\delta(K-\bar K)$ from now on.

Zoom In Zoom Out Reset image size

For $m<\bar K$ (which is met for any real system, i.e., each node has more connections than the quorum needed to activate), and low levels of noise x, eq. (4) has three solutions, two stable and one unstable, even without noise $(x=0)$. The system can be found in two regimes: one with very low activity, where $\Phi\sim p_0$, and the system is essentially inactive, and one where all the nodes are active, $\Phi=1$, signaling bursting behavior. There exists a critical noise value $x^*$, for which the low-activity solution and the unstable solution merge, hence defining a discontinuous phase transition. For any $x>x^*$, the whole system is active. This transition to $\Phi=1$ is equivalent to the emergence of a burst in real neuronal networks, as we will see later on. $\Phi=1$ is an absorbing state, and the inclusion of a mechanism similar to short-term synaptic depression is needed to get the inverse transition (from bursting to non-bursting). Note that if the system size is finite, for $x<x^*$ there exists a finite probability per unit time of transition between the states, so the system will eventually reach the absorbing state, i.e., it will burst in a finite time.

In the low-activity regime, the steady state of the system is characterized by the presence of cascades of activity, or avalanches, defined as a temporal sequence of causally connected activations. Following the formalism of branching processes and studies of criticality in similar systems [13,20,21], we can define a local branching ratio as

$$\begin{equation} \sigma = \bar K\sum_{n=0}^{\bar K-1}\binom{\bar K-1}{n}\phi^n(1-\phi)^{\bar K-n-1}{\cal P}(m-n-1, x),\quad \end{equation} \tag{ 5 }$$

for a network with fixed connectivity $\bar K$, where σ denotes the average number of nodes that a given active node will activate in a posterior time step. If we impose $p_0=0$, which corresponds to having separation of time-scales between random activations and activity propagation, $\sigma=1$ marks the presence of sustained activity, i.e., the active phase. There exists a critical level of noise, $x^\dagger$, where $\sigma=1$ and $x^\dagger$ marks a continuous, second-order phase transition. For $x<x^\dagger$, there is an absorbing state with $\Phi=0$, the inactive phase, and any cascade of activity will eventually die out. From eq. (5) we obtain in the inactive phase that $\sigma=\bar K{\cal P}(m-1, x)$, hence $x^\dagger$ is easily obtained from the relation $1/\bar K={\cal P}(m-1, x^\dagger)$, i.e., $p_1=1/\bar K$, as in most kinds of percolation on a Bethe lattice [22]. In fact, the model belongs to the same universality class as directed percolation.

The presence of spontaneous activity, $p_0\neq0$, however, destroys the lower absorbing state (x > 0) and consequently the transition [17]. Regardless, one can still define a critical noise value $x^\dagger$ with maximum susceptibility $\chi=\textrm{d}\Phi/\textrm{d}p_0$ that defines a nonequilibrium Widom line [23,24]. The mean-field phase diagram for a given network connectivity is shown in fig. 2, for the particular case of $p_0=0$.

Zoom In Zoom Out Reset image size

Although the previous mean-field approximation is only valid for infinite random graphs, a similar phase diagram can be constructed for real networks. However, several differences exist. In a finite network, even for the case of a random graph, given that the system is driven by noise, it can reach the upper absorbing state $(\Phi=1)$ even for low levels of noise. It only needs a fluctuation in the number of active nodes Φ big enough to reach the unstable solution of eq. (4) (see fig. 1). In this situation the lifetime of the fluctuations of the system, as well as the temporal correlation in the activity of the nodes, depend strongly on the precise network topology, as we will now see.

To show the effects network of topology on the dynamics we will restrict ourselves to the study of the metric networks presented in ref. [3] that reproduce the structure of dissociated neuronal cultures. These networks are obtained as follows: neuronal bodies (nodes) are placed at random in a substrate, usually a square of size L with periodic boundary conditions until a desired density is reached $(\rho=100\text{--}1000\ \text{neurons/mm}^2)$. Then, an axon is grown from each body as a biased random walk that results in almost straight axons with a characteristic length $\ell_a\sim1\ \text{mm}$ and, finally, the dendritic tree of a neuron is modeled as a circular area around its body with a characteristic radius of $r_d\sim0.15~\textrm{mm}$. A connection between two neurons $i\to j$ is created with probability α whenever the axon of neuron i intersects with the dendritic tree of neuron j.

This procedure generates networks with in- and out-degree distributions characterized by the distributions of dendritic and axonal sizes, respectively (we usually choose a Gaussian one for the first and a Rayleigh for the latter), and with a clustering coefficient CC that depends on α and ρ. Typically $\langle CC\rangle\sim0.25$. These networks present connectivity correlations that decay exponentially with a characteristic length of $\ell_c=0.26\ \text{mm}$, and can be seen as random graphs at small distances $r<0.15\ \text{mm}$ but are highly directional at long distances due to the particular morphology of neurons, with all their output connections found in a narrow area around their axons. To study the impact that the underlying metric correlations have on the network dynamics, for each constructed network we generate a maximally entropic conjugate network with the constraint of keeping the same in- and out-degree distributions via a random swap of links between unconnected neuronal pairs.

As we have discussed previously, in finite systems with spontaneous activity (which is the case for any real system), the continuous phase transition does not really exist and the discontinuous transition is not marked by a sharp line and instead is associated to a characteristic first-passage time, that we will call the bursting time, i.e., the average time it takes the system to reach the bursting phase, i.e., the upper absorbing state. In fig. 3 we characterize the transition by the smallest spontaneous firing rate $\omega_0=p_0/\Delta t$ required for the system to reach the bursting phase within a fixed, arbitrarily long time (5000 s). This figure illustrates the important effect of network topology in the dynamics, where the presence of clustering and higher-order correlations in the network structure favors the presence of the bursting phase, shifting the transition line to smaller spontaneous firing rates.

Zoom In Zoom Out Reset image size

This effect is more dramatically observed in the dependence of the bursting time on the spontaneous firing rate for a given quorum threshold m (see fig. 4). For high spontaneous firing rates, the bursting time is insensitive to the topology. For low rates, however, it has a strong dependence, increasing at a much faster rate in the random case than in the metric one. The range of frequencies where the system can reach the bursting phase in a realistic time is largely expanded in the metric network.

Zoom In Zoom Out Reset image size

The strong dependence of the bursting time on the metric correlations of the network is an essential feature of the noise focusing mechanism, as pointed out in ref. [3]. The non-linear behavior of the induced firing probability function ${\cal P}$ results in lower requirements of simultaneous inputs to reach the bursting phase in regions of high input connectivity and clustering, for which the noise is effectively stronger. In a network embedded in a two-dimensional space, once the bursting phase is reached within a finite region the active phase will spread to the whole system (as long as $\bar K/2>m$ [9]), regardless of its size. An inactive neuron near the neighborhood to the active phase will have, on average, half of its input connections active, since they will come from neurons inside the active phase. Hence as long as $\bar K/2>m$, neurons in the neighborhood of the active phase will become active, extending the boundary of the active phase even further.

The SQP model provides an appropriate description of neuronal dynamics in terms of percolation concepts, providing a simplified and unified framework to account for both bursting dynamics and neuronal avalanches (power-law statistics of cascades of induced firings). The two types of dynamics coexist within the same model, and bursting behavior can appear with either critical or supercritical avalanches.

For simplicity, and with the aim of identifying general mechanisms that go beyond the case of neuronal systems, we have not considered the effects of synaptic depression, facilitation, inhibition or variable synaptic weights, but these are straightforward additions to the model, since all of them can be associated to changes in ${\cal P}$. Note that with the addition of synaptic depression or inhibition, the model can easily exhibit self-organized criticality, as in ref. [25].

Neuronal avalanches are characteristic of many neuronal systems, and have been described in a myriad of systems within the theory of branching processes [13,26–29], and recently with directed percolation [30]. But even more ubiquitous is the presence of bursts, which are a direct consequence of the integrate-and-fire dynamics. For sufficient activity the system will exhibit a discontinuous transition, not just super-critical avalanches, as predicted by simpler models. Critical neuronal avalanches are believed to be a desirable feature [31], and there have been many attempts to describe mechanisms by which neuronal systems self-organize towards that behavior [21,25,32]. System-wide bursts, on the other hand, are often linked to epileptic behavior and run-away excitation [33], and neuronal systems have to implement mechanisms to prevent, stop or reduce them. In particular, networks with balanced excitatory and inhibitory connectivity often show smaller and heterogeneous bursting activity, both in in vitro and in vivo preparations.

The behavior presented here is not exclusive of neuronal systems, as any process running on a network of integrate-and-fire elements with non-linear summation of probabilities should operate in a similar way. Our model could easily be adapted to disease-spreading processes for diseases where the probability of contagion can quickly rise or without separation of time-scales [24]; or to rumor spreading with the illusion of truth effect [34], where the bombardment of false information in a short time window can lead to false beliefs and further spreading.

We thank J. Davidsen and R. Williams-Garcia for useful discussions. We thank financial support of MINECO under project PID2019-108842GB-C21 and from Generalitat de Catalunya under project 2017-SGR-1061.