Dynamic regulation of resource transport induces criticality in multilayer networks of excitable units
Abstract
Past work has shown that the function of a network of excitable units can be enhanced if the network is in the ‘critical regime’, where excitations are, one average, neither damped nor amplified. In this Letter, we show that resource transport dynamics can robustly maintain the network dynamics in the critical regime. More specifically, we consider the example of a neural network with neurons (nodes) and synapses (edges). We propose a model where synapse strengths are regulated by metabolic resources distributed by a secondary network of glial cells. We find that this twolayer network robustly preserves the critical state and produces powerlaw distributed avalanches over a wide range of parameters. In addition, the glial cell network protects the system against the destabilizing effect of local variations in parameters and heterogeneity in network structure. For homogeneous networks, we derive a reduced 3dimensional map which reproduces the behavior of the full system.
Networks of excitable units are found in varied disciplines such as social science pei:tang:zheng:2015 (), neuroscience larremore:shew:restrepo:2011 (), epidemiology newmann:2003 (), genetics newmann:2003 (), etc. Previous studies show that some aspects of network function can be optimized when the network operates in the ‘critical regime’, between low and high firing rates, as in neural networks shew:plenz:2013 (), or at the ‘edge of chaos’, between order and disorder, as in gene networks kauffman:1969 (). In particular, for neural networks, criticality results in potential information handling benefits shew:plenz:2013 (). A natural question is what mechanisms can lead such complex and distributed systems to operate in the critical regime, which typically occurs in a very small region of parameter space. In this Letter, we propose a general mechanism based on the regulation of the excitable network dynamics by a resource which enables the interactions between the excitable elements and that is transported across a secondary network. As a concrete case, we will focus on the case of neural networks, where metabolic resources that facilitate the transmission of neural excitations are transported across a secondary glial network Giaume:2010fj (); Fields:2015yq (); Suzuki:2011rt (); giaume:et:al:2010 (). We find that resourcetransport dynamics induces critical neural network dynamics characterized by neuronal avalanches following a powerlaw distribution for a wide range of model parameters. Finally, we show that in the presence of heterogeneity in local parameters or network structure, resourcetransport amongst the glial cells can be essential to maintain criticality.
Following Ref. virkar:et:al:2016 (), we consider a network model with two layers: a weighted directed neural network, and an unweighted undirected glial network which transports and regulates the supply of resources needed for the functioning of the neural layer (see Fig. 1).
Neural network dynamics The neural network consists of excitable nodes that represent neurons, labeled , and directed edges on each of which a synapse is located and labeled . We will also indicate a synapse pointing from node to node by the indices . At each discrete time step, , neuron is in either the quiescent state () or the active state (). We define as the adjacency matrix whose entry denotes the weight of the synapse on the edge from neuron to neuron at time . The state of neuron , , is updated probabilistically based on the sum of its synaptic input from active presynaptic neurons in the previous time step,
(1) 
As in Refs. larremore:et:al:2014 (); virkar:et:al:2016 (), the model transfer function probability is piecewise linear; for , for , and for .
At time each synapse is assumed to have a supply of a metabolic resource, some of which is consumed every time the presynaptic neuron, , fires. While in this paper we do not focus on a particular resource, we note that could represent various metabolites that are transported diffusively among the glial cells, such as glucose and lactate rouach:et:al:2008 (). Reflecting the increasing synaptic firing ability with increasing resource, we assume that the weight in synapse from neuron to neuron is proportional to the amount of resource in the synapse, virkar:et:al:2016 (). Finally, for simplicity, we consider only excitatory neurons and assume that there is no learning (these were considered in virkar:et:al:2016 ()). Thus, the synaptic weight changes are caused only by the dynamics of resource transport.
The second layer of our model, the unweighted and undirected glial network, consists of glial cells labeled . Each glial cell holds an amount of resource at time step . Resources diffuse between the glial cells that are connected to each other. We define a symmetric glial adjacency matrix such that entry if glial cell is connected to glial cell and otherwise. Each glial cell serves a set of synapses by supplying resource to them. Hence we define a matrix with entries if glial cell serves synapse and otherwise.
Resourcetransport dynamics: Resource diffuses between glia through their connection network (characterized by the adjacency matrix ) and between glia and the synapses they serve (via the glialneural connection network characterized by the adjacency matrix ). Our model for the evolution of the amount of resource at glial cell and the amount of resource at synapse is
(2) 
(3) 
where is the rate of diffusion between glial cells, and is the rate of diffusion between glia and synapses. Moreover, we enforce , i.e., if Eq. (3) yields , then we replace it by . The model parameter on the right hand side of Eq. (Dynamic regulation of resource transport induces criticality in multilayer networks of excitable units) denotes the amount of resource added to each glial cell at each time step (e.g., supplied by capillary blood vessels). For simplicity, we assume first that each glial cell has the same (the effect of heterogeneous values of will be discussed later). The last two terms in Eq. (Dynamic regulation of resource transport induces criticality in multilayer networks of excitable units) are the amount of resource transported to glial cell , respectively, from its neighboring glial cells and from the synapses that it serves.
The term proportional to in Eq. (3) denotes the amount of resource gained (if ) or lost (if ) from glial cell that serves synapse . If the presynaptic neuron fires at time step (), then synapse consumes an amount of resource , thus decreasing the resource at synapse by this amount.
We now describe and present the results of numerical experiments on our model, Eqs. (1)(3). Our main goal is to show that resource transport dynamics robustly regulates the operation of the neural network in the critical regime. In the neural model used here, the critical regime is characterized by the largest eigenvalue of the neural synapse matrix , , being one larremore:et:al:2014 (); virkar:et:al:2016 (). Therefore we will consider as one criterion for criticality. However, a more practical definition of criticality, applicable more generally, is a powerlaw distribution of the sizes of activity bursts, or neuronal avalanches. We will also verify that the model robustly produces powerlaw distributed neuronal avalanches.
In our numerical experiments, we will consider an ErdösRenyi network structure for both the neural and glial networks. The neural network is described using an adjacency matrix, , such that with probability we have an entry that represents a link from node to node , and otherwise we have an entry that represents the absence of such a link. At time , we set the resource at each synapse , , and draw from a uniform distribution over . By choosing the value of , we can set the initial largest eigenvalue of , i.e., , to a desired value, and test whether the subsequent evolution of the model results in .
Layer  Nodes  Adjacency  Prob. of  Mean  Mean no. 

