Backtracking activation impacts the criticality of excitable networks
aa r X i v : . [ n li n . AO ] D ec Backtracking activation impacts the criticality of excitablenetworks
Renquan Zhang and Guoyi Quan
School of Mathematical Sciences, Dalian University of Technology, Dalian 116024, China
Jiannan Wang ∗ Research Institute of Frontier Science,Beihang University, Beijing 100191, ChinaKey Laboratory of Mathematics InformaticsBehavioral Semantics, Ministry of Education, China
Sen Pei † Department of Environmental Health Sciences, Mailman School of Public Health,Columbia University, New York, NY 10032, USA bstract Networks of excitable elements are widely used to model real-world biological and social systems.The dynamic range of an excitable network quantifies the range of stimulus intensities that canbe robustly distinguished by the network response, and is maximized at the critical state. In thisstudy, we examine the impacts of backtracking activation on system criticality in excitable networksconsisting of both excitatory and inhibitory units. We find that, for dynamics with refractory statesthat prohibit backtracking activation, the critical state occurs when the largest eigenvalue of theweighted non-backtracking (WNB) matrix for excitatory units, λ ENB , is close to one, regardless ofthe strength of inhibition. In contrast, for dynamics without refractory state in which backtrackingactivation is allowed, the strength of inhibition affects the critical condition through suppressionof backtracking activation. As inhibitory strength increases, backtracking activation is graduallysuppressed. Accordingly, the system shifts continuously along a continuum between two extremeregimes – from one where the criticality is determined by the largest eigenvalue of the weightedadjacency matrix for excitatory units, λ EW , to the other where the critical state is reached when λ ENB is close to one. For systems in between, we find that λ ENB < λ EW > Keywords: excitable networks, criticality, backtracking activation, non-backtracking matrix, excitatory-inhibitory networks, dynamic range ∗ [email protected] † [email protected] . INTRODUCTION Excitable networks have been used to model a range of phenomena in biological andsocial systems including signal propagation in neural networks [1–6], information processingin brain networks [7–9], epidemic spread in human population and information diffusion insocial networks [10–13]. The collective dynamics of excitable nodes enable the networkedsystem to distinguish stimulus intensities varied by several orders of magnitude, character-ized by a large dynamic range in response to external stimuli. In previous studies, it wasfound that, for a number of excitable network models, the dynamic range is maximized atthe critical state [1, 14–17]. As a result, understanding the condition of criticality is essentialfor improving the performance and functionality of excitable networks.The critical condition for excitable networks composed of only excitatory nodes has beenextensively studied. In homogeneous random networks, the critical state corresponds to aunit branching ratio [1]. For more general network structures, the criticality for dynamicswithout refractory state is characterized by the unit largest eigenvalue of the weighted adja-cency matrix [14]. More recently, it was shown that for dynamics with refractory states, thecritical state is governed by the largest eigenvalue of the weighted non-backtracking (WNB)matrix [17]. In these studies, the largest eigenvalue of the weighted adjacency matrix orWNB matrix is used to define the critical state of excitable networks. However, when in-hibitory nodes are introduced, it is unclear how criticality is related to the largest eigenvaluesof these two matrices.In this study, we explore the critical condition of excitable networks consisting of bothexcitatory (E) and inhibitory (I) nodes. Inhibition presents in many real-world systems andplays a critical role in model dynamics and functions [18–23]. For instance, the introductionof inhibitory nodes into an excitable network operating near the critical state leads to self-sustained network activity [22], and inhibitory connectivity may be essential in maintaininglong-term information storage in volatile cortex [23]. In order to elucidate the relationshipbetween criticality and the largest eigenvalues of the weighted adjacency matrix and WNBmatrix, we study an excitatory-inhibitory (EI) network model equipped with a threshold-likeactivation rule [22]. Specifically, we focus on the impact of backtracking activation paths onthe critical condition.We first analyze the model dynamics of EI networks in two extreme conditions where3acktracking activation is allowed without restrictions or entirely prohibited. We find that,in the former case, the critical state is better characterized by the largest eigenvalue of theweighted adjacency matrix for excitatory nodes, λ EW , while in the latter case, the criticalityis more related to the largest eigenvalue of the WNB matrix for excitatory nodes, λ ENB .For EI models with refractory states that preclude backtracking activation, the criticalstate is achieved when λ ENB is close to one. For EI models without refractory state (i.e.,with only resting and excited states), however, the analytical form of the critical conditionbecomes intractable. We show that, qualitatively, the system gradually shifts from theformer case to the latter case as the strength of inhibition increases: for negligible inhibition, λ EW is closer to one at the critical state; for strong inhibition, λ ENB is closer to one atthe critical state; for moderate inhibition, we find λ ENB < λ EW >
2. MATERIALS AND METHODSA. The excitable network model
We consider excitable networks consisting of both excitatory (E) and inhibitory (I)nodes [22]. Contrary to the function of excitatory nodes, the effect of inhibitory nodesis to decrease the activation probability of their neighbors once they are activated. In anetwork with N nodes, we use s i ( t ) to represent the state of node i at time t . Both types ofnodes can be in one of m + 1 states: the resting state s i ( t ) = 0, the excited state s i ( t ) = 1,and refractory states s i ( t ) = 2 , , · · · , m . At each discrete time t , a resting node can beactivated by an external stimulus with a probability η , or activation propagation from itsneighbors independently. Specifically, the signal input strength from a node j to a neigh-boring node i , denoted by a ij , satisfies a ij > j is excitatory and a ij < j is inhibitory. If node j and node i are not connected in the network, we set a ij = 0. Theweighted adjacency matrix A = { a ij } N × N thus fully describes the network structure as well4s the signal input strength between all pairs of nodes. If a resting node i is not activatedby the external stimulus, its activation probability in the next time step t + 1 is calculatedby summing inputs from all neighbors through a transfer function σ ( · ): s i ( t + 1) = 1 with probability σ N X j =1 a ij τ ( s j ( t )) ! . (1)Here σ ( · ) is a piecewise linear function: σ ( x ) = 0 for x ≤ σ ( x ) = x for 0 < x <
1; and σ ( x ) = 1 for x ≥ τ ( · ) is a characteristic function: when s j ( t ) = 1, τ ( s j ( t )) = 1; otherwise, τ ( s j ( t )) = 0. A node can be activated by its neighbors if the net input is positive. Onceactivated, node i will transit to refractory states deterministically. That is, s i ( t +1) = s i ( t )+1if 1 ≤ s i ( t ) < m and s i ( t +1) = 0 if s i ( t ) = m . Note that, if m = 1, there will be no refractorystate and the node will directly return to the resting state after activation.In this study, we use undirected networks in which signals can be transmitted in bothdirections, and assume that the number of E nodes N e is larger than the number of Inodes N i . Empirical studies on real-world excitable systems reveals that the majority ofneural inputs are excitatory [24]. For ease of analysis, we rearrange node indices so thatnodes with index 1 ≤ i ≤ N e are excitatory and the rest are inhibitory. We consider bothhomogeneous and heterogeneous networks. For homogeneous network structure, we firstgenerate two Erd˝os-R´enyi (ER) random networks consisting of N e excitatory nodes and N i inhibitory nodes. Within each network, each pair of nodes is connected with a probability α . We then randomly connect E nodes and I nodes with a probability β . For heterogeneousnetwork structure, two scale-free (SF) networks of E nodes and I nodes with a power-lawdegree distribution P ( k ) ∝ k − γ are generated using the configuration model [25]. These twonetworks are then connected by randomly linking E nodes and I nodes with a probability β . In other words, networks of same types of nodes are constructed using either an ERmodel with a pairwise linking probability α or a configuration model for SF networks witha power-law exponent γ . Different types of nodes are connected randomly using a pairwiselinking probability β . We assume the absolute values of link weights | a ij | are distributeduniformly within a range, and the effect of inhibitory nodes is solely represented by theconnections from I to E nodes.In real-world excitable systems, the majority of units are excitatory. For instance, incortex, approximately 20% of neuron inputs are inhibitory [24]. Together with the fact thatthe inhibitory nodes need to be activated first before their release of inhibitory signals, the5ffect of I-I interactions is quite nominal in response to weak external stimuli. In somemodeling studies, the I-I interactions were even neglected [26]. In addition, we focus on therelative strength of same-type and cross-type interactions. It would be sufficient to keep oneconstant and vary the other (i.e., β ). Considering these factors, we decided to study theeffect of a varying β . B. Dynamic range and criticality
The dynamic range of an excitable network measures the range of stimulus intensities thatare distinguishable based on network response [1]. For a given stimulus intensity η ∈ (0 , F is defined as F = lim T →∞ T T X t =0 S t , (2)where S t is the fraction of excited nodes at time t . The response F ( η ) increases monotonicallywith a growing intensity of external stimulus η in a nonlinear fashion (see figure 1). For astrong stimulus intensity η →
1, the response F ( η ) will saturate and retain at a maximumvalue F max = 1 / ( m + 1). For a negligible stimulus intensity η →
0, the minimum response F depends on the state of the excitable system. In subcritical state, F ( η ) is a linearfunction of η for η →
0, i.e., F ( η ) ∝ η , with F →
0. At the critical state, F stillapproaches to zero but the function F ( η ) becomes nonlinear: F ( η ) ∝ η / (see figure 1(a)).The exponent 1/2 is called the Steven’s exponent, which characterizes the criticality of thecollective dynamics [1, 27]. In supercritical state, the excitation caused by external stimuluscan be self-sustained. Therefore, F becomes a positive number.The dynamic range of an excitable network is defined based on the function F ( η ). Inparticular, we define dynamic range ∆ as∆ = 10 log η high η low , (3)where η high and η low are stimulus intensities corresponding to network responses F high and F low (here F x = F + x ( F max − F ) for x ∈ [0 , η . and η . ,discarding stimuli that are too close to saturation or too weak to be distinguished from F . Previous studies have demonstrated that dynamic range is maximized at the criticalstate of an excitable system [1, 14]. Without forcing of external stimuli, excitation activity6 -4 -2 η -4 -3 -2 -1 F -4 -2 η F (b)(a) m=1/2 FIG. 1. Network response F in response to external stimulus intensity η . We show the responsecurves for an EI network constructed by connecting two ER networks ( N e = 3 , N i = 2 , α = 3 × − , β = 1 × − , m = 4). The absolute link weights | a ij | are drawn uniformly from 0.1to 0.2. We multiply link wights with different values to adjust the system to stay in the subcritical,critical and supercritical states. The network response F is shown in logarithmic scale in (a) andlinear scale in (b). At the critical state, we show that F ( η ) ∝ η / in (a). The dynamic range ofthe system is calculated based on the response curve. will eventually die out in subcritical state but become self-sustained above the critical point.This feature allows us to identify the critical point using the maximization of dynamic range.The criticality of the proposed model is closely related to percolation. As pointed out in aprevious study [28], activation or disease transmission in networks can be mapped to a bondpercolation process. In models with only excitatory units, an absorbing state of activationextinction (or a negligible giant component in percolation) exists for a low transmissionprobability. The criticality of the system corresponds to the transition point of the stabilityof this absorbing state, which is determined by the branching ratio [1], defined as the averagenumber of secondary activations induced by one excited node. In disease transmission, thisquantity is referred to as the reproductive number. If the branching ratio is below one, theabsorbing state is stable as all activations eventually die out; in contrast, if the branchingratio is larger than one, the absorbing state becomes unstable and any initial activation will7mplify and converge to a non-zero absorbing state. This stability change is equivalent tothe emergence of a non-zero giant component in percolation. The introduction of inhibitoryunits reduces the transmission probability of activation, thus delays the emergence of anon-zero absorbing state. In particular, inhibitory inputs can be viewed as a means ofregulation of the system, and how inhibition is imposed on excitatory units could affectwhen the stability transition occurs. The phase transition in systems with both excitatoryand inhibitory units therefore depends on the competition and interplay between excitatoryand inhibitory signals. Depending on whether an excited node can be activated immediatelyafter its last excitation, the effect of inhibitory units is different. This study aims to explorehow inhibitory units would change the condition of criticality.
3. RESULTS
The number of refractory states in model dynamics determines whether backtrackingactivation is permitted. Backtracking activation describes the following phenomenon ofdynamical resonance: an excitatory node i activated at time t increases the activationprobability of its excitatory neighbors at time t + 1, which in turn increases the activationprobability of node i at time t + 2. This behavior is only possible when there is no refractorystate (i.e., m = 1) so that excited nodes can directly return to the resting state at time t + 1. For dynamics with refractory states (i.e., m > t will enterrefractory states at time t + 1 thus cannot be activated again at time t + 2. Following thisdynamical rule, any backtracking activation is prohibited. A. Dynamics without backtracking activation
We first analyze the simpler case where backtracking is precluded by the existence ofrefractory states (i.e., m > i to j ( i → j ), wecreate a “cavity” at node j by “virtually” removing it from the network, and examine theprobability of node i being activated in the absence of node j , denoted by p ti → j at time t . Theprocedure of creating a virtual cavity at node j blocks the backtracking path i → j → i ,8nd therefore excludes the contribution via the consecutive activation i → j → i to theactivation probability of node i . This framework precisely depicts the model dynamics withrefractory states.For sparse networks without too many short loops, the probabilities p ti → j for neighboringnodes are mutual independent. Under this condition, the probability p ti → j for each node i can be recursively written as follows: p t +1 i → j = (1 − m − X l =0 p t − li → j ) η + (1 − η ) σ X k ∈ ∂i \ j a ik p tk → i . (4)Here ∂i \ j is the set of neighbors of node i excluding j . The probability that node i isexcited at time t + 1, denoted by p t +1 i , is calculated by putting node j back to the network: p t +1 i = (1 − m − X l =0 p t − li ) " η + (1 − η ) σ X k ∈ ∂i a ik p tk → i ! . (5)Note that, although Eqs (4)-(5) are derived for locally tree-like sparse networks, it has beenfound that results based on the sparseness assumption work well even for some networkswith dense clusters [32].
1. Analysis in the case of negligible inhibition
The piecewise transfer function σ ( · ) imposes a threshold-like activation rule that dependson the collective dynamics of all neighbors. Because the value of net input P k ∈ ∂i \ j a ik p tk → i is unknown, it becomes complicated to expand the right-hand-side of Eq (4) except forsome extreme cases. Here, we consider a special case where the cross-type interaction β isnegligible, i.e., β →
0. Under this extreme condition, excitatory and inhibitory nodes ineffect form two nearly separate communities. In particular, inhibitory nodes are unlikely tobe activated in response to weak stimuli as they almost only receive signals from inhibitorypeers. As a consequence, it is suffice to consider only excitatory nodes to compute networkresponse.In the steady state, denote the limiting probabilities as lim t →∞ p ti → j = p i → j and lim t →∞ p ti = p i . For excitatory nodes, Eqs (4)-(5) becomes p i → j = (1 − mp i → j ) η + (1 − η ) X k ∈ ∂ E i \ j a ik p k → i , (6)9 i = (1 − mp i ) " η + (1 − η ) X k ∈ ∂ E i a ik p k → i , (7)where 1 ≤ i ≤ N e , 1 ≤ j ≤ N e and ∂ E i is the set of excitatory neighbors of node i .To solve the self-consistent equations, we introduce two auxiliary variables: G i → j ( η, p → ) = η +(1 − η ) P k ∈ ∂ E i \ j a ik p k → i , G i ( η, p → ) = η +(1 − η ) P k ∈ ∂ E i a ik p k → i . We rearrange Eqs (6)-(7)and obtain p i → j = G i → j ( η, p → ) mG i → j ( η, p → ) + 1 (8)and p i = G i ( η, p → ) mG i ( η, p → ) + 1 . (9)For η = 0 (that is, without external forcing), Eq (8) has a trivial solution: p i → j = 0 forall links i → j . The stability of this solution determines the critical state of the system.If the solution is stable, the network activity triggered by a weak stimulus will eventuallydisappear; otherwise, the response will maintain at a nonzero level.The stability of the zero solution depends on the Jacobian matrix c M E = {M Ek → l,i → j } defined on all pairs of links k → l and i → j between E nodes. Specifically, we have M Ek → l,i → j = ∂p i → j ∂p k → l (cid:12)(cid:12)(cid:12)(cid:12) { η =0 ,p i → j =0 } = ∂G i → j ∂p k → l ( mG i → j + 1) − mG i → j ∂G i → j ∂p k → l ( mG i → j + 1) = ∂G i → j ∂p k → l (cid:12)(cid:12)(cid:12)(cid:12) { η =0 ,p i → j =0 } . (10)Here G i → j = 0 when η = 0 and p i → j = 0 for all i → j . According to the definition of G i → j ,the partial derivative of G i → j is given by ∂G i → j ∂p k → l (cid:12)(cid:12)(cid:12)(cid:12) { η =0 ,p i → j =0 } = (cid:26) a lk if l = i and j = k ,0 otherwise. (11)The elements of c M E are M Ek → l,i → j = a lk if l = i and j = k and 0 otherwise. Notethat, M Ek → l,i → j is non-zero only if the links k → l and i → j are consecutive ( l = i ) butnot backtracking ( j = k ). The weighted non-backtracking (WNB) matrix, or Hashimotomatrix [33], has recently found applications in several problems in network science [30, 34–40]. Because the stability of the zero solution is determined by the largest eigenvalue λ ENB of c M E , the system reaches the critical state if λ ENB = 1.10 . Numerical results for dynamics with inhibition
For the general case where inhibition cannot be neglected, it is challenging to derive theanalytical condition of criticality from Eq (4). As a result, we have to use numerical methodsto find the critical state. In particular, we are interested in how the largest eigenvalue ofthe WNB matrix for E nodes λ ENB changes with the inhibitory strength β at the criticalstate. We treat the increasing level of inhibition as a perturbation to the special case of β = 0, and examine to what extend the critical condition λ ENB = 1 will remain valid. Inorder to tune the system to the critical state, for a fixed inhibitory strength β , we randomlydraw absolute link weights | a ij | from a uniform distribution between 0.1 and 0.2, and thenmultiply | a ij | with a varying constant until the dynamic range of the system is maximized(i.e., the critical state is reached). The largest eigenvalue λ ENB of c M E is then computed usinga power method [41]. In figure 2, we show the relationship between dynamic range and thelargest eigenvalues of four matrices (i.e., the weighted adjacency matrix for all nodes, theweighted non-backtracking matrix for all nodes, the weighted adjacency matrix for E nodes,and the weighted non-backtracking matrix for E nodes). The curves present the largesteigenvalues of different matrices at the critical state where dynamic range is maximized.We first analyze homogeneous network structure. Without loss of generality, we assumethere are 3 refractory states ( m = 4). For ER networks with N e = 3 ,
000 excitatory nodesand N i = 2 ,
000 inhibitory nodes, we set the within-type connection probability α = 3 × − .An increasing level of inhibitory strength β = 1 × − , 5 × − , 1 × − , 2 × − and3 × − are tested. For each β , we slowly tune the link weights to drive the system to thecritical state, and record the largest eigenvalue of the WNB matrix for E nodes λ ENB . Weperform 300 realizations of this procedure, and report the distributions of λ ENB in figure 3.For comparison, we also computed the largest eigenvalue of the weighted adjacency matrixfor E nodes λ EW at criticality.Interestingly, even with non-negligible inhibition, λ ENB consistently distributes around oneat the critical state for an increasing inhibitory strength β . In contrast, λ EW distributes wellabove one. We note that the variation of λ ENB and λ EW is attributed to the finite networksize and numerical inaccuracy, as pointed out in a previous study on excitable networkswith only E nodes [17]. The numerical results in figure 3 indicate that the criticality of EInetworks with refractory states occurs when λ ENB is close to one, regardless of the strength11 .5 1 1.5 2 λ D y na m i c r ange m=4 λ W λ NB λ EW λ ENB λ D y na m i c r ange m=1 (a) (b) FIG. 2. The relationship between dynamic range and the largest eigenvalues of four matrices fordynamics with 3 refractory states (a) and without refractory state (b). Here, λ W , λ NB , λ EW and λ ENB are the largest eigenvalues of the weighted adjacency matrix for all nodes, the weighted non-backtracking matrix for all nodes, the weighted adjacency matrix for E nodes, and the weightednon-backtracking matrix for E nodes, respectively. We perform the experiment on an EI networkconstructed using two ER networks ( N e = 3 , N i = 2 , α = 3 × − , β = 1 × − ), andvary link weights to change the state of the system. For each setting of link weights, we calculatethe dynamic range and the corresponding largest eigenvalues. The setting that maximizes thedynamic range corresponds to the critical state. We use this numerical method to find the criticalstate of an EI network. At criticality, we find that λ ENB is close to one for m = 4 and λ EW is closeto one for m = 1. λ W ( λ NB ) is always smaller than λ EW ( λ ENB ) due to the existence of inhibitorynodes. of inhibition β . A closer inspection of figure 3 reveals that the average value of λ ENB isslighted larger than one and slowly increases with β . This slight shift of λ ENB indicates thatinhibition does impact the critical condition but its impact is very limited.According to model dynamics, the function of I nodes is passive – they need to be firstactivated before they can release inhibition signals. Without external stimuli, the conductionof inhibitory signals proceeds as follows: a set of excited E nodes activate an I node at time t ; at time t + 1, the excited I node exerts inhibitory signals to all its neighbors, among whichthe excited E nodes enter refractory state. With the presence of nodes in refractory state,we hypothesize that the inhibition effect is weakened. To demonstrate this, we plot F as12 .9 1 1.1 1.2 λ P r obab ili t y β =0.0001 λ WE λ NBE λ β =0.0005 λ β =0.001 λ P r obab ili t y β =0.002 λ β =0.003 λ F (a)(d) (e)(b) (f)(c) FIG. 3. Distributions of λ EW and λ ENB at the critical state of homogeneous EI networks withrefractory states ( m = 4) (a-e). Networks are constructed by connecting two ER networks ( N e =3 , N i = 2 , α = 3 × − ), and varying the cross-type link probability β . For differentsettings of β , λ ENB is consistently distributed near one. The relationship between F and λ EW and λ ENB is shown in (f) for β = 1 × − (up triangle), 5 × − (square), 1 × − (down triangle),2 × − (circle), and 3 × − (diamond). The transition point from F = 0 to F > β . a function of λ ENB and λ EW for increasing β in figure 3(f). Although the inhibitory strength β intensifies, the F curve does not change significantly, especially near the transition pointfrom F = 0 to F >
0. This result directly shows that, for dynamics with refractory states,the impact of inhibitory nodes is rather limited near the critical state.We performed the same analysis in SF networks with N e = 6 , N i = 4 ,
000 and thepower-law exponent γ = 3. Results in figure 4 show that, consistent with the results forER networks, λ ENB is close to one at the critical state. The effect of inhibitory nodes nearthe critical state is also nominal as the F curves are almost identical for different values ofinhibitory strength β . 13 =0.0001 λ P r obab ili t y λ WE λ NBE β =0.0005 λ β =0.001 λ β =0.002 λ P r obab ili t y β =0.003 λ λ F (a) (b) (c)(d) (f)(e) FIG. 4. Distributions of λ EW and λ ENB at the critical state of heterogeneous EI networks withrefractory states ( m = 4) (a-e). Networks are constructed by connecting two SF networks ( N e =6 , N i = 4 , γ = 3), and varying the cross-type link probability β . For different settingsof β , λ ENB is consistently distributed near one. The relationship between F and λ EW and λ ENB isshown in (f) for β = 1 × − (up triangle), 5 × − (square), 1 × − (down triangle), 2 × − (circle), and 3 × − (diamond). The transition point from F = 0 to F > β . B. Dynamics with backtracking activation
We now explore the more complicated dynamics in which backtracking activation is al-lowed. In this case, nodes have only two states – resting and excited. For each node i , denote p ti as the probability that node i is excited at time t . According to the model dynamics, theevolution of p ti is described by p t +1 i = (1 − p ti ) " η + (1 − η ) σ N X k =1 a ik p tk ! . (12)Backtracking activation is properly represented in Eq (12): if we expand p tk on the right-hand-side of Eq (12) in terms of the activation probability at time t − p t +1 i becomesexplicitly dependent on p t − i . This implies, the activation probability of each E node ata given time can contribute to the probability of its re-activation two time-steps later (as14ong as at least one of its E neighbors are activated), which exactly depicts the effect ofbacktracking activation.
1. Analysis in the case of negligible inhibition
Similar with our analysis of dynamics with refractory states, we first explore the extremecase where the cross-type linking probability β →
0. In this case, we only consider thenetwork of excitatory nodes. The stationary activation probability p i = lim t →∞ p ti satisfies p i = (1 − p i ) " η + (1 − η ) N e X k =1 a ik p k . (13)Without external stimuli, the system has a trivial solution p i = 0 for 1 ≤ i ≤ N e . Thestability of the zero solution is determined by the largest eigenvalue λ EW of the weightedadjacency matrix for excitatory nodes A E = { a ij } N e × N e . As a result, the critical state ischaracterized by λ EW = 1 as β →
2. Numerical results for dynamics with inhibition
We perturb the extreme case β → β , which introduces more inhibitory nodes connected to excitatory nodes. Withoutrefractory states, an “excitatory → inhibitory → excitatory” feedback loop appears: a groupof excited E nodes activate an inhibitory node; the excited I node then releases inhibitorysignals and decreases the activation probability of the E nodes who just activated it and nowreturned to the resting state. The inhibitory signals (negative inputs) impose a thresholdfor the re-activation of those E nodes. As a consequence, contributions from certain back-tracking paths may not be realized. This phenomenon is caused by the threshold-like featureof the transfer function σ ( · ). If the contribution of a backtracking path is lower than thethreshold imposed by inhibitory nodes, it may never contribute to the activation probabilityas σ ( x ) > x >
0. Following this line of reasoning, Eq (13) thus over-estimates the effect of backtracking activation when more inhibitory nodes are connected toexcitatory nodes. A stronger inhibitory strength β will suppress more backtracking activa-tions, which drives the dynamics of EI networks closer to the opposite extreme case wherebacktracking activation is entirely prohibited, described by Eqs (6)-(7).15e therefore hypothesize that, for a weak inhibitory strength β , λ EW is close to one atthe critical state; whereas for a strong inhibitory strength, λ ENB is close to one at the criticalstate. Varying the cross-type linking probability β modulates the system shifting betweenthese two extreme regimes. For an intermediate inhibitory strength β , we hypothesize that λ ENB < λ EW > λ EW and λ ENB atthe critical state for ER networks is shown in figure 5. In agreement with our hypothesis,as β increases, λ EW shifts from near one to above one, and λ ENB shifts from below one tonear one. The same phenomenon is also observed for SF networks in figure 6. In oder toexamine the effect of inhibitory nodes, we plot the F curve as a function of λ EW and λ ENB in figure 5(f) and figure 6(f). In contrast to dynamics with refractory states, introductionof more inhibitory links effectively reduces F , thus strongly impacts the critical state ofthe system. Such impact is reflected by the change of the transition point above which F becomes non-zero.Ideally, it would be desirable to show that the number of instances of backtracking ac-tivation decreases with an increasing inhibitory strength β . However, as the activation ofa node is collectively determined by a group of nodes, it is difficult to disentangle such in-teraction and identify definitively which backtracking path is responsible for the activation.Despite that, the impact of inhibitory nodes can be reflected by the threshold values thatthey impose on excitatory nodes. We calculate the average input from I nodes to E nodes inER and SF networks. Specifically, for a given stimulus intensity η , we compute the averageinput of resting E nodes from their excited inhibitory neighbors at each time step, and thenaverage over all time steps. In figure 7(a) and (c), we show that the average threshold indeedincreases monotonically with β . In addition, a stronger external stimulus η leads to a higheraverage threshold due to a larger number of excited nodes.We further calculate the fraction of excitatory links connected to resting E nodes whoseweights are lower than the threshold. To be specific, for each resting E node, we findits excited excitatory neighbors and focus on the links connected to them. These links arepotential candidates of backtracking activation, i.e., the actual backtracking activation pathsbelong to this set of links. Among these links, we calculate the proportion whose weights are16 .8 1 1.2 λ P r obab ili t y β =0.0001 λ WE λ NBE λ β =0.0005 λ β =0.001 λ P r obab ili t y β =0.002 λ β =0.003 λ F (b)(d) (e) (c)(f)(a) FIG. 5. Distributions of λ EW and λ ENB at the critical state of homogeneous EI networks withoutrefractory state ( m = 1) (a-e). Networks are constructed by connecting two ER networks ( N e =3 , N i = 2 , α = 3 × − ), and varying the cross-type link probability β . At the criticalstate, we show that in general λ EW > λ ENB <
1. As β increases, λ ENB becomes closer to oneand λ EW shifts away from one. The relationship between F and λ EW and λ ENB is shown in (f) for β = 1 × − (up triangle), 5 × − (square), 1 × − (down triangle), 2 × − (circle), and3 × − (diamond). The transition point from F = 0 to F > β . lower than the threshold of the E node. The contribution from such below-threshold linksare likely to be negated by the threshold. Therefore, the fraction of below-threshold linkscan partly reflect the magnitude of backtracking suppression. We present an illustrationfor computing this below-threshold fraction in figure 8. The mean fraction values averagedover all resting E nodes in all time steps for ER and SF networks are shown in figure 7(b)and (d). For both network structures, the fraction of below-threshold links increases as β grows with the presence of different levels of external stimuli. This analysis agrees with ourhypothesis and partially explains the transition between the two extreme cases.17 =0.0001 λ P r obab ili t y λ WE λ NBE β =0.0005 λ β =0.001 λ β =0.002 λ P r obab ili t y β =0.003 λ λ F (b)(a) (c)(d) (e) (f) FIG. 6. Distributions of λ EW and λ ENB at the critical state of heterogeneous EI networks withoutrefractory state ( m = 1) (a-e). Networks are constructed by connecting two SF networks ( N e =6 , N i = 4 , γ = 3), and varying the cross-type link probability β . At the critical state, weshow that in general λ EW > λ ENB <
1. As β increases, λ ENB becomes closer to one and λ EW shifts away from one. The relationship between F and λ EW and λ ENB is shown in (f) for β = 1 × − (up triangle), 5 × − (square), 1 × − (down triangle), 2 × − (circle), and 3 × − (diamond).The transition point from F = 0 to F > β . C. Simulations in synthetic neural networks
We further validate our findings in synthetic neural networks that have more realisticstructures. As it is difficult to find a real-world neural network dataset that contains bothexcitatory and inhibitory neurons, we have to construct a synthetic network using networkgeneration models [26, 42, 43]. In particular, networks of neurons in brain follow a clustered,distance-dependent connection pattern [43]. We generate networks with this organizationalpattern using a distance-dependent method employed in previous studies [43]. Specifically,3 ,
000 excitatory nodes and 2 ,
000 inhibitory nodes are placed uniformly on the surface of aunit sphere (figure 9(a)). The degree of each node was generated from a normal distribution.Nodes were then connected according to a distance-dependent probability P ∝ /d , where d is the geodesic distance between two nodes on the spherical surface. We assign the synaptic18 β × -3 T h r e s ho l d ER-ER η =0.01 η =0.1 η =0.2 η =0.5 β × -3 B e l o w - t h r e s ho l d f r a c t i on ER-ER β × -3 T h r e s ho l d SF-SF η =0.01 η =0.1 η =0.2 η =0.5 β × -3 B e l o w - t h r e s ho l d f r a c t i on SF-SF (a) (b)(c) (d)
FIG. 7. The threshold of resting E nodes imposed by their inhibitory neighbors for EI networksgenerated using ER ( N e = 3 , N i = 2 , α = 3 × − ) (a) and SF ( N e = 6 , N i = 4 , γ = 3) (c) networks. The fraction of below-threshold links for resting E nodes is reported in (b) and(d). For an increasing level of inhibition strength β , we tune the system to the critical state, andcalculate the threshold values and fraction of below-threshold links for different stimulus intensities η . As β and η increase, both threshold and below-threshold fraction increase. strength of all links from a uniform distribution between 0.1 and 0.2, and multiply a factor c to the weights of cross-type links in order to adjust the strength of inhibition. The valueof c reflects the inhibition strength in the system.We run simulations of model dynamics without refractory state ( m = 1) for c = 1. Wevary link weights, and calculate the dynamic range λ EW and λ ENB for each weight setting.The relationship between the dynamic range and λ EW and λ ENB is shown in figure 9(b).Similar with the results in random networks, at the critical state, we find that λ EW > λ ENB <
1. In the inset, we show the values of λ EW and λ ENB at the critical state for anincreasing inhibition strength c . As c grows, at the critical state, λ EW shifts away from oneand λ ENB gets closer to one. This result further corroborates our hypothesis that the system19
IG. 8. An example to show the calculation of threshold and fraction of below-threshold linksfor a resting E node i . Here, the resting E node has a total input | a ik + a ik | for its activated Ineighbors. This value is defined as the threshold. Among the 4 links connected to its activated Eneighbors, 2 links have weights below the threshold ( a ij < | a ik + a ik | , a ij < | a ik + a ik | ). Thefraction of below-threshold links is calculated as 2 / . lies between two extremes with ( λ EW ≈
1) and without ( λ ENB ≈
1) backtracking activation.
4. CONCLUSION
In this study, we explore the impact of backtracking activation on the criticality of ex-citable networks with both excitatory and inhibitory nodes. We find that, for dynamicswith refractory state that precludes backtracking activation, the critical state occurs whenthe largest eigenvalue of the WNB matrix for excitatory nodes is close to one. However, fordynamics without refractory state, the introduction of inhibitory nodes affects backtrackingactivation and the critical condition of the system. The EI model with inhibition essentiallyprovides an intermediate system between two extreme cases in which backtracking activationis allowed or prohibited. For the dynamics with a medium inhibitory strength, λ EW and λ ENB can be viewed as the upper and lower bound of the critical condition: at the critical state, λ EW > λ ENB <
1. In practice, this criterion can be used to assess whether a systemmay be at the critical state. If a system resides in a state where λ EW < λ ENB >
1, wecan assert that the system is not close to the critical state. Our results imply that a precisedescription of model dynamics is essential in theoretical analysis of phase transitions. These20 -11 0.5-0.5 0.5 00 0 -0.50.5 -0.5 -11 -1 λ D y na m i c R ange λ EW λ ENB c λ (b)(a) FIG. 9. (a) Generation of the synthetic neural network. 3,000 excitatory nodes (blue) and 2,000inhibitory nodes (red) are uniformly placed on the surface of a unit sphere. (b) Relationshipbetween the dynamic range and λ EW and λ ENB for a synthetic neural network ( N e = 3 , N i =2 , c = 1, the degree distribution follows a normal distribution with a mean of 9 and a varianceof 1). We vary link weights to change the state of the system, and calculate the dynamic range, λ EW and λ ENB for each setting. At the critical state where the dynamic range is maximized, we find λ ENB < λ EW >
1, which agrees with our hypothesis. Inset shows the values of λ EW and λ ENB at the critical state for an increasing inhibition strength c . findings highlight the important role of backtracking activation in excitable dynamics.Critical behavior is common in biological systems [44]. Besides the commonly addressedneuronal networks, fingerprint of criticality have been reported for calcium singallization inmyocytes [45], excitable beta cells [46], oocytes [47], etc. The operation of these biologicalsystems near critical states may be crucial for their proper functioning. Certain inhibitorymechanisms exist in cells to regulate the dynamics of calcium signaling, which could bepotentially relevant to our findings in this study. Further, several experimental studies onsubcellular [48] as well as on the cellular level [49] echo our findings that the refractoryperiods are important in excitable dynamics. Another relevant field is the study on pace-maker activities induced by intracellular calcium waves, which has been found essential inthe interstitial cell of Cajal (ICC) in the gastrointestinal tract and cardiac pacemaker cellsin the heart. It has been shown that refractory phases are crucial to prevent backtracking ofactivations in systems guided by pacemaker activities [50, 51]. Findings in this study may21nd applications in these biological systems in future works. It also merits further studywhether imposing inhibition on influencers of excitable dynamics would result in efficientregulation of model dynamics [52, 53]. DATA ACCESSIBILITY
No real-world data were used in this study.
AUTHORS’ CONTRIBUTIONS
R.Z. and S.P. designed the study. R.Z. and G.Q. wrote the code, ran simulations andperformed the analysis. S.P. wrote the first draft of the manuscript. R.Z., G.Q. and J.W.reviewed and edited the manuscript.
COMPETING INTERESTS
We declare we have no competing interests.
FUNDING
This work was supported by the National Natural Science Foundation of China (11801058),the Fundamental Research Funds for the Central Universities (DUT18RC(4)066) and BeijingNatural Science Foundation (1192012, Z180005). [1] Kinouchi O, Copelli M. 2006 Optimal dynamical range of excitable networks at criticality.
Nat. Phys. , 348-351. (doi:10.1038/nphys289)[2] Copelli M, Oliveira RF, Roque AC, Kinouchi O.2005 Signal compression in the sensory pe-riphery. Neurocomputing , 691-696.(doi:https://doi.org/10.1016/j.neucom.2004.10.099)[3] Gollo LL, Kinouchi O, Copelli M. 2009 Active dendrites enhance neuronal dynamic range.
PLoS Comput. Biol. , e1000402. (doi:10.1371/journal.pcbi.1000402)
4] Gollo LL, Copelli M, Roberts JA. 2016 Diversity improves performance in excitable networks.
PeerJ , e1912. (doi:10.7717/peerj.1912)[5] Wang CY, Wu ZX, Chen MZ. 2017 Approximate-master-equation approach forthe Kinouchi-Copelli neural model on networks. Phys. Rev. E , 012310 .(doi:10.1103/physRevE.95.012310)[6] Kinouchi O, Brochini L, Costa AA, Campos JGF, Copelli M. 2019 Stochastic oscilla-tions and dragon king avalanches in self-organized quasi-critical systems. Sci. Rep. , 3874.(doi:10.1038/s41598-019-40473-1)[7] Gautam SH, Hoang TT, McClanahan K, Grady SK, Shew WL. 2015 Maximizing sensorydynamic range by tuning the cortical state to criticality. PLoS Comput. Biol. , e1004576.(doi:10.1371/journal.pcbi.1004576)[8] Williams-Garc´ıa RV, Moore M, Beggs JM, Ortiz G. 2014 Quasicritical brain dynamics on anonequilibrium Widom line. Phys. Rev. E , 062714.(doi:10.1103/PhyRevE.90.062714)[9] Marro J, Mejias JF, Pinamonti G, Torres JJ. 2013 Signal transmission competing with noisein model excitable brains. AIP Conference Proceedings , 85. (doi:10.1063/1.4776504)[10] Karrer B, Newman ME. 2011 Competing epidemics on complex networks.
Phys. Rev. E ,036106. (doi:10.1103/PhysRevE.84.036106)[11] Van Mieghem P. 2012 Epidemic phase transition of the SIS type in networks. EPL , 48004.(doi:10.1209/0295-5075/97/48004)[12] Dodds PS, Harris KD, Danforth CM. 2013 Limited imitation contagion on randomnetworks: Chaos, universality, and unpredictability. Phys. Rev. Lett. , 158701.(doi:10.1103/PhysRevLett.110.158701)[13] Pei S, Tang S, Zheng Z. 2015 Detecting the influence of spreading in social networks withexcitable sensor networks.
PLoS ONE , e0124848. (doi:10.1371/journal.pone.0124848)[14] Larremore DB, Shew WL, Restrepo JG. 2011 Predicting criticality and dynamicrange in complex networks: effects of topology. Phys. Rev. Lett. , 058101.(doi:10.1103/PhysRevLett.106.058101)[15] Larremore DB, Shew WL, Ott E, Restrepo JG. 2011 Effects of network topology, transmissiondelays, and refractoriness on the response of coupled excitable systems to a stochastic stimulus.
Chaos , 025117. (doi:10.1063/1.3600760)
16] Copelli M, Kinouchi O. 2005 Intensity coding in two-dimensional excitable neural networks.
Physica A , 431-442.(doi:https://doi.org/ 10.1016/j.physa.2004.10.043)[17] Zhang R, Pei S. 2018 Dynamic range maximization in excitable networks.
Chaos , 013103.(doi:10.1063/1.4997254)[18] Adini Y, Sagi D, Tsodyks M. 1997 Excitatory-inhibitory network in the visual cortex: Psy-chophysical evidence. Proc. Natl. Acad. Sci. USA , 10426. (doi:10.1073/pnas.94.19.10426)[19] Park C, Terman D. 2010 Irregular behavior in an excitatory-inhibitory neuronal network. Chaos , 023122.(doi:10.1063/1.3430545)[20] Folias SE, Ermentrout GB. 2011 New patterns of activity in a pair of interacting excitatory-inhibitory neural fields. Phys. Rev. Lett. , 228103. (doi:10.1103/PhysRevLett.107.228103)[21] Pei S, Tang S, Yan S, Jiang S, Zhang X, Zheng Z. 2012 How to enhance the dy-namic range of excitatory-inhibitory excitable networks.
Phys. Rev. E , 021909.(doi:10.1103/PhysRevE.86.021909)[22] Larremore DB, Shew WL, Ott E, Sorrentino F, Restrepo JG. 2014 Inhibition causesceaseless dynamics in networks of excitable nodes. Phys. Rev. Lett. , 138103.(doi:10.1103/PhysRevLett.112.138103)[23] Mongillo G, Rumpel S, Loewenstein Y. 2018 Inhibitory connectivity defines the realm ofexcitatory plasticity.
Nat. Neurosci. , 1463-1470. (doi:10.1038/s41593-018-0226-x)[24] Soriano J, Mart´ınez MR, Tlusty T, Moses E. 2008 Development of input connections in neuralcultures. Proc. Natl. Acad. Sci. USA , 13758-13763. (doi:10.1073/pnas.0707492105)[25] Molloy M, Reed B. 1995 A critical point for random graphs with a given degree sequence.
Random Struct. Algor. , 161-180. (doi:10.1002/rsa.3240060204)[26] Stefanescu RA, Jirsa VK. 2008 A low dimensional description of globally coupled heteroge-neous neural networks of excitatory and inhibitory neurons. PLoS Comput. Biol. , e1000219.(doi:10.1371/journal.pcbi.1000219)[27] Copelli M, Roque AC, Oliveira RF, Kinouchi O. 2002 Physics of psychophysics: Stevensand Weber-Fechner laws are transfer functions of excitable media. Phys. Rev. E , 060901.(doi:10.1103/PhysRevE.65.060901)[28] Newman M. E. J. 2002 Spread of epidemic disease on networks. Phys. Rev. E , 016128.(doi:10.1103/PhysRevE.66.016128)
29] Mezard M, Montanari A. 2009
Information, Physics, and Computation . Oxford, UK: OxfordUniversity Press.[30] Morone F, Makse HA. 2015 Influence maximization in complex networks through optimalpercolation.
Nature
Sci. Rep. , 45240.(doi:10.1038/srep45240)[32] Melnik S, Hackett A, Porter MA, Mucha PJ, Gleeson JP. 2011 The unreasonable ef-fectiveness of tree-based theory for networks with clustering. Phys. Rev. E , 036112.(doi:10.1103/PhysRevE.83.036112)[33] Hashimoto K. 1989 Zeta functions of finite graphs and representations of p-adic groups. In Automorphic Forms and Geometry of Arithmetic Varieties (eds K Hashimoto, Y Namikawa) ,pp. 211-280. Academic Press.[34] Martin T, Zhang X, Newman ME. 2014 Localization and centrality in networks.
Phys. Rev.E , 052808. (doi:10.1103/PhysRevE.90.052808)[35] Karrer B, Newman ME, Zdeborov´a L. 2014 Percolation on sparse networks. Phys. Rev. Lett. , 208702. (doi:10.1103/PhysRevLett.113.208702)[36] Hamilton KE, Pryadko LP. 2014 Tight lower bound for percolation threshold on an infinitegraph.
Phys. Rev. Lett. , 208701. (doi:10.1103/PhysRevLett.113.208701)[37] Teng X, Pei S, Morone F, Makse HA. 2016 Collective influence of multiple spreaders eval-uated by tracing real information flow in large-scale social networks.
Sci. Rep. , 36043.(doi:0.1038/srep36043)[38] Wang J, Pei S, Wei W, Feng X, Zheng Z. 2018 Optimal stabilization of Boolean networksthrough collective influence. Phys. Rev. E , 032305. (doi:10.1103/PhysRevE.97.032305)[39] Wang J, Zhang R, Wei W, Pei S, Zheng Z. 2019 On the stability of multilayer Boolean networksunder targeted immunization. Chaos , 013133. (doi:10.1063/1.5053820)[40] Aleja D, Criado R, del Amo AJG, P´erez ´A, Romance M. 2019 Non-backtracking PageRank:From the classic model to hashimoto matrices. Chaos, Solitons & Fractals , 283-291.(doi:10.1016/j.chaos.2019.06.017)[41] Saad Y. 2011
Numerical Methods for Large Eigenvalue Problems.
Philadelphia, US: SIAM.
42] Markram H, Toledo-Rodriguez M, Wang Y, Gupta A, Silberberg G, Wu C. 2004 Interneuronsof the neocortical inhibitory system.
Nat. Rev. Neurosci. , 793. (doi:10.1038/nrn1519)[43] Wiles L, Gu S, Pasqualetti F, Parvesse B, Gabrieli D, Bassett DS, Meaney DF.2017 Autaptic connections shift network excitability and bursting. Sci. Rep. , 44006.(doi:10.1038/srep44006)[44] Mu˜noz MA. 2018 Colloquium: Criticality and dynamical scaling in living systems. Rev. Mod.Phys. , 031001. (doi:10.1103/RevModPhys.90.031001)[45] Nivala M, Ko CY, Nivala M, Weiss JN, Qu Z. 2012 Criticality in intracellular calcium signalingin cardiac myocytes. Biophys. J. , 2433-2442. (doi:10.1016/j.bpj.2012.05.001)[46] Stoˇzer A, Markoviˇc R, Dolenˇsek J, Perc M, Slak Rupnik M, Marhl M, Gosak M. 2019 Het-erogeneity and delayed activation as hallmarks of self-organization and criticality in excitabletissue.
Front. Psychol. , 869. (doi:10.3389/fphys.2019.00869)[47] Lopez L, Piegari E, Sigaut L, Ponce Dawson S. 2012 Intracellular calcium signalsdisplay an avalanche-like behavior over multiple lengthscales. Front. Psychol. , 350.(doi:10.3389/fphys.2012.00350)[48] Thul R, Coombes S, Roderick HL, Bootman MD. 2012 Subcellular calcium dynamics ina whole-cell model of an atrial myocyte. Proc. Natl. Acad. Sci. USA , 2150-2155.(doi:10.1073/pnas.1115855109)[49] Schuster S, Marhl M, H¨ofer T. 2002 Modelling of simple and complex calcium oscilla-tions: From singlecell responses to intercellular signalling.
Eur. J. Biochem. , 1333-1355.(doi:10.1046/j.0014-2956.2001.02720.x)[50] Gosak M, Marhl M, Perc M. 2009 Pacemaker-guided noise-induced spatial periodicity in ex-citable media.
Physica D , 506-515. (doi:10.1016/j.physd.2008.11.007)[51] Means SA, Sneyd J. 2010 Spatio-temporal calcium dynamics in pacemaking units of the in-terstitial cells of Cajal.
J. Theor. Biol. , 137-152. (10.1016/j.jtbi.2010.08.008)[52] Pei S, Wang J, Morone F, Makse HA. 2019 Influencer identification in dynamical complexsystems.
J. Complex Netw.
In Press.[53] Pei S, Morone F, Makse HA. 2018 Theories for influencer identification in complex net-works. In
Complex Spreading Phenomena in Social Systems pp. 125-148. Springer, Cham.(doi:10.1007/978-3-319-77332-2 8)pp. 125-148. Springer, Cham.(doi:10.1007/978-3-319-77332-2 8)