e058101

Predicting Criticality and Dynamic Range in Complex Networks: Effects of Topology Daniel B. Larremore,1,* Woodrow L. Shew,2 and Juan G. Restrepo1 1Department of Applied Mathematics, University of Colorado, Boulder, Colorado 80309, USA 2National Institutes of Health, National Institute of Mental Health, Bethesda, Maryland 20892, USA (Received 30 July 2010; published 31 January 2011) The collective dynamics of a network of coupled excitable systems in response to an external stimulus depends on the topology of the connections in the network. Here we develop a general theoretical approach to study the effects of network topology on dynamic range, which quantifies the range of stimulus intensities resulting in distinguishable network responses. We find that the largest eigenvalue of the weighted network adjacency matrix governs the network dynamic range. When the largest eigenvalue is exactly one, the system is in a critical state and its dynamic range is maximized. Further, we examine higher order behavior of the steady state system, which predicts that networks with more homogeneous degree distributions should have higher dynamic range. Our analysis, confirmed by numerical simulations, generalizes previous studies in terms of the largest eigenvalue of the adjacency matrix. DOI: 10.1103/PhysRevLett.106.058101 PACS numbers: 87.18.Sn, 05.45.�a, 87.19.lj Numerous natural [1,2] and social [3] systems are accu- rately described as networks of interacting excitable nodes. The collective dynamics of such excitable networks often defy naive expectations based on the dynamics of the single nodes which constitute the network. Recent experi- ments in neural networks [4] suggest that their dynamic range (the range of stimuli over which there is significant variation in the collective response of the network) is maximized in a critical regime in which neuronal ava- lanches [5] occur, confirming earlier theoretical predic- tions [2]. It has been argued [2,4] that this critical regime occurs when the expected number of excited nodes pro- duced by one excited node is one. However, this criterion is invalid for networks with broad degree distributions [6,7]. A general understanding of how dynamic range and criticality depend on network structure remains lacking. In this Letter, we present a unified theoretical treatment of stimulus-response relationships in excitable networks, which holds for diverse networks including those with random, scale-free, degree-correlated, and as- sortative topologies. As a tractable model of an excitable network, here we consider the Kinouchi-Copelli model [2], which consists of N coupled excitable nodes. Each node i can be in one of m states xi. The state xi ¼ 0 is the resting state, xi ¼ 1 is the excited state, and there may be additional refractory states xi ¼ 2; 3; . . . ; m� 1. At discrete times t ¼ 0; 1; . . . the states of the nodes xti are updated as follows: (i) If node i is in the resting state, xti ¼ 0, it can be excited by another excited node j, xtj ¼ 1, with probability Aij, or indepen- dently by an external process with probability �. The network topology and strength of interactions between the nodes is described by the connectivity matrix A ¼ fAijg. In this model, � is considered the stimulus strength. (ii) The nodes that are excited or in a refractory state, xti � 1, will deterministically make a transition to the next refractory state if one is available, or otherwise return to the resting state (i.e., xtþ1i ¼ xti þ 1 if 1 � xti < m� 1, and xtþ1i ¼ 0 if xti ¼ m� 1). An important property of excitable networks is the dynamic range, which is defined as the range of stimuli that is distinguishable based on the system’s response F. Following [2], we quantify the network response with the average activity F ¼ hfit, where h�it denotes an average over time and ft is the fraction of excited nodes at time t. To calculate a system’s dynamic range, we first determine a lower stimulus threshold �low below which the change in the response is negligible and an upper stimulus thresh- old �high above which the response saturates. Dynamic range (�), measured in decibels, is defined as � ¼ 10log10�high=�low. To analyze the dynamics of this system, we denote the probability that a given node i is excited at time t by pti. For simplicity, we will consider from now on only two states, resting and excited (m ¼ 2) [8]. Then, the update equation for pti is ptþ1i ¼ ð1�ptiÞ � �þ ð1��Þ � 1�Y N j ð1�ptjAijÞ �� ; (1) which can be obtained by noting that 1� pti is the probability that node i is resting at time t, and the term in large parentheses is the probability that it makes a transition to the excited state. We assumed that the events of neighbors of node i being excited at time t are statisti- cally independent. As noted before [3,9–11], this approxi- mation yields good results even when the network has a non-negligible amount of short loops. In Ref. [2], the response F was theoretically analyzed as a function of the external stimulation probability � using a mean-field approximation in which connection strengths were considered uniform, Aij ¼ �=N for all i; j. It was shown that at the critical value � ¼ 1, the PRL 106, 058101 (2011) P HY S I CA L R EV I EW LE T T E R S week ending 4 FEBRUARY 2011 0031-9007=11=106(5)=058101(4) 058101-1 � 2011 American Physical Society network response F changes its qualitative behavior. In particular, lim�!0F ¼ 0 if �< 1 and lim�!0F > 0 if �> 1. In addition, the dynamic range of the network was found to be maximized at � ¼ 1. The parameter � is defined in Refs. [2,4] as an average branching ratio, written here as � ¼ 1N P i;jAij ¼ hdini ¼ hdouti, where dini ¼ P jAij and d out i ¼ P jAji are the in and out degrees of node i, respectively, and h�i is an average over nodes. For the network topology studied by Ref. [2], � ¼ 1marks the critical regime in which the expected number of excited nodes is equal in consecutive time steps. While � ¼ 1 successfully predicts the critical regime for Erdo˝s-Re´nyi random networks [2], this prediction fails in networks with a more heterogeneous degree distribution [6,7]. Perhaps more importantly, previous theoretical analyses [2,6,7] do not account for features that are commonly found in real networks, such as community structure, correlation be- tween in and out degree of a given node, or correlation between the degree of two nodes at the ends of a given edge [12]. Here, we will generalize the mean-field criterion � ¼ 1 to account for complex network topologies. To begin, we note that lim�!0F ¼ 0 corresponds to the fixed point ~p ¼ 0 of Eq. (1) with � ¼ 0. To examine the linear stability of this fixed point, we set � ¼ 0 and lin- earize around pti ¼ 0, assuming pti to be small, obtaining ptþ1i ¼ PN j p t jAij. Assuming p t i ¼ ui�t yields �ui ¼ XN j ujAij: (2) Thus, the stability of the solution ~p ¼ 0 is governed by the largest eigenvalue of the network adjacency matrix �, with � < 1 being stable and � > 1 being unstable. Therefore, the critical state described in previous literature, occurring at various values of hdi, should universally occur at � ¼ 1. Importantly, since Aij � 0, the Perron-Frobenius theorem guarantees that � is real and positive [13]. Other previous studies in random networks have also investigated spectral properties of A to gain insight on the stability of dynamics in neural networks [14] and have shown how � could be changed by modifying the distribution of synapse strengths [15]. An important implication of Eq. (2) is that, when p and � are small enough, p should be almost proportional to the right eigenvector u corresponding to �, so we write pi ¼ Cui þ �i, where C is a proportionality constant and �i is a small error term. To first order, the constant C is related to F since, neglecting �, F ¼ hfit ¼ 1N X i pi � 1N X i Cui ¼ Chui: (3) The linear analysis allowed us to identify � ¼ 1 as the point at which the network response becomes nonzero as �! 0. In what follows, we use a weakly nonlinear analysis to obtain approximations to the response Fð�Þ when � is small. As we will show, these approximations depend only on a few spectral properties of A. Assuming Aijpj � 1 (which is valid near the critical regime if each node has many incoming connections), we approximate the product term of Eq. (1) with an exponential, obtaining in steady state pi ¼ ð1�piÞ � �þ ð1��Þ � 1� exp � �X j pjAij ��� ; (4) which we expand to second order using pi ¼ Cui þ �i and Au ¼ �u, Cui þ �i ¼ ðA�Þi þ �ð1� CuiÞ þ ð1� �Þ�Cui � � �þ 1 2 �2 � C2u2i : (5) To eliminate the error term �i from Eq. (5), we multiply by vi, the ith entry of the left eigenvector corresponding to �, and sum over i. We use the fact that vTA� ¼ �vT�, where vT denotes the transpose of v, and neglect the resulting small term ð1� �ÞPivi�i close to � ¼ 1, obtaining Chuvi ¼ �ðhvi � ChuviÞ þ ð1� �ÞC�huvi � � �þ 1 2 �2 � C2hvu2i: (6) This equation is quadratic in C [and therefore in F, via Eq. (3)] and linear in �, and may be easily solved for either. For � ¼ 0 the nonzero solution for F is F�¼0 ¼ ð�� 1Þ½�þ ½ð1=2Þ�2� huvihui hu2vi : (7) A crude nonperturbative approximation can be obtained by repeating this process without expanding Eq. (4), which yields the linear equation for � Chuvi¼X i ð1�CuiÞf�þð1��Þ½1�expð��CuiÞ�g: (8) Before numerically testing our theory, we will explain how it relates to previous results. For a network with correlations between degrees at the ends of a randomly chosen edge (assortative mixing by degree [12]), measured by the correlation coefficient � ¼ hdini doutj ie=hdindouti [16], with h�ie denoting an average over edges, the largest eigenvalue may be approximated by � � �hdindouti=hdi [16]. In the absence of assortativity, when � ¼ 1, � � hdindouti=hdi. If, in addition, there are no correlations between din and dout (node degree correlations) or if the degree distribution is sufficiently homogenous, then hdindouti � hdi2 and the approximation reduces to � � hdi. In the case of Ref. [2], � � hdi applies, and in the case of Refs. [6,7], � � hdindouti=hdi applies. We test our theoretical results via direct simulation of the Kinouchi-Copelli model on six categories of directed networks with N ¼ 10 000 nodes: (category 1) Random networks with no node degree correlation between din and dout; (category 2) random networks with maximal degree correlation, din ¼ dout; (category 3) random networks with moderate correlation between din and dout; (category 4) networks with power-law degree distribution PRL 106, 058101 (2011) P HY S I CA L R EV I EW LE T T E R S week ending 4 FEBRUARY 2011 058101-2 with power-law exponents � 2 ½2:0; 6:0�, with and without node degree correlations; (category 5) networks con- structed with hdi ¼ 1, and assortativity coefficient � varying in ½0:7; 1:3�; (category 6) networks with weights which depend on the degree of the node from which the edge originates, Aij ¼ �=douti . We created networks in multiple steps: first, we created binary networks (Aij 2 f0; 1g) with target degree distribu- tions as described below; next, we assigned a weight to each link, drawn from a uniform distribution between 0 and 1; finally, we calculated � for the resulting network and multiplied A by a constant to rescale the largest eigen- value to the targeted eigenvalue. The initial binary net- works in categories 1–3 were Erdo˝s-Re´nyi random networks, with hdi ¼ 10 [17]. Maximal degree correlation resulted from creating undirected binary networks and then forcing Aij ¼ Aji for i < j. Moderate degree correlation resulted from making undirected binary networks but al- lowing Aij � Aji when weights were assigned. The initial binary networks of categories 4–6 were constructed using the configuration model [18]. We enforced a minimum degree of 10 and a maximum degree of 200. In creating category 5 networks, we initially created one scale-free network with power-law exponent � ¼ 2:5 and � ¼ 1. Then we modified this original network by choosing two links at random and swapping them if the resulting swap would change the assortativity in the direction desired. This process was repeated until a desired value of � was achieved. Using this method we maintain exactly the same degree distribution and mean degree, yet modify � by virtue of � / �. In the six network types tested, results of simulations unanimously confirm the hypothesis that criticality occurs only for largest eigenvalue � ¼ 1. We present representa- tive results in Fig. 1(a), noting that each line and set of points corresponds to a single network realization, imply- ing that the effect of the largest eigenvalue on criticality is robust for individual systems. Figure 1(a) shows the response F as a function of stimulus � for scale-free networks with exponent � ¼ 2:5, constructed with no correlation between in and out degree, highlighting the significant difference between the regimes of � < 1 and � > 1, with the critical data corresponding to � ¼ 1. The lines were obtained by using Eqs. (3) and (8). Figure 1(b) shows � as a function of �, using �high ¼ 1 and �low ¼ 0:01, with the maximum occurring at � ¼ 1. Similar re- sults showing criticality and maximum dynamic range at � ¼ 1 are obtained for networks of all categories 1–5. Figure 2 shows F�!0 for networks of categories 3–5, confirming the transition predicted by the leading order analysis in Eq. (2). The symbols show the result of direct numerical simulation of the Kinouchi-Copelli model, the solid lines were obtained by iterating Eq. (1), and the dashed lines were obtained from Eq. (7). Figure 2(a) shows that criticality occurs at � ¼ 1 (indicated by a vertical arrow) rather than at hdi ¼ 1 for a category 3 random network. Figure 2(b) shows that criticality occurs at � ¼ 1 for scale-free networks (category 4). Correlations between din and dout affect the point at which � ¼ 1 occurs (vertical arrows). In Fig. 2(c), the mean degree was fixed at hdi ¼ 1, while �was changed by modifying the assortative coefficient �. As predicted by the theory, there is a tran- sition at � ¼ 1 even though the mean degree is fixed. We now explore the question of what network topology will best enhance dynamic range. We consider the follow- ing approximate measure of dynamic range �, obtained by setting �high to one in the definition of �, � ¼ 10log101=��, where �� is the stimulus value corre- sponding to a lower threshold response F�. Since dynamic range is maximized at criticality, we set � ¼ 1, solve Eq. (6) for ��, substitute it into the definition of � using Eq. (3), retaining the leading order behavior to get �max ¼ 10log10 2 3F2� � 10log10 hvu 2i hvihui2 : (9) Since the entries of the right (left) dominant eigenvector are a first order approximation to the in degree (out degree) of the corresponding nodes [19], the second term suggests that maximum dynamic range should increase (decrease) as the degree distribution becomes more homogenous (het- erogeneous). For example, consider the case of an undir- ected, uncorrelated network, in which vi ¼ ui � di. The second term is then approximately �10log10ðhd3i=hdi3Þ, which is maximized when di is independent of i. This 0 0.5 1 1.5 2 20 25 30 35 40 dy na m ic ra ng e, ∆ largest eigenvalue, λ 10−5 10−4 10−3 10−2 10−1 100 10 −5 10−4 10−3 10−2 10−1 10 0 λ > 1 λ = 1 λ < 1 stimulus, η re sp on se , F (a) (b) FIG. 1 (color online). (a) Response F vs stimulus � for power- law networks with exponent � ¼ 2:5 and no correlation between din and dout. Equation (8) (lines) captures the behavior of the simulation (circles), particularly for low levels of � and F. (b) Dynamic range � is maximized at � ¼ 1 in both simulation results (circles) and Eq. (6) (line). PRL 106, 058101 (2011) P HY S I CA L R EV I EW LE T T E R S week ending 4 FEBRUARY 2011 058101-3 corroborates the numerical findings in Refs. [2,7] that random graphs enhance dynamic range more than more heterogeneous scale-free graphs, and that the heterogeneity of the degree distribution affects dynamic range [7]. To test our result, we simulate scale-free networks with different power-law exponents � 2 ½2:0; 6:0�, yet with � ¼ 1 to maximize dynamic range in each case. Results of simula- tion (circles) plotted against the prediction of Eq. (9) (line) are shown in Fig. 3. In summary, we analytically predict and numerically confirm that criticality and peak dynamic range occur in networks with largest eigenvalue � ¼ 1. This result holds for diverse network topologies including random, scale- free, assortative, and/or degree-correlated networks, and for networks in which edge weights are related to nodal degree, thus generalizing previous work. Moreover, we find that homogeneous (heterogeneous) network topolo- gies result in higher (lower) dynamic range. Previous demonstrations of how � governs network dynamics in many other models (see [19] and references therein) sug- gest that the generality of our findings may extend beyond the particular model studied here. Taken together with related experimental findings [4], our results are consistent with the hypotheses that (1) real brain networks operate with � � 1, and (2) if an organism benefits from large dynamic range, then evolutionary pressures may act to homogenize the network topology of the brain. We thank Ed Ott and Dietmar Plenz for useful discus- sions. The work of W. L. S. was supported by the Intramural Research Program of the National Institute of Mental Health. J. G. R. was supported by NSF Grant No. DMS-0908221. *larremor@colorado.edu [1] L. L. Gollo et al., PLoS Comput. Biol. 5, e1000402 (2009). [2] O. Kinouchi et al., Nature Phys. 2, 348 (2006). [3] S. Gomez et al., Europhys. Lett. 89, 38 009 (2010). [4] W. L. Shew et al., J. Neurosci. 29, 15 595 (2009). [5] J.M. Beggs et al., J. Neurosci. 23, 11 167 (2003). [6] M. Copelli et al., Eur. Phys. J. B 56, 273 (2007). [7] A. Wu, X. J. Xu, and Y.H. Wang, Phys. Rev. E 75, 032901 (2007). [8] Our approach is easily generalized to include more re- fractory states. We also note that, in analogy to Ref. [9], our method can be generalized to include transmission delays and asynchronous updating. This will be discussed in a forthcoming publication. [9] E. Ott and A. Pomerance, Phys. Rev. E 79, 056111 (2009). [10] J. G. Restrepo, E. Ott, and B. R. Hunt, Phys. Rev. Lett. 100, 058701 (2008). [11] A. Pomerance et al., Proc. Natl. Acad. Sci. U.S.A. 106, 8209 (2009). [12] M. E. J. Newman, Phys. Rev. E 67, 026126 (2003). [13] C. R. MacCluer, SIAM Rev. 42, 487 (2000). [14] R. T. Gray et al., Neurocomputing 70, 1000 (2007). [15] K. Rajan and L. F. Abbott, Phys. Rev. Lett. 97, 188104 (2006). [16] J. G. Restrepo, E. Ott, and B. R. Hunt, Phys. Rev. E 76, 056119 (2007). [17] P. Erdo˝s et al., Publ. Math. 6, 290 (1959). [18] M. E. J. Newman, SIAM Rev. 45, 167 (2003). [19] J. G. Restrepo, E. Ott, and B. R. Hunt, Phys. Rev. Lett. 97, 094102 (2006). 2 3 4 5 6 30 32 34 36 38 power law exponent, γ dy na m ic ra ng e, Λ m AX FIG. 3. For power-law degree distributions with � ¼ 1, peak dynamic range increases monotonically with network homoge- neity, as measured by power-law exponent �. Simulations (circles) agree well with our predictions [Eq. (9); line]. 0 0.1 0.2 0.3 0 0.05 0.15 0.25 degree uncorrelated degree correlated 0 0.01 0.03 0.05 mean degree, d largest eigenvalue, λ re sp on se , F η 0 re sp on se , F η 0 re sp on se , F η 0 0 0.5 1 1.5 0 0.5 1 1.5 0.8 1 1.2 mean degree, d (a) random, degree correlated (b) scale free (c) scale free, assortative FIG. 2 (color online). F�!0 obtained from direct numerical simulation of the Kinouchi-Copelli model (symbols) plotted against hdi (a),(b) and � (c). Blue solid lines result from iterating Eq. (1) and green dashed lines result from Eq. (7). Small arrows