Backtracking activation impacts the criticality of excitable networks

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 can be robustly distinguished by the network response, and is maximized at the critical state. In this study, we examine the impacts of backtracking activation on system criticality in excitable networks consisting of both excitatory and inhibitory units. We find that, for dynamics with refractory states that prohibit backtracking activation, the critical state occurs when the largest eigenvalue of the weighted non-backtracking (WNB) matrix for excitatory units, $\lambda^E_{NB}$, is close to one, regardless of the strength of inhibition. In contrast, for dynamics without refractory state in which backtracking activation is allowed, the strength of inhibition affects the critical condition through suppression of backtracking activation. As inhibitory strength increases, backtracking activation is gradually suppressed. Accordingly, the system shifts continuously along a continuum between two extreme regimes -- from one where the criticality is determined by the largest eigenvalue of the weighted adjacency matrix for excitatory units, $\lambda^E_W$, to the other where the critical state is reached when $\lambda_{NB}^E$ is close to one. For systems in between, we find that $\lambda^E_{NB}<1$ and $\lambda^E_W>1$ at the critical state. These findings, confirmed by numerical simulations using both random and synthetic neural networks, indicate that backtracking activation impacts the criticality of excitable networks.


INTRODUCTION
Excitable networks have been used to model a range of phenomena in biological and social systems including signal propagation in neural networks [1][2][3][4][5][6], information processing in brain networks [7][8][9], epidemic spread in human population and information diffusion in social networks [10][11][12][13]. The collective dynamics of excitable nodes enable the networked system to distinguish stimulus intensities varied by several orders of magnitude, characterized by a large dynamic range in response to external stimuli. In previous studies, it was found that, for a number of excitable network models, the dynamic range is maximized at the critical state [1, [14][15][16][17]. As a result, understanding the condition of criticality is essential for improving the performance and functionality of excitable networks.
The critical condition for excitable networks composed of only excitatory nodes has been extensively studied. In homogeneous random networks, the critical state corresponds to a unit branching ratio [1]. For more general network structures, the criticality for dynamics without refractory state is characterized by the unit largest eigenvalue of the weighted adjacency matrix [14]. More recently, it was shown that for dynamics with refractory states, the critical 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 or WNB matrix is used to define the critical state of excitable networks. However, when inhibitory nodes are introduced, it is unclear how criticality is related to the largest eigenvalues of these two matrices.
In this study, we explore the critical condition of excitable networks consisting of both excitatory (E) and inhibitory (I) nodes. Inhibition presents in many real-world systems and plays a critical role in model dynamics and functions [18][19][20][21][22][23]. For instance, the introduction of inhibitory nodes into an excitable network operating near the critical state leads to selfsustained network activity [22], and inhibitory connectivity may be essential in maintaining long-term information storage in volatile cortex [23]. In order to elucidate the relationship between criticality and the largest eigenvalues of the weighted adjacency matrix and WNB matrix, we study an excitatory-inhibitory (EI) network model equipped with a threshold-like activation rule [22]. Specifically, we focus on the impact of backtracking activation paths on the critical condition.
We first analyze the model dynamics of EI networks in two extreme conditions where backtracking 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 the weighted adjacency matrix for excitatory nodes, λ E W , while in the latter case, the criticality is more related to the largest eigenvalue of the WNB matrix for excitatory nodes, λ E N B . For EI models with refractory states that preclude backtracking activation, the critical state is achieved when λ E N B 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 condition becomes intractable. We show that, qualitatively, the system gradually shifts from the former case to the latter case as the strength of inhibition increases: for negligible inhibition, λ E W is closer to one at the critical state; for strong inhibition, λ E N B is closer to one at the critical state; for moderate inhibition, we find λ E N B < 1 and λ E W > 1 at the critical state. Using numerical simulations in both random and synthetic networks, we verify that a larger inhibitory strength tends to suppress more backtracking activation, which explains the transition between these two regimes. Our findings highlight the impact of backtracking activation, a form of dynamical resonance, on the criticality of excitable networks, and may provide new insight into the study of similar dynamical processes in networked systems.

A. 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 nodes is to decrease the activation probability of their neighbors once they are activated. In a network with N nodes, we use s i (t) to represent the state of node i at time t. Both types of nodes 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, 3, · · · , m. At each discrete time t, a resting node can be activated by an external stimulus with a probability η, or activation propagation from its neighbors independently. Specifically, the signal input strength from a node j to a neighboring node i, denoted by a ij , satisfies a ij > 0 if node j is excitatory and a ij < 0 if node j is inhibitory. If node j and node i are not connected in the network, we set a ij = 0. The weighted adjacency matrix A = {a ij } N ×N thus fully describes the network structure as well as the signal input strength between all pairs of nodes. If a resting node i is not activated by the external stimulus, its activation probability in the next time step t + 1 is calculated by summing inputs from all neighbors through a transfer function σ(·): (1) Here σ(·) is a piecewise linear function: σ(x) = 0 for x ≤ 0; σ(x) = x for 0 < x < 1; and σ(x) = 1 for x ≥ 1. τ (·) 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. Once activated, node i will transit to refractory states deterministically. That is, Note that, if m = 1, there will be no refractory state 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 both directions, and assume that the number of E nodes N e is larger than the number of I nodes N i . Empirical studies on real-world excitable systems reveals that the majority of neural inputs are excitatory [24]. For ease of analysis, we rearrange node indices so that nodes with index 1 ≤ i ≤ N e are excitatory and the rest are inhibitory. We consider both homogeneous and heterogeneous networks. For homogeneous network structure, we first generate two Erdős-Rényi (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 heterogeneous network structure, two scale-free (SF) networks of E nodes and I nodes with a power-law degree distribution P (k) ∝ k −γ are generated using the configuration model [25]. These two networks 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 ER model with a pairwise linking probability α or a configuration model for SF networks with a power-law exponent γ. Different types of nodes are connected randomly using a pairwise linking probability β. We assume the absolute values of link weights |a ij | are distributed uniformly within a range, and the effect of inhibitory nodes is solely represented by the connections from I to E nodes.
In real-world excitable systems, the majority of units are excitatory. For instance, in cortex, approximately 20% of neuron inputs are inhibitory [24]. Together with the fact that the inhibitory nodes need to be activated first before their release of inhibitory signals, the effect of I-I interactions is quite nominal in response to weak external stimuli. In some modeling studies, the I-I interactions were even neglected [26]. In addition, we focus on the relative strength of same-type and cross-type interactions. It would be sufficient to keep one constant and vary the other (i.e., β). Considering these factors, we decided to study the effect of a varying β.

B. Dynamic range and criticality
The dynamic range of an excitable network measures the range of stimulus intensities that are distinguishable based on network response [1]. For a given stimulus intensity η ∈ (0, 1), the network response F is defined as where S t is the fraction of excited nodes at time t. The response F (η) increases monotonically with a growing intensity of external stimulus η in a nonlinear fashion (see figure 1). For a strong stimulus intensity η → 1, the response F (η) will saturate and retain at a maximum value F max = 1/(m + 1). For a negligible stimulus intensity η → 0, the minimum response F 0 depends on the state of the excitable system. In subcritical state, F (η) is a linear function of η for η → 0, i.e., F (η) ∝ η, with F 0 → 0. At the critical state, F 0 still approaches to zero but the function F (η) becomes nonlinear: F (η) ∝ η 1/2 (see figure 1(a)).
The exponent 1/2 is called the Steven's exponent, which characterizes the criticality of the collective dynamics [1,27]. In supercritical state, the excitation caused by external stimulus can be self-sustained. Therefore, F 0 becomes a positive number.
The dynamic range of an excitable network is defined based on the function F (η). In particular, we define dynamic range ∆ as ∆ = 10 log 10 η high η low , where η high and η low are stimulus intensities corresponding to network responses F high and In this study, we use η 0.9 and η 0.1 , discarding stimuli that are too close to saturation or too weak to be distinguished from  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 a previous study [28], activation or disease transmission in networks can be mapped to a bond percolation process. In models with only excitatory units, an absorbing state of activation extinction (or a negligible giant component in percolation) exists for a low transmission probability. The criticality of the system corresponds to the transition point of the stability of this absorbing state, which is determined by the branching ratio [1], defined as the average number of secondary activations induced by one excited node. In disease transmission, this quantity is referred to as the reproductive number. If the branching ratio is below one, the absorbing state is stable as all activations eventually die out; in contrast, if the branching ratio is larger than one, the absorbing state becomes unstable and any initial activation will amplify and converge to a non-zero absorbing state. This stability change is equivalent to the emergence of a non-zero giant component in percolation. The introduction of inhibitory units reduces the transmission probability of activation, thus delays the emergence of a non-zero absorbing state. In particular, inhibitory inputs can be viewed as a means of regulation of the system, and how inhibition is imposed on excitatory units could affect when the stability transition occurs. The phase transition in systems with both excitatory and inhibitory units therefore depends on the competition and interplay between excitatory and inhibitory signals. Depending on whether an excited node can be activated immediately after its last excitation, the effect of inhibitory units is different. This study aims to explore how inhibitory units would change the condition of criticality.

RESULTS
The number of refractory states in model dynamics determines whether backtracking activation is permitted. Backtracking activation describes the following phenomenon of dynamical resonance: an excitatory node i activated at time t increases the activation probability of its excitatory neighbors at time t + 1, which in turn increases the activation probability of node i at time t + 2. This behavior is only possible when there is no refractory state (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 > 1), nodes excited at time t will enter refractory states at time t + 1 thus cannot be activated again at time t + 2. Following this dynamical rule, any backtracking activation is prohibited.

A. Dynamics without backtracking activation
We first analyze the simpler case where backtracking is precluded by the existence of refractory states (i.e., m > 1). To account for the dynamics without backtracking, we formulate the model evolution in a message-passing framework [17], which is frequently used in statistical physics and network science [29][30][31]. For a link from i to j (i → j), we create a "cavity" at node j by "virtually" removing it from the network, and examine the probability of node i being activated in the absence of node j, denoted by p t i→j at time t. The procedure of creating a virtual cavity at node j blocks the backtracking path i → j → i, and therefore excludes the contribution via the consecutive activation i → j → i to the activation probability of node i. This framework precisely depicts the model dynamics with refractory states.
For sparse networks without too many short loops, the probabilities p t i→j for neighboring nodes are mutual independent. Under this condition, the probability p t i→j for each node i can be recursively written as follows: Here ∂i \ j is the set of neighbors of node i excluding j. The probability that node i is excited at time t + 1, denoted by p t+1 i , is calculated by putting node j back to the network: Note that, although Eqs (4)-(5) are derived for locally tree-like sparse networks, it has been found that results based on the sparseness assumption work well even for some networks with dense clusters [32].

Analysis in the case of negligible inhibition
The piecewise transfer function σ(·) imposes a threshold-like activation rule that depends on the collective dynamics of all neighbors. Because the value of net input k∈∂i\j a ik p t k→i is unknown, it becomes complicated to expand the right-hand-side of Eq (4) except for some extreme cases. Here, we consider a special case where the cross-type interaction β is negligible, i.e., β → 0. Under this extreme condition, excitatory and inhibitory nodes in effect form two nearly separate communities. In particular, inhibitory nodes are unlikely to be activated in response to weak stimuli as they almost only receive signals from inhibitory peers. As a consequence, it is suffice to consider only excitatory nodes to compute network response.
In the steady state, denote the limiting probabilities as lim t→∞ p t i→j = p i→j and lim t→∞ p t i = p i . For excitatory nodes, Eqs (4)-(5) becomes 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: We rearrange Eqs (6)- (7) and obtain and For η = 0 (that is, without external forcing), Eq (8) has a trivial solution: p i→j = 0 for all 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 eventually disappear; otherwise, the response will maintain at a nonzero level.
The stability of the zero solution depends on the Jacobian matrix M E = {M E k→l,i→j } defined on all pairs of links k → l and i → j between E nodes. Specifically, we have 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 The elements of M E are M E k→l,i→j = a lk if l = i and j = k and 0 otherwise. Note that, M E k→l,i→j is non-zero only if the links k → l and i → j are consecutive (l = i) but not backtracking (j = k). The weighted non-backtracking (WNB) matrix, or Hashimoto matrix [33], has recently found applications in several problems in network science [30,[34][35][36][37][38][39][40]. Because the stability of the zero solution is determined by the largest eigenvalue λ E N B of M E , the system reaches the critical state if λ E N B = 1.

Numerical results for dynamics with inhibition
For the general case where inhibition cannot be neglected, it is challenging to derive the analytical condition of criticality from Eq (4). As a result, we have to use numerical methods to find the critical state. In particular, we are interested in how the largest eigenvalue of the WNB matrix for E nodes λ E N B changes with the inhibitory strength β at the critical state. We treat the increasing level of inhibition as a perturbation to the special case of β = 0, and examine to what extend the critical condition λ E N B = 1 will remain valid. In order to tune the system to the critical state, for a fixed inhibitory strength β, we randomly draw absolute link weights |a ij | from a uniform distribution between 0.1 and 0.2, and then multiply |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 λ E N B of M E is then computed using a power method [41]. In figure 2, we show the relationship between dynamic range and the largest eigenvalues of four matrices (i.e., the weighted adjacency matrix for all nodes, the weighted 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 largest eigenvalues of different matrices at the critical state where dynamic range is maximized.
We first analyze homogeneous network structure. Without loss of generality, we assume there are 3 refractory states (m = 4). For ER networks with N e = 3, 000 excitatory nodes and N i = 2, 000 inhibitory nodes, we set the within-type connection probability α = 3×10 −3 .
An increasing level of inhibitory strength β = 1 × 10 −4 , 5 × 10 −4 , 1 × 10 −3 , 2 × 10 −3 and 3 × 10 −3 are tested. For each β, we slowly tune the link weights to drive the system to the critical state, and record the largest eigenvalue of the WNB matrix for E nodes λ E N B . We perform 300 realizations of this procedure, and report the distributions of λ E N B in figure 3. For comparison, we also computed the largest eigenvalue of the weighted adjacency matrix for E nodes λ E W at criticality.
Interestingly, even with non-negligible inhibition, λ E N B consistently distributes around one at the critical state for an increasing inhibitory strength β. In contrast, λ E W distributes well above one. We note that the variation of λ E N B and λ E W is attributed to the finite network size and numerical inaccuracy, as pointed out in a previous study on excitable networks with only E nodes [17]. The numerical results in figure 3 indicate that the criticality of EI networks with refractory states occurs when λ E N B is close to one, regardless of the strength Here, λ W , λ N B , λ E W and λ E N B are the largest eigenvalues of the weighted adjacency matrix for all nodes, the weighted nonbacktracking matrix for all nodes, the weighted adjacency matrix for E nodes, and the weighted non-backtracking matrix for E nodes, respectively. We perform the experiment on an EI network constructed using two ER networks (N e = 3, 000, N i = 2, 000, α = 3 × 10 −3 , β = 1 × 10 −3 ), and vary link weights to change the state of the system. For each setting of link weights, we calculate the dynamic range and the corresponding largest eigenvalues. The setting that maximizes the dynamic range corresponds to the critical state. We use this numerical method to find the critical state of an EI network. At criticality, we find that λ E N B is close to one for m = 4 and λ E W is close to one for m = 1. λ W (λ N B ) is always smaller than λ E W (λ E N B ) due to the existence of inhibitory nodes.
of inhibition β. A closer inspection of figure 3 reveals that the average value of λ E N B is slighted larger than one and slowly increases with β. This slight shift of λ E N B indicates that inhibition 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 first activated before they can release inhibition signals. Without external stimuli, the conduction of 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 which the 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 0 as

B. Dynamics with backtracking activation
We now explore the more complicated dynamics in which backtracking activation is allowed. In this case, nodes have only two states -resting and excited. For each node i, denote p t i as the probability that node i is excited at time t. According to the model dynamics, the evolution of p t i is described by Backtracking activation is properly represented in Eq (12): if we expand p t k on the righthand-side of Eq (12) in terms of the activation probability at time t − 1, p t+1 i becomes explicitly dependent on p t−1 i . This implies, the activation probability of each E node at a given time can contribute to the probability of its re-activation two time-steps later (as long as at least one of its E neighbors are activated), which exactly depicts the effect of backtracking activation.

Analysis in the case of negligible inhibition
Similar with our analysis of dynamics with refractory states, we first explore the extreme case where the cross-type linking probability β → 0. In this case, we only consider the network of excitatory nodes. The stationary activation probability p i = lim t→∞ p t i satisfies Without external stimuli, the system has a trivial solution p i = 0 for 1 ≤ i ≤ N e . The stability of the zero solution is determined by the largest eigenvalue λ E W of the weighted adjacency matrix for excitatory nodes A E = {a ij } Ne×Ne . As a result, the critical state is characterized by λ E W = 1 as β → 0.

Numerical results for dynamics with inhibition
We perturb the extreme case β → 0 by gradually increasing the cross-type linking probability β, which introduces more inhibitory nodes connected to excitatory nodes. Without refractory states, an "excitatory→inhibitory→excitatory" feedback loop appears: a group of excited E nodes activate an inhibitory node; the excited I node then releases inhibitory signals and decreases the activation probability of the E nodes who just activated it and now returned to the resting state. The inhibitory signals (negative inputs) impose a threshold for the re-activation of those E nodes. As a consequence, contributions from certain backtracking paths may not be realized. This phenomenon is caused by the threshold-like feature of the transfer function σ(·). If the contribution of a backtracking path is lower than the threshold imposed by inhibitory nodes, it may never contribute to the activation probability as σ(x) > 0 only if the net input x > 0. Following this line of reasoning, Eq (13) thus overestimates the effect of backtracking activation when more inhibitory nodes are connected to excitatory nodes. A stronger inhibitory strength β will suppress more backtracking activations, which drives the dynamics of EI networks closer to the opposite extreme case where backtracking activation is entirely prohibited, described by Eqs (6)- (7).
We therefore hypothesize that, for a weak inhibitory strength β, λ E W is close to one at the critical state; whereas for a strong inhibitory strength, λ E N B is close to one at the critical state. Varying the cross-type linking probability β modulates the system shifting between these two extreme regimes. For an intermediate inhibitory strength β, we hypothesize that λ E N B < 1 and λ E W > 1 at the critical state. We verify this hypothesis using numerical simulations in both homogeneous and heterogeneous networks.
We performed the same analysis as in figure 3 and figure 4, except using a different model dynamics with only resting and excited states. The distributions of λ E W and λ E N B at the critical state for ER networks is shown in figure 5. In agreement with our hypothesis, as β increases, λ E W shifts from near one to above one, and λ E N B shifts from below one to near one. The same phenomenon is also observed for SF networks in figure 6. In oder to examine the effect of inhibitory nodes, we plot the F 0 curve as a function of λ E W and λ E N B in figure 5(f) and figure 6(f). In contrast to dynamics with refractory states, introduction of more inhibitory links effectively reduces F 0 , thus strongly impacts the critical state of the system. Such impact is reflected by the change of the transition point above which F 0 becomes non-zero.
Ideally, it would be desirable to show that the number of instances of backtracking activation decreases with an increasing inhibitory strength β. However, as the activation of a node is collectively determined by a group of nodes, it is difficult to disentangle such interaction 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 that they impose on excitatory nodes. We calculate the average input from I nodes to E nodes in ER and SF networks. Specifically, for a given stimulus intensity η, we compute the average input of resting E nodes from their excited inhibitory neighbors at each time step, and then average over all time steps. In figure 7(a) and (c), we show that the average threshold indeed increases monotonically with β. In addition, a stronger external stimulus η leads to a higher average threshold due to a larger number of excited nodes.
We further calculate the fraction of excitatory links connected to resting E nodes whose weights are lower than the threshold. To be specific, for each resting E node, we find its excited excitatory neighbors and focus on the links connected to them. These links are potential candidates of backtracking activation, i.e., the actual backtracking activation paths belong to this set of links. Among these links, we calculate the proportion whose weights are  The transition point from F 0 = 0 to F 0 > 0 is significantly affected by the strength of inhibition β.

C. Simulations in synthetic neural networks
We further validate our findings in synthetic neural networks that have more realistic structures. As it is difficult to find a real-world neural network dataset that contains both excitatory and inhibitory neurons, we have to construct a synthetic network using network generation models [26,42,43]. In particular, networks of neurons in brain follow a clustered, distance-dependent connection pattern [43]. We generate networks with this organizational pattern 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 a unit 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 ∝ 1/d 2 , where d is the geodesic distance between two nodes on the spherical surface. We assign the synaptic We run simulations of model dynamics without refractory state (m = 1) for c = 1. We vary link weights, and calculate the dynamic range λ E W and λ E N B for each weight setting. The relationship between the dynamic range and λ E W and λ E N B is shown in figure 9(b). Similar with the results in random networks, at the critical state, we find that λ E W > 1 and λ E N B < 1. In the inset, we show the values of λ E W and λ E N B at the critical state for an increasing inhibition strength c. As c grows, at the critical state, λ E W shifts away from one and λ E N B gets closer to one. This result further corroborates our hypothesis that the system for a resting E node i. Here, the resting E node has a total input |a ik 2 + a ik 3 | for its activated I neighbors. This value is defined as the threshold. Among the 4 links connected to its activated E neighbors, 2 links have weights below the threshold (a ij 3 < |a ik 2 + a ik 3 |, a ij 4 < |a ik 2 + a ik 3 |). The fraction of below-threshold links is calculated as 2/4 = 0.5.
lies between two extremes with (λ E W ≈ 1) and without (λ E N B ≈ 1) backtracking activation.

CONCLUSION
In this study, we explore the impact of backtracking activation on the criticality of excitable networks with both excitatory and inhibitory nodes. We find that, for dynamics with refractory state that precludes backtracking activation, the critical state occurs when the largest eigenvalue of the WNB matrix for excitatory nodes is close to one. However, for dynamics without refractory state, the introduction of inhibitory nodes affects backtracking activation and the critical condition of the system. The EI model with inhibition essentially provides an intermediate system between two extreme cases in which backtracking activation is allowed or prohibited. For the dynamics with a medium inhibitory strength, λ E W and λ E N B can be viewed as the upper and lower bound of the critical condition: at the critical state, λ E W > 1 and λ E N B < 1. In practice, this criterion can be used to assess whether a system may be at the critical state. If a system resides in a state where λ E W < 1 or λ E N B > 1, we can assert that the system is not close to the critical state. Our results imply that a precise description of model dynamics is essential in theoretical analysis of phase transitions. These between the dynamic range and λ E W and λ E N B for a synthetic neural network (N e = 3, 000, N i = 2, 000, c = 1, the degree distribution follows a normal distribution with a mean of 9 and a variance of 1). We vary link weights to change the state of the system, and calculate the dynamic range, λ E W and λ E N B for each setting. At the critical state where the dynamic range is maximized, we find λ E N B < 1 and λ E W > 1, which agrees with our hypothesis. Inset shows the values of λ E W and λ E N B 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 addressed neuronal networks, fingerprint of criticality have been reported for calcium singallization in myocytes [45], excitable beta cells [46], oocytes [47], etc. The operation of these biological systems near critical states may be crucial for their proper functioning. Certain inhibitory mechanisms exist in cells to regulate the dynamics of calcium signaling, which could be potentially relevant to our findings in this study. Further, several experimental studies on subcellular [48] as well as on the cellular level [49] echo our findings that the refractory periods are important in excitable dynamics. Another relevant field is the study on pacemaker activities induced by intracellular calcium waves, which has been found essential in the interstitial cell of Cajal (ICC) in the gastrointestinal tract and cardiac pacemaker cells in the heart. It has been shown that refractory phases are crucial to prevent backtracking of activations in systems guided by pacemaker activities [50,51]. Findings in this study may find applications in these biological systems in future works. It also merits further study whether imposing inhibition on influencers of excitable dynamics would result in efficient regulation 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 and performed the analysis. S.P. wrote the first draft of the manuscript. R.Z., G.Q. and J.W.
reviewed and edited the manuscript.