The role intrinsic statistical fluctuations play in creating avalanches – patterns of complex bursting activity with scale-free properties – is examined in leaky Markovian networks. Using this broad class of models, we develop a probabilistic approach that employs a potential energy landscape perspective coupled with a macroscopic description based on statistical thermodynamics. We identify six important thermodynamic quantities essential for characterizing system behavior as a function of network size: the internal potential energy, entropy, free potential energy, internal pressure, pressure, and bulk modulus. In agreement with classical phase transitions, these quantities evolve smoothly as a function of the network size until a critical value is reached. At that value, a discontinuity in pressure is observed that leads to a spike in the bulk modulus demarcating loss of thermodynamic robustness. We attribute this novel result to a reallocation of the ground states (global minima) of the system's stationary potential energy landscape caused by a noise-induced deformation of its topographic surface. Further analysis demonstrates that appreciable levels of intrinsic noise can cause avalanching, a complex mode of operation that dominates system dynamics at near-critical or subcritical network sizes. Illustrative examples are provided using an epidemiological model of bacterial infection, where avalanching has not been characterized before, and a previously studied model of computational neuroscience, where avalanching was erroneously attributed to specific neural architectures. The general methods developed here can be used to study the emergence of avalanching (and other complex phenomena) in many biological, physical and man-made interaction networks.

Networks of noisy interacting components arise in diverse scientific disciplines. Here, we develop a mathematical framework to study the underlying causes of a bursting phenomenon in network activity known as avalanching. As prototypical examples, we study a model of disease spreading in a population of individuals and a model of brain activity in a neural network. Although avalanching is well-documented in neural networks, thought to be crucial for learning, information processing, and memory, it has not been studied before in disease spreading. We employ tools originally used to analyze thermodynamic systems to argue that randomness in the actions of individual network components plays a fundamental role in avalanche formation. We show that avalanching is a spontaneous behavior, brought about by a phenomenon reminiscent to a phase transition in statistical mechanics, caused by increasing randomness as the network size decreases. Our work demonstrates that a previously suggested balanced feed-forward network structure is not necessary for neuronal avalanching. Instead, we attribute avalanching to a reallocation of the global minima of the network's stationary potential energy landscape, caused by a noise-induced deformation of its topographic surface.