matrix  an edge  degree  of edges  
Neural  
1000  0.05  50  50000  
Glial  
1000  0.05  50  25000 
The glial network is undirected and unweighted and its adjacency matrix is given by a matrix, , such that with probability we have an entry that represents an undirected link between nodes and . Motivated by experiments azevedo:et:al:2009 () showing , for specificity, we take the number of glia and neurons to be equal, . Consistent with this, and the additional experimental finding that all incoming synapses of a given neuron are served by the same glial cell halassa:et:al:2007 (), we further assume that each glial cell serves a unique neuron. We set the initial resource for each glial cell to be equal to an uniform value, (note that may be different from ). For all numerical experiments we use the values of the parameters shown in Table 1 and assume, for simplicity, that the entries of matrices and are independent of each other. Unless mentioned otherwise, the parameters for resourcetransport dynamics are set as , , and . Although we chose these values somewhat arbitrarily, we will show that our findings are fairly robust to changes in model parameters, and that we can determine the parameter ranges for which our model results in from a simplified 3dimensional map.
In the first experiment, we show that starting with different initial conditions , the resourcetransport dynamics causes the system to selforganize to the critical state corresponding after a transient period. In Fig. 2 we show for the three different initial conditions (black circles, red triangles, and blue squares, respectively). In the three cases, approaches and subsequently remains close to (this will be quantified in Fig. 4).
As discussed above, we are interested in whether the dynamics of the neuronal network reproduces experimental signatures of critical behavior, in particular powerlaw distributed avalanches of activity. To do this, following larremore:et:al:2014 (), we define a measure of activity, , and define an avalanche as the excursion of activity above a threshold , i.e., for , and for ). We define the size of the avalanche as , the number of spikes (excitations) over a single excursion.
To investigate the robustness of our model to changes in parameters, we fix and vary and logarithmically roughly from to keeping the ratio constant. Using the threshold , we calculate avalanche size distributions, , for each parameter setting. We then fit a powerlaw model using standard techniques clauset:shalizi:newman:2009 () based on maximumlikelihood estimation (MLE) and a hypothesis test that generates a value using the KolmogorovSmirnov (KS) statistic. Since we have finitesize effects, in addition to the lower size cutoff used in clauset:shalizi:newman:2009 (), we introduce an upper cutoff , i.e., we test the plausibility of a powerlaw model where we condition on avalanche sizes in a range Langlois:2014ij (); Shew:2015pb (). We accept as plausible power laws only those distributions for which this range spans at least three decades.
In Fig. 3 we show size distributions for various values of . The blue curves correspond to those distributions for which we get a plausible powerlaw fit (under level of significance) with such that (exponents range from to ), and the red curves correspond to distributions for which a powerlaw fit is rejected. We find that there is a large range of values of the parameters (small values of and for which the system results in powerlaw distributed avalanches, but that for some parameter choices (larger values of and ) the distribution of avalanche sizes no longer satisfies our KolmogorovSmirnov powerlaw test.
In order to gain an understanding of the mechanisms that lead to the critical regime and to determine conditions on the model parameters that result in critical behavior, we assume a homogeneous network structure for both the neural and the glial networks, i.e., the outdegree of each neuron in the neuronal network is . In this case, we can approximate the largest eigenvalue of the neural network adjacency matrix by the mean degree, , where we assumed that the weights are uncorrelated with the resource and defined the average resource in glial cells as (recall we are taking ). Summing Eq. (Dynamic regulation of resource transport induces criticality in multilayer networks of excitable units) over and Eq. (3) over and , we obtain
(4)  
(5) 
The above equations are coupled to the average activity , which is a stochastic variable determined by Eq. (1). In order to obtain a tractable map, we model in two different ways: in the first one, we neglect stochastic effects and use the deterministic approximation:
(6) 
This approximation is based on the fact that, except for values of close to or , the expectation of calculated from Eq. (1) is . This approximation neglects the nonlinear effects that keep below , and thus should be interpreted only as a guide to determine the fixed points and their stability in the limit of vanishing stochastic effects [i.e., when the expected number of terms in the sum in Eq. (1) is large]. We refer to Eqs. (4), (5) and (6) as the 3D map without noise. While this approximation neglects the effects of noise and is valid only close to the fixed point, it is useful since the stability properties of the fixed point underlie the robustness of the critical state to changes in model parameters. A more realistic model for includes a stochastic noise term and a mechanism to enforce :
(7) 
Here, is a Gaussian noise term with zero mean and standard deviation , as estimated in larremore:et:al:2014 (), while represents an external stimulus taken to be with probability and otherwise, introduced to prevent from decaying to zero. Effectively, this stimulus excites one neuron every time step with probability . We refer to Eqs. (4), (5) and (7) as the 3D map with noise. This variant of the map is useful to predict the behavior of the macroscopic variable of the full model. As an example, in Fig. 2 we show with dotted () and dashed () lines the predictions of the evolution of obtained from the 3D map with noise. The predictions agree very well with direct simulations of the full model and, being a three dimensional map, are much more computationally efficient. The resource dynamics, although not shown, also is well predicted by the 3D map with noise.
While the 3D map with noise allows us to reproduce the behavior of the full model, the simplicity of the 3D map without noise allows us to derive parameter constraints that must be satisfied for a stable critical state. In particular, the system of Eqs. (4)(6) has the fixed point
(8) 
The critical state is a fixed point of the deterministic map. Its stability is determined by whether the eigenvalues of the Jacobian of the map (4)(6) evaluated at the fixed point (8) are inside the complex unit circle. Applying the RouthHurwitz criterion, we find that the fixed point is stable (which we interpret as robustness of the critical state) if
(9)  
(10)  
(11)  
(12) 
In addition, since , we have the additional inequality
(13) 
which represents the constraint that the amount of resource supplied to a glial cell per time step can not exceed the amount that can be consumed at the synapses it serves.
To demonstrate the usefulness of these inequalities, we verify which of the curves in Figure 3 correspond to parameters which satisfy these inequalities. Parameters that satisfy (don’t satisfy) the inequalities correspond to distributions which follow (don’t follow) a powerlaw. We note, however, that in other cases, even though the fixed point is linearly unstable, can oscillate so close about that the avalanche size distribution is still a plausible power law. Conversely, for moderately large networks finite size effects can lead to fluctuations of about that are significant enough to cause deviations from power law behavior, as noted in Ref. LABEL:quasicriticality_reference. To illustrate this, we plot in Fig. 4 the standard deviation in the time series of , which we denote as , as a function of (as in Fig. 2, we take and ). Since we find that fluctuates about , serves as a measure of the deviation of the system from . The red triangles correspond to simulations of the full model, the blue circles to the 3D map with noise, and the dashed line to the 3D map without noise. The shaded grey region indicates parameter values for which the linear stability analysis predicts the fixed point to be unstable. We observe that the 3D map with noise captures the deviations from very well until these become large slightly past the onset of instability, i.e., approximately when . The 3D map without noise, neglecting fluctuations, fails to capture the small deviations from that occur before the onset of instability. To relate these findings with the distribution of avalanche sizes, we indicate with an arrow the value of such that values yield avalanche size distributions that have plausible powerlaw fits with exponent close to (the same information that was used to color the curves in Fig. 2). The map without noise thus predicts roughly the onset of instability and, correspondingly, of avalanche size distributions that are not powerlaw distributed.


So far, our results have been independent of resource transport in the glial network [e.g., does not appear in Eqs. (4) and (5), and numerical experiments show that our results so far hold if we set ], and thus one might wonder if local supply and consumption of resources might be sufficient to lead the system towards a global critical state. However, in the next numerical experiment we show that, when there are heterogeneities (in the network structure, in the supply and consumption rates and , etc), the diffusion of resources between glial cells can compensate for this effects and prevent destabilization of the critical state. We consider the particular case of heterogeneous source rates, where now each glial cell has its own . As an example, we draw the from a Gaussian distribution with mean and standard deviation , so that approximately of them do not satisfy the inequality (13), . In the absence of resource transport, resource accumulates in these glial cells and the associated synapses, bringing the network to the supercritical state. However, when resource is allowed to diffuse, the critical state is maintained. This is shown in Fig. 5 (a), where we plot as a function of time with (red triangles) and (blue circles). Figure 5 (b) shows that the total resource accumulates indefinitely when , and reaches a steady state when .
To summarize, in this Letter we have shown that resourcetransport dynamics can stabilize the dynamics of excitable units so that they operate at a critical state characterized by powerlaw distributed avalanche sizes. We found that for a large range of parameters, the system selforganizes to a critical state that is characterized by powerlaw distributed avalanche sizes with an exponent value near the characteristic exponent found in various experimental studies. We found that resource transport dynamics protects the system against the destabilizing effects of parameter heterogeneity. Finally, we emphasize that, although we presented our model in the context of neuronal networks, where the resource represents metabolites, our results could be applicable to other networks of excitable elements whose interaction efficacy depends on the availability of a shared resource.
References
 [1] S. Pei, S. Tang, and Z. Zheng. PLoS ONE, 10, e0124848 (2015).
 [2] D. B. Larremore, W. L. Shew, and J. G. Restrepo. Phys. Rev. Lett., 106, 058101 (2011).
 [3] M. E. J. Newmann. SIAM Rev., 45, 167 (2003).
 [4] W. L. Shew and D. Plenz. Neuroscientist, 19, 88 (2013).
 [5] S. A. Kauffman. J. Theo. Biol., 22, 437 (1969).
 [6] L. Roux D. Holcman C. Giaume, A Koulakoff and N. Rouach. Nature Reviews Neuroscience, 11(2):87–99
 [7] D. H. Woo R. D. Fields and P. J. Basser. Neuron, 86(2):374–386
 [8] O. Bozdagi G. W. Huntley R. H. Walker P. J. Magistretti A. Suzuki, S. A. Stern and C. M. Alberini. Cell, 144(5):810–823
 [9] C Giaume, A. Koulakoff, L. Roux, D. Holcman, and N. Rouach. Nat. Rev. Neurosci., 11, 87 (2010).
 [10] Y. S. Virkar, W. L. Shew, J. G. Restrepo, and E. Ott. Phys. Rev. E, 94, 042310 (2016).
 [11] D. B. Larremore, W. L. Shew, E. Ott, F. Sorrentino, and J. G. Restrepo. Phys. Rev. Lett., 112, 138103 (2014).
 [12] N. Rouach, A. Koulakoff, V. Abudara, K. Willecke, and C. Giaume. Science, 322, 1551 (2008).
 [13] F. A. C. Azevedo et al. J. Comp. Neurol., 513, 532 (2009).
 [14] M. M. Halassa, T. Fellin, H. Takano, JH. Dong, and P. G. Haydon. J. Neurosci., 27, 6473 (2007).
 [15] A. Clauset, C. R. Shalizi, and M. E. J. Newmann. SIAM Rev., 51, 661 (2009).
 [16] D. Cousineau D. Langlois and J.P. Thivierge. Physical Review E, 89(1):012709, 2014.
 [17] J. Pobst Y. Karimipanah N. C. Wright W. L. Shew, W. P. Clawson and R. Wessel. Nature Physics, 11(8):659–663