Funding: This material is based upon work supported by the National Science Foundation ( http://www.nsf.gov ) under Grant Number CCF-1217213. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of the National Science Foundation.

Introduction

An important problem in many scientific disciplines is understanding how extrinsic and intrinsic factors enable a complex physical system to exhibit a bursting behavior that leads to avalanching [1], [2]. Avalanching is a form of spontaneous behavior characterized by irregular and isolated bursts of activity that follow a scale-free distribution typical to systems near criticality. In the brain, this mode of operation is thought to play a crucial role in information processing, memory, and learning [2]–[6].

Although avalanche dynamics have been extensively studied in vitro [7] and in vivo [8], [9] for cortical neural networks, it is not clear what causes avalanching. A recent in silico attempt to address this issue [10] was based on approximating the dynamics of a Markovian model of nonlinear interactions between noisy excitatory and inhibitory neurons by Gaussian fluctuations around the macroscopic system behavior using the linear noise approximation (LNA) method of van Kampen [11]. This led to the conclusion that the cause of neural avalanches is a balanced feed-forward (BFF) network structure. We argue here that the Gaussian approximation used to arrive at this conclusion is not appropriate for studying avalanching, thus leading to deficient results. As a consequence, understanding the underlying causes of avalanching in silico is still an open problem.

To address this challenge, we introduce a theoretical framework that allows us to examine the role of intrinsic noise in inducing critical behavior that leads to avalanching. Although the idea that noise may induce avalanching has been proposed more that a decade ago [12], our framework leads to a novel understanding of the underlying causes of avalanching in a particular class of complex networks. We focus on a general Markovian network model, which we term leaky Markovian network (LMN), with binary-valued state dynamics. These dynamics are described by a time-dependent probability distribution that evolves according to a well-defined master equation [13] (see Methods for details). It turns out that a LMN is a continuous-time stochastic Boolean network model with a state-dependent asynchronous node updating scheme (we provide details in Text S1). LMNs can model a number of natural and man-made systems of interacting species, such as genetic, neural, epidemiological, and social networks.

Recent work has clearly demonstrated the importance of stochastically modeling physical systems using Markovian networks. The main reason is that intrinsic noise produced by these networks may induce behavior not accounted for by deterministic models [14]–[17]. Examples of such behavior include the emergence of noise-induced modes, stochastic transitions between different operational states, and “stabilization” of existing modes.

In this paper, we study the effect of intrinsic noise on avalanching by using a LMN model. We do so by employing the notion of potential energy landscape [13], [18], [19] and by establishing a connection between statistical thermodynamics and the kinetics of bursting. We quantify the landscape by calculating logarithms of the ratios between the stationary probabilities of individual states and the stationary probability of the most probable state. To reduce computational complexity, we follow a coarse graining approach that transforms the original LMN model into another (non-binary) LMN model with appreciably smaller state-space. To accomplish this task, we partition the nodes of the LMN into homogeneous subpopulations and characterize system behavior by using the dynamic evolution of the fractional activity process, which quantifies the fraction of active nodes (nodes with value 1) in each subpopulation. Moreover, we parameterize the LMN in terms of the network size , where is the net number of nodes in the network and is a normalizing constant such that can be approximately considered to be continuous-valued. We refer the reader to the Methods section for details.

The behavior of the fractional activity process is fundamentally affected by . In general, the strength of stochastic fluctuations (intrinsic noise) in the activity process may be thought of as the probability of moving uphill on a fixed potential energy surface, which decays exponentially with increasing . At sufficiently large network sizes , the LMN operates around a ground state of the potential surface located at a fixed point predicted by the macroscopic equations associated with the LNA method, which we assume to be unique, nonzero, and stable (see Methods and Text S1 for details). In this case, a new mode of operation is introduced in the system, as the network size decreases, in the form of a potential well in the topographic surface of the energy landscape, located at the inactive state . This is a “noise-induced” mode, since it appears at small network sizes at which the fractional activity process is subject to appreciable intrinsic fluctuations.

We show in this paper that noise-induced deformation of the stationary potential energy landscape is the underlying cause of avalanching in LMNs. For sufficiently large network sizes, the potential energy landscape can be approximated by a quadratic surface centered at . In this case, the LMN operates within the potential well associated with this mode, except for rare and brief random excursions away from that mode. As a consequence, the fractional activity process will fluctuate in a Gaussian-like manner around the macroscopic mode. At smaller network sizes, the fractional activity process is characterized by a bistable behavior between the macroscopic and noise-induced modes, spending most time within the potential well associated with the macroscopic mode, at which the potential energy surface attains its global minimum, while occasionally jumping inside the potential well associated with the noise-induced mode at . As a consequence, the fractional activity dynamics take on a bursting behavior characterized by long periods of appreciable activity followed by short periods of minimal (almost zero) activity. When the network size decreases further, the noise-induced mode becomes the main stable operating point (i.e., the point at which the potential energy surface attains its global minimum), whereas the macroscopic mode becomes shallower and eventually disappears. In this case, the system is trapped within the potential well associated with the noise-induced mode, except for random and brief excursions away from that mode. As a consequence, the fractional activity process will still exhibit bursting, but now characterized by long periods of minimal (almost zero) activity followed by short bursts of appreciable activity.

Thermodynamic analysis reveals critical behavior in LMNs (we provide details in the Methods section and Text S1). By employing a number of statistical thermodynamic quantities, such as internal and free potential energies, entropy, internal pressure, pressure and bulk modulus (inverse compressibility), we effectively summarize the stochastic behavior of a LMN as its size decreases to zero. We also use these summaries to quantify network robustness and the stability of a given state. In agreement with the classical theory of phase transitions, the previous thermodynamic quantities evolve smoothly as a function of until a critical network size is reached. At this size, a discontinuity is observed in the system pressure, which produces a spike in the bulk modulus demarcating loss of thermodynamic robustness. Critical behavior is caused by reallocation of the ground states (global minima) of the potential energy landscape due to noise-induced deformation of its topographic surface. In particular, observed critical behavior produces two distinct phases: one in which the fixed point predicted by the macroscopic equations associated with the LNA method constitutes the ground state of the potential energy landscape and one in which the ground state is reallocated to the noise-induced mode at . We conclude that avalanching is a complex mode of operation that dominates system dynamics at near-critical and subcritical network sizes due to deformations of the potential energy landscape as the network size decreases to zero, caused by appreciable levels of intrinsic noise.

It is important to mention here that our work provides a novel stochastic perspective to the well-known phenomenon of self-organized criticality (SOC) [20]; i.e., the spontaneous emergence of critical behavior without tuning system parameters whose values are influenced by external factors. We approach SOC from the perspective of nonlinear Markovian dynamics and directly associate self-organization properties of a complex network with the existence of a unique probability distribution at steady-state, which leads to a unique stationary potential energy landscape. Our work demonstrates that SOC may emerge as a consequence of two interweaved adaptive processes that may take place on separate timescales [21], [22]: a short timescale convergence to dynamic equilibrium (stationarity), during which the topological structure of a neural network subsystem is kept relatively fixed (is quasistatic), and longer timescale alterations in topological properties that lead to changes in the size of the network. For example, the neural activities modeled by our LMNs occur on a short timescale compared to the timescale of neural development, where it is well known that programmed cell death plays a major role [23]. This large reduction in the number of neurons could presumably serve to bring neural subsystems within proximity of their critical sizes in accordance with their underlying connectivity structures. On the other hand, an adaptive interplay between the short timescale neural dynamics and the quasistatic topological network structure may lead to a longer timescale topological self-organization that enlarges or contracts the neural subsystem with an objective to keep the system robustly close to criticality [22]. We show that such long timescale alterations can result in a spontaneous reallocation of the ground states in the stationary potential energy landscape due to a noise-induced deformation of its topographic surface. Our approach associates SOC with observable stochastic multistability, which is directly related to the phenomenon of phase transition, thus bridging the gap between “self-organized” and “classical” criticality (see also [24]).

Finally, we would like to point out that, in the neural network literature, it is commonly said that avalanching occurs in the supercritical regime near the critical point. On the other hand, we show in this paper that avalanching occurs in a subcritical regime near the critical point. To avoid confusion, the reader must keep in mind that these statements do not contradict each other. Correctly using the terms “subcritical” and “supercritical” depends on the parameter employed to take a system from one regime to the other. Contrary to existing works that study criticality in terms of functional parameters (e.g., in terms of the firing rate of all neurons in a neural network), we study in this paper criticality in terms of a structural parameter, the system size, which is inversely related to the strength of intrinsic noise. In this case, avalanching occurs at network sizes smaller but near a critical value, which forces us to use the aforementioned terminology.