PLoS Computational Biology
Public Library of Science
Scale free topology as an effective feedback system
Volume: 16 , Issue: 5
Doi: 10.1371/journal.pcbi.1007825

Table of Contents




Nature abounds with complex networks of interacting elements—from the proteins in our cells, through neural networks in our brains, to species interacting in ecosystems. In all of these fields, the relation between network structure and dynamics is an important research question. A recurring feature of natural networks is their heterogeneous structure: individual elements exhibit a huge diversity of connectivity patterns, which complicates the understanding of network dynamics. To address this problem, we devised a simplified approximation for complex structured networks which captures their dynamical properties. Separating out the largest “hubs”—a small number of nodes with disproportionately high connectivity—we represent them by a single node linked to the rest of the network. This enables us to borrow concepts from control theory, where a system’s output is linked back to itself forming a feedback loop. In this analogy, hubs in heterogeneous networks implement a feedback circuit with the rest of the network. The analogy reveals how these hubs can coordinate the network and drive it more easily towards stable states. Our approach enables analyzing dynamical properties of heterogeneous networks, which is difficult to achieve with existing techniques. It is potentially applicable to many fields where heterogeneous networks are important.

Rivkind, Schreier, Brenner, Barak, and Morozov: Scale free topology as an effective feedback system


Complex networks, their structure and dynamics, have become an indispensable tool for investigating biological systems [13]. Networks provide a mathematical model for a system of many interacting components; such models are a cornerstone of computational neuroscience and are widely used in cell biology to describe genetic networks, protein interactions, metabolism and more. Motivated at least partly by these applications, the mathematical theory of network analysis has gained general and fundamental interest and has advanced tremendously in the past decades. Topics of interest include network structure and topology; dynamic behaviour such as fixed points, their stability and their scaling with network size; robustness of these dynamic properties to noise and to evolution; and network controllability and information processing properties [47].

One important property of biological networks that has raised much interest is their heterogeneous topology. Analyses of metabolic, protein interaction, gene regulatory and neural networks all show heterogeneous connectivity, including heavy tail distributions and modular structure [814]. Heterogeneous networks, such as those with a broad connectivity distribution, are generally more difficult to analyze; inferring their detailed topology requires exceedingly high statistics. Power-law distributions often provide a good approximation to such networks, but statistical difficulties have led to some debate concerning their adequacy [15].

Understanding how heterogeneous connectivity profiles relate to properties of the associated dynamical system is a non-trivial theoretical challenge [1620]. In particular, heterogeneity makes these networks non-tractable by mean field methods [20]. In the current work we propose a different view on this problem, which interprets network topology as an effective feedback system. We develop a novel approximation for network ensembles with a broad connectivity distribution, based on the dominant role of a handful of hubs in a finite network. Identifying such hubs [21], understanding their impact on network functionality [22], and advancing the related theory [23] are all active research areas. Approximating a few hubs at the tail-end of the distribution by a single effective hub, we map the problem of heterogeneous network dynamics to that of a homogeneous network coupled to a single external node: each scale free network is approximated by a single lumped-hub connected to a homogeneous network of the remaining nodes (the bulk). Despite the simplification, this mapping retains properties of the original network ensemble. The new network is parameterized by hub strength and bulk connectivity.

We focus on the probability of random network ensembles to converge to attractors, and how this probability is modulated by ensemble topology. Convergence to fixed points is one of the simplest dynamical behaviors, and hence a suitable starting point for analysis. Moreover, stable fixed points have traditionally been linked to fundamental functional properties of networks [2426].

The lumped-hub approximation captures much of the behavior of heterogeneous network ensembles in terms of their probability of convergence to fixed points. The simplified ensemble can be analyzed by a combination of mean field approximations and feedback analysis. For networks with outgoing hubs, this analysis predicts a phase transition between converging and diverging dynamics as a function of the balance between hub strength and bulk connectivity, which is verified by numerical simulations. We show that the prediction holds on average for scale-free networks in a range of parameters relevant to real networks. Deviations from the theory stemming from large variability in scale-free ensembles as well as corrections to the mean field assumption are discussed.

Our results offer a novel framework for analyzing heterogeneous networks, which links insights from two different angles—internal dynamics of complex networks and the response to external input of simpler networks. The success of our mapping suggests that the role of outgoing hubs, in remarkable contrast to incoming ones [27], can be considered analogous to that of an external input in overcoming recurrent activity, suppressing chaotic dynamics, and ultimately driving the system to a stable fixed point. It highlights the crucial role of a small number of hubs in a finite network: their existence and relative strength is predictive of a given network’s dynamics, possibly more so than other quantitative features of the distribution from which the connectivity was sampled.


Dynamics with scale-free networks: Motivating observations

In what follows we analyze a model of nonlinear dynamics, often used to describe biological interactions such as in neuronal or genetic regulatory networks [24, 28]. Binary dynamic variables, si (1 ≤ iN), approximate the activity of N individual elements as being in one of two states: expressed (+1) or repressed (−1) genes, firing or quiescent neurons. The effect of elements on one another is described by a weighted sum of all incoming links:

where W is the N × N connectivity matrix with real entries giving the strengths of interactions, drawn from some random network ensemble, and

This is a special case of the widely studied Boolean networks [24, 29, 30], with the Boolean function chosen here to be the sign of the inputs’ weighted sum (for significance and alternative definitions of the sign function, see Materials and methods).

The ensemble from which W is drawn is a crucial ingredient of Eq (1). To study the effects of topology on the dynamics, we define W = TJ as a Hadamard (element-wise) product between a random topology matrix T and a random interaction strength matrix J. The 0/1 adjacency matrix T defines a directed graph, whose edges are sampled from specified distributions of incoming and outgoing connections (see Materials and methods). The strengths J are i.i.d. Gaussian variables with zero mean and standard deviation of unity.

For networks described by a directed graph, in general, incoming and outgoing degree distributions, Pin(k), Pout(k), need not be the same. In one case of special interest, gene regulatory networks, an empirical observation of such dissimilarity was reported. Outgoing connections were found to have a broad distribution consistent with a power-law: Pout(k) ∼ kγ , a scale-free (SF) distribution, while incoming connections were much more narrowly distributed around the average [31]. These statistical properties represent the biological observation that while any given gene is regulated by at most a handful of others, some transcription factors can regulate the expression of up to hundreds of genes in the cell. These ‘master regulators’ are of much interest in cell biology, and have been at the focus of many specific molecular studies [32]. In neuroscience, evidence was also found for a hierarchy of connectivity, with scale-free topology suggested in several contexts [10, 11]. In particular, “leader neurons” were identified and studied [33].

In a recent work describing gene regulatory networks by random interaction matrices, it was found that among various network topologies, those drawn from ensembles with outgoing hubs have a high probability to converge to attractors under exploratory adaptation [27]. Ensembles that showed efficient ability for adaptation also exhibited a high probability of convergence to a fixed point attractor in their intrinsic dynamics.

Particularly noteworthy in this context is the observation that for the transpose random network ensemble, where incoming connections are broadly distributed whereas outgoing connection are not, the abundance of fixed point attractors is decreased dramatically. This result is illustrated for the model of Eq (1) in Fig 1A (open circles). For each of the two ensembles—Scale-Free Out degree distribution (SFO) and the transposed SFI—the fraction of simulations converging to a fixed point after a given time interval was computed; this is used as an estimate of the convergence probability across the ensemble. The result is presented as a function of network size N for both ensembles. For SFO networks a considerable convergence probability is maintained up to the maximal network size tested, N = 5000, whereas for SFI networks this probability decreases rapidly and is negligible for networks with N > 1000. These results are consistent with those of [27], where a model with continuous dynamics was used, extending the observation on the two transposed ensembles to Boolean dynamics.

Dynamic properties of Scale-Free-Out (SFO, blue) vs. Scale-Free-In (SFI, red) network ensembles.
Fig 1
(A) Convergence probability is shown as a function of network size N. Filled circles show convergence to frozen cores larger than 90% of the network; open circles show convergence to fixed points. Each data point represents an ensemble average over 600 networks with random T, J and initial conditions. (B) Distribution of frozen core sizes (FC) relative to network size, across the SFO ensemble for N = 1500. The distribution is further zoomed near FC = 1 justifying the notion of quasi fixed point, defined in the text and elaborated in Materials and methods. Scale free parameters for both panels are γ = 2.2 and kmin = 1, Binomial parameters are N and p = k/N with k ∼ 5 (exact value of k is determined by the realization of scale free distribution. See Materials and methods for details).Dynamic properties of Scale-Free-Out (SFO, blue) vs. Scale-Free-In (SFI, red) network ensembles.

When simulating these Boolean dynamics, we noticed that networks often converge to a state where the vast majority of the nodes are fixed—“frozen”, while the rest continue to change (Fig 1B). From a biological perspective these are of interest as partially fixed states, and mathematically they have been studied in the context of general Boolean networks [30]. In what follows we define a quasi-fixed-point (QFP) as a state in which most (> 90%) of the nodes are frozen. Fig 1A (filled circles) demonstrates that the probability of convergence to a QFP behaves qualitatively similar to the corresponding probability for fixed points: the remarkable difference between SFO and SFI ensembles and its dependence on network size remain the same. This shows that the constraining effect of outgoing hubs on network dynamics is not a unique feature of fixed points, but is also present for weaker notions of convergence.

Our empirical observations rely on convergence of dynamics—and hence is affected both by the existence of fixed points and by their stability. In Materials and methods we show that networks from both SFO and SFI ensembles typically have O(1) fixed points. This suggests that the reported differences between ensembles stem from stability of fixed points. In the following sections we develop an approach which enables to better understand the stability properties of the two ensembles.

The lumped-hub approximation

A finite realization of a power-law connectivity distribution is dominated by a handful of hubs connected to a macroscopic fraction of the network. Fig 2A shows an empirical histogram of the outgoing connections in an SFO network with N = 1500, where the largest hubs each connect to ∼ 1000 nodes. To capture the properties of such a network we divide it into a small group of leading hubs, and the “bulk”—the rest of the network. The m largest hubs are lumped into one node, while the remaining nodes are substituted by a homogeneous network with binomially distributed incoming and outgoing degrees, preserving the average connectivity.

The lumped-hub approximation.
Fig 2
(A) In a finite network with SFO degree distribution, the m largest hubs (gray shading) are connected to a macroscopic fraction of the network. The histogram shows the number of outgoing connections from each node. (B) In the lumped-hub approximation of the same network, these are substituted by one effective hub and the rest of the nodes are approximated by a Binomial network. The connection strengths of the single hub to the bulk are determined by summing the strengths of the m lumped hubs. The mean degree of the bulk, kb, is retained in the mapping. Here N = 1500.The lumped-hub approximation.

We call this mapping the “lumped hub approximation”; the resulting network is characterized by the bulk mean connectivity kb and the pattern of connections ui from the lumped-hub to each node i in the bulk. These depend both on the original network topology and the lumping parameter m.

The connections ui are derived by summing the contributions of all the hubs that were lumped according to:

where we assume that hubs in the original SFO network are represented by rows 1, …, m of the connectivity matrix. In this construction, if more than one of the largest m hubs in the original network were connected to node i, then the connection between the lumped hub and this node is drawn from a distribution with appropriately scaled variance. We refer to this setting as Exact Pattern approximation, since the connectivity from the lumped hubs is fully preserved (as opposed to coarser approximations discussed below).

Despite the crudeness of the lumping approximation, we find that it preserves dynamic properties of scale free networks at the ensemble level. Fig 3 shows that the probability to converge to a QFP in the approximated networks follows the same dependence on network size as the original scale free ensembles. In particular, the significant difference between SFO and SFI, mapped to either an outgoing or incoming lumped hub, is captured by the approximation. These results are presented for a specific SFO ensemble, characterized by a fixed set of parameters (see figure caption for details).

Convergence probability of SFO/SFI lumped-hub approximation.
Fig 3
Probability of convergence to a Quasi Fixed Point (QFP), for SFO and SFI (blue and red filled circles respectively; same data as in Fig 1A). Their corresponding lumped-hub approximations are shown in blue and red open symbols with lumping parameter m = 4. Statistics and parameters as in Fig 1.Convergence probability of SFO/SFI lumped-hub approximation.

We next ask whether our lumping approximation can reproduce the range of different behaviors across SFO ensembles with various power-law distributions. To answer this question we map a broad range of SFO ensembles to lumped-hub networks and compare their dynamic properties. Fig 4 (top left) shows a scatter-plot of convergence probabilities for SFO networks vs. their lumped-hub approximations. Each point represents a specific SFO topology realization T, and the convergence probability is estimated as an average over realizations of connection strength J. It is seen that the convergence probability is largely retained in the mapping across a range of scale free parameter γ . In addition, the bottom-left panel of Fig 4 shows that the distribution of frozen core (FC) sizes, exhibiting a bi-modal shape dependent on ensemble parameters, is also very well captured by the lumping approximation. Comparisons of the autocorrelation and dimensionality of non-converging trajectories are given in Materials and methods, also supporting the similarity of the approximation to the original networks. These results demonstrate that our lumped-hub approximation describes dynamic properties of the SFO network ensembles across a broad range of parameters.

Lumped-hub approximation captures SFO network properties for a range of power-laws.
Fig 4
Within the lumped-hub approximation, we show three coarse grained versions of the hub outgoing connectivity pattern. (Left) Exact pattern of connectivity is maintained. (Middle) Number of nonzero connections is maintained, but their weights are replaced by a Gaussian distribution. (Right) A single Gaussian distribution describes all outgoing connections. (Top) Probability of converging to QFP for SFO networks and their respective lumped-hub approximations with increasingly coarse graining. For each SFO parameter γ (color code and legend), 125 topology realizations T were drawn, and the convergence probability was estimated from 50 realizations of the connection strength J and initial conditions. On each colored cluster of points, longer lines denote the direction of the linear fit to each γ, and the size of the perpendicular segment is the remaining variance after the fit. The identity line is added for reference (dashed black). Data points were jittered by a small uniform random noise (in range of ±0.01) to avoid their overlap due to discrete sampling. (Bottom) Cumulative distributions of frozen core sizes for SFO networks and their lumped-hub approximations for two power-law exponents: γ = 2.2, 2.8 (color coded as in the top panels). All scale free distributions have kmin = 1 (see Materials and methods for more details).Lumped-hub approximation captures SFO network properties for a range of power-laws.

Feedback analysis for the lumped-hub network ensemble

We have seen that despite its coarseness, the lumped hub approximation somehow captures an important feature of scale free networks that strongly affects their dynamics. Can it be harnessed to gain a deeper understanding of these networks?

In particular we are interested in the significantly larger convergence probability in the SFO ensemble compared to the corresponding SFI. The essence of this difference becomes intuitively clear in the lumped-hub approximation. In lumped SFO networks, the hub coherently drives a macroscopic fraction of the bulk nodes. In contrast, in lumped SFI networks the hub receives a large number of inputs but only drives a small fraction (≈ kb ) of the nodes. Previous work showed that intrinsically chaotic dynamics of homogeneous random networks can be suppressed by a coherent external drive [34]. Although here the hub is not external to the network, we argue that a similar suppression effect underlies convergence in SFO networks.

The lumped-hub network is described by two coupled equations of the bulk and the hub:

where s′ and W′ denote the bulk activity and connectivity respectively, the scalar h is the hub activity and the vectors u, v are the connections between the bulk and the lumped hub.

Our results so far relied on the outgoing connectivity from the lumped hub maintaining the Exact Pattern of the original connectivity, according to Eq (3). To facilitate the analysis, we consider an ensemble that retains statistics of the lumped hub connectivity rather than its exact pattern:

This ensemble is characterized by three parameters: mean degree of the bulk kb; the hub sparseness α, representing the fraction of its nonzero outgoing connections; and the variance σh2 representing its connection strengths. Assuming a Gaussian distribution of connections, this defines the Sparse Gaussian ensemble (Fig 4, center). This approximation is shown to predict both convergence probabilities (top) and frozen core statistics (bottom) nearly as well as the lumped hub with the Exact Pattern of connectivity.

Notably, a coarser model, which we refer to as the Dense Gaussian approximation, and where sparseness is not preserved, turns out to be dramatically less accurate. This ensemble, where the hub is connected to the entire network by outgoing connections drawn from a Gaussian distribution with variance ασh2, tends to systematically overestimate convergence probabilities (Fig 4, right). In what follows, we therefore analyze the Sparse Gaussian ensemble to assess the influence of the lumped hub on network dynamics.

Consider first an auxiliary setting—the “open loop” setting—where the hub is held at a fixed value. Here the conditions on column mean [20] are fulfilled and Mean Field analysis can be applied. Effectively, connections from the bulk of the network to the hub are discarded, and it remains connected only through its outgoing links (Fig 5A). The hub then acts as an input to the rest of the network, a simple situation that can be analyzed by Mean Field Theory. Formally, we replace the dynamic variable h by a fixed value (+ 1 without loss of generality), namely Eq (5) is substituted by h(t + 1) = 1. In this open loop setting, the regularizing effect of the fixed input competes with the recurrent dynamics of the bulk, and we ask under what conditions the system converges to a fixed point. Later, we relate the open-loop results back to a fully recurrent network including the hub as a dynamic variable [35].

Feedback analysis for the lumped-hub network ensemble.
Fig 5
(A) Conceptual scheme of feedback analysis. The lumped-hub network contains a single hub (gray) connected to a homogeneous network representing the bulk (green). In the open loop setting, the hub is clamped to a fixed value while its incoming connections from the bulk are disabled (×). (B) Mean Field Theory predicts a threshold value of σh above which chaos in the bulk network is suppressed and the system is driven to a stable fixed point. These transition lines are plotted in the (kb, σh) plane for various values of hub sparseness α. (C) The Probability of converging to a QFP for lumped-hub network ensembles computed as a function of hub strength σh, for open-loop (crossmarks, converging to 1) and closed-loop (circles, converging to 1/2). Two values of average bulk connectivity and sparseness are shown: kb = 3, α = 0.5 (blue) and kb = 5, α = 1 (red). MFT-predicted transition is depicted by vertical lines of the corresponding color. The solid lines are the open loop values divided by two. The numerical results confirm both the transition point and the factor of half due to closing the loop. N = 1500, and statistics is obtained from 239 realizations. (D) Same as (C), for the open loop case and as a function of kb for fixed (σh, α) with N = 5000. Blue: (σh = 1.0,α = 0.25), red: (σh = 2.12,α = 0.5). Statistics is obtained from 180 realizations. Gray regions in (B), (D) denote the limit of validity of MFT (See section Limits to the lumped hub approximation).Feedback analysis for the lumped-hub network ensemble.

To determine whether or not the system converges to a fixed point under input, we consider the time-evolution of the Hamming distance d(t) between two bulk trajectories s1(t) and s2(t). In the thermodynamic limit N → ∞, this evolution is deterministic and given by some function f(d ) (see Ref. [36] and Materials and methods for derivations):

The external drive h can suppress chaotic activity and drive the network to a stable state s′* if and only if d * = 0 is a stable fixed point of the mapping Eq (7). This, in turn, requires the derivative of d at the origin, denoted by ∇d, to be smaller than one. We show in Materials and Methods that, within Mean Field Theory (MFT), the condition for this can be expressed as a threshold on the hub strength,
with σcrit depending on the average bulk connectivity kb and the sparseness α . Fig 5B shows this critical line for several values of α in the (σh, kb) plane. Importantly, the lines are distinct for different hub sparsity α . As already seen in Fig 4 (right), ignoring sparsity leads to an oversimplified and inaccurate model.

Fig 5C illustrates the accuracy of our theory. Blue and red × symbols show simulated probability of convergence rising from zero to one as a function of σh for two sets of (kb, α ), with the critical values predicted by Eq (8) (vertical lines in corresponding colors). These results demonstrate both the qualitative prediction—a (smoothed) transition from zero to one in convergence probability, and the correct location of the transition at the predicted σcrit. The width of the transition is affected by two factors—the slope of the critical lines in the (σh, kb ) plane (Fig 5B) and the finite size of the system. We demonstrate these effects by using a larger network and traversing the plane along kb instead of σh , resulting in a much sharper transition (Fig 5D). We therefore conclude that MFT adequately predicts the open-loop convergence probability. Note that the critical driving strength increases indefinitely for large kb , in agreement with previous results [36, 37].

Returning now to the closed-loop network dynamics Eqs (4) and (5), where feedback to the hub is restored and it is a dynamic variable with incoming and outgoing connections, we ask whether a stable attractor under external drive is consistently maintained as an attractor of the recurrent dynamics. Intuitively one may argue that, with probability 1/2 with respect to the random connections in the ensemble, the clamped value will be consistent with the hub’s input; the driven steady-state will then also be an attractor of the full recurrent dynamics. This intuitive argument only refers to existence of a closed-loop fixed point, and does not guarantee its stability. Therefore, beyond the open-loop stability criterion Eq (8), the closed-loop setting requires to consider the recurrent dynamics of the hub. As long as the hub does not flip from its clamped value (h changes from +1 to −1), open and closed loop system dynamics coincide.

Consider a state s′, a small distance d = δ > 0 away from the fixed point s′*. For a stable fixed point we expect exponential convergence with a rate ∇d. Taking into account the finite size of the network, we estimate that convergence to a distance of less than one node (d < N−1) takes approximately

Preserving this stability in the full closed-loop network requires that the hub does not flip its sign during the time of convergence:

In Materials and Methods we show that the typical hub flipping time Tflip is approximately:

which implies that Eq (10) holds for ∇d far enough from the phase transition. Furthermore, it is argued that even in the proximity of the critical value ∇d ≈ 1 , the corrections due to closing the loop are smaller than the error of the mean field approximation made in the open loop analysis. Taken together, these arguments show that all effects of closing the loop will reduce the probability of the system converging to a stable fixed point by half. Fig 5C shows this is indeed the case. Circles depict closed-loop simulations of lumped-hub networks, and the matching lines show the smoothed result for open-loop divided by two.

Applicability of feedback analysis to finite SFO networks

We now return to evaluate the accuracy of our mean field approximation to the original scale free network ensembles. Making the connection between the theoretical results and empirical simulations entails two steps: first, for each scale free network realization we define the lumped-hub parameters (kb, σh, α); second, these parameters are used to compute ∇d in MFT (Materials and Methods), predicting convergence for values smaller than one. Importantly, the first step is not a one-to-one transformation between the parameters of the scale free ensemble and those of the lumped-hub ensemble: due to the large heterogeneity of scale free networks, different realizations generally result in widely varying lumped-hub parameters and consequently in broadly distributed values of ∇d . An example of this distribution is shown for one SFO ensemble in Fig 6A. In particular, lumped-hub parameters across the ensemble are such that ∇d is distributed on both sides of the phase transition (dashed vertical line, ∇d = 1). This wide distribution for fixed generative parameters, implies a lack of sharp phase transition when these parametes are varied.

Comparing lumped-hub MFT prediction with SFO simulations.
Fig 6
(A) Variability of lumped-hub parameters in a single scale free network ensemble: histogram shows the distribution of ∇d as computed from the lumped-hub parameters, for a scale free ensemble with γ = 2.4, kmin = 1. Black dashed line: MFT predicted phase transition, ∇d = 1. (B) Transition in convergence probability PQFP (blue symbols) computed for networks with different kmin values in range 0.5 ≤ kmin ≤ 2 and γ = 2.2, 2.4, 2.6. For each set of parameters 256 network realizations of size N = 1500 were simulated. Results were then binned by ∇d width of 0.05. Bins with less than ten samples were omitted from the plot (restricting the standard error to 0.16 per point). A short red line shows the data for the single ensemble of panel A. (C) Effect of network size on average lumped hub parameters k, σh, and α (sparseness; inset). Increasing network size causes a gradual shift towards the unstable phase (see Fig 5B), with kb affected most strongly. Colors show different γ. (D) Effect of network size on convergence probability: comparing simulation, lumped-hub approximation and theory. Empirical data points represent 256 Monte-Carlo realizations for scale free with γ = 2.4, kmin = 1.Comparing lumped-hub MFT prediction with SFO simulations.

Pooling together realizations of different power-law parameters opens a broader range of ∇d such that comparison to the theory is more meaningful. We simulated many SFO networks with a range of underlying power-law parameters and tested the prediction of a transition as a function of ∇d . Fig 6B compares the theoretical prediction for the convergence probability (dotted black line) with the pooled simulation results (blue symbols); despite the pooling of a highly variable set of simulations, ∇d is revealed as a crucial determinant of convergence. The smooth sigmoid obtained follows the transition around the predicted value of ∇d = 1. A decrease in convergence probability is found for very low ∇d values, in contrast to the prediction; this effect will be discussed below. For one SFO ensemble, the simulations follow the prediction but cover only a small range of ∇d (red symbols).

Empirical convergence probabilities for scale free networks also depend on network size (Fig 1). This dependence can be partially explained by the lumped-hub approximation: for a given scale free ensemble, we can follow the lumped-hub parameters averaged over the ensemble as a function of network size N . Fig 6C shows that the average (kb, σh, α) values progress towards the unstable regime as network size increases, with kb being influenced most strongly; these effects are consistent for several scale free ensembles (colors). The implications of this trend can be seen in the dependence of convergence probability on network size, shown in Fig 6D. Simulation results of scale free networks and their corresponding lumped-hub approximations are shown together with the theoretical prediction, showing that most of the N-dependency is captured correctly.

Limits to the lumped hub approximation

The lumped hub approximation relies on the possibility of decomposing the network into a hub and bulk nodes. In turn, the phase transition predicted for lumped-hub networks relies on this decomposition and its relation to homogeneous networks driven by an external input. This approximation holds for a range of parameter values of practical interest. To understand the limits of validity of the approximation, we consider several conditions that are required to hold.

First, the most highly connected nodes are being lumped together to one effective hub. The approximation thus neglects their internal dynamics, which should therefore be of little importance for convergence of the system. Ideally there should exist an interval of lumping numbers m across which ∇d remains largely unchanged. Fig 7A shows that this is the case for SFO networks with γ ≥ 2.4 (a change of roughly 0.1 compared to a width of more than 0.2 in Fig 6B). For γ = 2.2 we see that ∇d depends more strongly on the choice of m , compromising the stability of the approximation. Indeed, comparing the approximation to direct simulation results in Fig 7B, shows a large discrepancy between scale free and lumped-hub behaviors for γ = 2.2.

Limits of the approximation.
Fig 7
(A) Sensitivity of lumped hub approximation to the number of leading nodes being lumped is shown for 1000 realisations of scale-free topology with fixed parameters kmin = 1, N = 1500, and for three values of γ = 2.2, 2.4, 2.6. For each realization, the ∇d at m = 1 is taken as reference, from which the subsequent decrease in ∇d is measured. (B) Fig 6D is replicated for γ = 2.2, 2.8. In the large γ case, the lumped hub model holds, while mean field theory becomes invalid due to extreme sparseness. For low γ it is the lumped hub model per se that becomes inaccurate, due to sensitivity to choice of m.Limits of the approximation.

Second, the network must remain globally connected. For very small mean connectivity, networks tend to split into a number of disconnected components. In this situation, the probability of an individual realization converging to a fixed point drops, while realizations with multiple fixed points emerge (See Materials and methods). We believe that this effect is responsible for the drop in convergence probability observed at extremely low ∇d (Fig 6B).

In addition, for the approximation to be relevant, the effective external drive by the hub must be strong enough. As indicated in Fig 5B, for a sparsely connected hub (α = 0.1) the transition between convergent and chaotic regimes is almost independent of hub strength (the regime where the hub becomes irrelevant is marked by a gray area at Fig 5B and 5D). This can explain the failure of the theory for low α values, that arise for high γ values, as exemplified for γ = 2.8 in Fig 7B.

Finally, by lumping m > 1 dominant hubs into a single node, we are limited to a single fixed point of the open-loop system (up to symmetry of the hub activation). If several combinations of the m hubs’ states can stabilize the open-loop system, the possibility of multiple fixed point pairs in a single realization emerges. According to [38] and our Materials and Methods, such a possibility should reduce the probability of convergence. This can partially explain the slight drop below the factor of 1/2 in Fig 6B.


Broad connectivity profiles are a ubiquitous feature of naturally occurring networks. Understanding the dynamics arising from such topology, however, remains a great challenge. Our main hypothesis in this work was that, for such networks, a small number of nodes has a disproportionate effect on the entire system dynamics. We suggest a simplified model—the lumped-hub approximation—that makes this hypothesis explicit by lumping these nodes into a single hub and replacing the rest of the network by a homogeneous bulk.

We showed that the lumped-hub approximation predicts the behavior of scale free networks at the ensemble level, explaining the change in convergence probability as a function of power-law exponent and network size. Furthermore, the approximation is also helpful at the single realization level, indicating that a large part of the variability between scale-free networks stems from the variability in their effective lumped-hub strength and sparseness, and in the average bulk connectivity. In this respect, our proposed three-parameter description of topology is more predictive than the generative scale free parameters.

Our analysis was based on scale-free connectivity, but should apply for other topologies with broadly distributed out degree, for which mean field does not converge even with high moment corrections [20]. In fact, our results indicate that the existence of a small number of hubs connected to a macroscopic fraction of the network nodes, rather than the precise shape of the distribution tail, is a crucial factor in shaping network dynamical properties. The existence of one such outgoing (but not incoming) hub was sufficient not only to induce increased convergence to stable attractors, but also to endow the ensemble with an extremely high heterogeneity among realizations.

The simplicity of the lumped-hub network ensemble, with a homogeneous bulk network and a single hub connected to a fraction α of other nodes, allows analytic treatment of the problem. Specifically we applied feedback analysis, based on interpreting the hub as an external drive to the bulk, and then devising a mean field theory to determine convergence of the system. We then closed the loop while requiring consistency, allowing us to analytically estimate the probability for convergence. This analysis explains simply and intuitively why networks with outgoing hubs, and not those with incoming hubs, tend to converge to fixed points. We showed that, perhaps surprisingly, it is the stability of fixed points, rather than their abundance that underlies this effect.

In homogeneous networks, the average activity is a natural quantity of interest and provides the basis for mean-field approximations [36, 37]. For heterogeneous networks a low-dimensional approximation is still desirable, but the choice is no longer obvious. If the network is composed of homogeneous sub-networks a version of mean-field theory can be applied [39, 40]. Alternatively, a degree-dependent mean-field approach can be developed [41, 42]. Weighted averages of activity, guided by topology, can also serve to reduce the dimensionality of the dynamics [18, 19]. Our approach can be thought of as splitting the network and then approximating each part as homogeneous. This allows on one hand to tailor the approximation to specific heterogeneous realizations, but on the other hand to treat the approximate problem with mean-field and feedback tools.

The analytic treatment was developed for binary dynamic variables, but the empirical convergence properties and their dependence on outgoing hubs are shared also for a continuous variable model [27]. Therefore it is plausible that these properties reflect more generally the ability of outgoing hubs to promote convergence and to stabilize network dynamics.

We note the comparison of our results with previous work on Boolean networks with a randomly chosen Boolean function at each node [43], which found that the phase transition did not differ between SFO and SFI networks. By specifying a threshold function in our model, the effect of a hub becomes coherent among all the nodes that it influences. In contrast, when randomizing over Boolean functions, the effect of the hub has no such coherence; it is this coherence which promotes convergence in the model studied here.

From a biological viewpoint, our results illustrate how outgoing hubs—master regulators—can act as global coordinators of dynamics in gene regulatory networks. This role is expected to be realized under strong perturbations, in addition to the usual role of hubs as regulators of local specific circuits. Indeed, recent experimental work in bacteria suggests an organizing role for hubs in the emergence of gene expression patterns following strong rewiring perturbations [44, 45]. This idea, suggesting a dual role for the structure of gene regulation networks depending on conditions, awaits further experimental and theoretical investigation.

Materials and methods

Suppression of chaos in the open-loop setting

We consider a lumped hub network ensemble, where the effective hub is clamped at a fixed value and acts as an input to the rest of the nodes. It is connected to the network by a vector u, covering a fraction α of the network, with connection strengths drawn from a Gaussian distribution with zero mean and standard deviation σh. Formally, we determine whether the input suppresses chaos by following the time evolution of the normalized Hamming distance between two trajectories s1(t) and s2(t), denoted by d(t). Suppose that d(t0) = 1/N, and assume, without loss of generality, that the states of interest s1, s2 differ at the first bit: s11s12 while sj1=sj2 for 2 ≤ jN. At time t0 + 1 the distance d will be determined by the number of rows in W for which the first element changes its sign:

Equivalently, one may consider the derivative at the origin of the normalized distance d, given by:
with the condition for stability being ∇d < 1. The inequality in curly brackets can be realized only if Wi1 is nonzero; moreover the other matrix element can be grouped according to the number of nonzero elements in each row. This grouping results in the following convenient way to rewrite the last equation:
with terms P1,2,3 defined as:
Here z1, z2 are i.i.d. normal Gaussian variables, as follows from the mean field approximation, and i is an arbitrary row. The terms P1 and P2 are determined by the network structure: P1 is simply the sparsity of the Binomial network which makes up the bulk nodes,
and the probability of a row in the network having exactly k′ nonzero entries is:
The only term which depends on the actual realization of weights is P3. To evaluate it we note that for a positive scalar η:
To derive this expression, one notes that in the |z1|, |z2| plane, we are interested in the portion of probability mass that lies below the line |z2|=|z1|η. Since the probability density in this plane only depends on the distance from the origin, we simply compute the angle between this line and the |z1 | axis and Eq (20) follows.

Combining all the terms together, and assuming large N, we arrive at:


For α = 1 and large kb, the sum over k′ can be approximated by its expected value, and tanh by its argument leading to:

This expression accounts for the asymptotically linear relation between k and σ at the phase transition curve in Fig 5B. Conversely, for sparse hubs with α < 1 there exists a kmax above which we have ∇d > 1 for any σh.

Quasi fixed points

Along with a possibility of converging to a stable fixed point another scenario exists where most network nodes are fixed but a small fraction is still toggling.

More specifically, a situation where all nodes except a small, interconnected clique are frozen becomes typical. We refer to such a regime as a quasi fixed point (QFP). A threshold to define a QFP was set empirically at a value of 0.9 where a rise in probability of frozen cores starts (Fig 1B).

Fixed points are always there, regardless of topology

Defined by equation Eq (1), our system turns out to be a special case of the Theorem in Ref. [38] which states that the expectation EM of number of fixed points M in random Boolean networks is one , subject to a condition on the distribution of random Boolean functions that holds in our setting. Specifically, to comply with [38] an ensemble of random Boolean functions ϕi(s) defining the network dynamics via si(t + 1) = ϕi(s(t)), must have a set of neutral links, the removal of which renders the network acyclic. Here a link ji between a pair of nodes is said to be neutral iff

for any Boolean function g . In our case of a threshold function, the condition Eq (23) is fulfilled immediately because
for Gaussian matrices W1 and W2 differing by inversion of a single element.

Unfortunately, the expectation, without any higher moments known, does not provide information about a typical case. Depending on the topology T, a typical realization of W might have M = O(1). Alternatively, there could be a few exceptional realizations with a huge number of fixed points while typically M = 0.

For example, for T = I there are M = 0 fixed points with probability 1 − 2N and M = 2N with probability 2N. Conversely if T is a cyclic graph (e.g. cyclic permutation) then with probability of half the network has two fixed points, and it has zero fixed points in the complementary case. Although it is not immediately clear which scenario is relevant for scale free networks, for lumped hub networks some insights are available.

For a lumped hub network with a strong hub, such that ∇d < 1 in the open loop setting, global convergence to a stable fixed point s′ = s* is expected. If the consistency condition associated with relaxation of clamping h(t ) in Eq (5) is met, then the closed loop system inherits the fixed point at s = s*. Moreover, another fixed point emerges at s = −s* corresponding to a drive of −1 in the open loop system. Since the aforementioned condition is met with probability half we have an expected number of 2×12=1 fixed points in agreement with [38].

In the opposite extreme case of zero drive from the hub, it follows from mean field theory that the network exhibits chaotic dynamics with any pair of state-space trajectories becoming asymptotically orthogonal. In particular, to comply with MFT, distinct fixed points should be either orthogonal to- or inverse of one another. We assume, without a rigorous proof, that orthogonality approximately implies independence. Namely, events of having a fixed point at s = s1 and at s = s2 are statistically independent for s1s2. Now, we recall that fixed points appear in pairs. Hence the number of pairs of fixed points is distributed binomially with parameter n=2N2=2N-1 of possible pairs ±s* of a state and its inverse, and a probability p = 2N that a particular pair of points is fixed. In the limit of N → ∞ this distribution converges to Poisson with parameter λ=pn=12:


To validate this theory numerically, we performed an exhaustive search for fixed points in small, fully connected networks with 16 ≤ N ≤ 22 nodes. For N > 22 an exhaustive search was not feasible due to computational constraints. The results, shown in Fig 8, depict an approximate match to the theory which is not perfect, with inaccuracies exceeding the standard error. We attribute these inaccuracies to the finite size of the networks, emphasizing that orthogonality of chaotic trajectories upon which our theory is built, is not achieved, even approximately, for values of N small enough for exhaustive search. Future work may include numerical experiments with larger and sparser networks using advanced search algorithms e.g. the method of [46].

Probability of number of fixed points in Boolean networks.
Fig 8
’+’—theoretical prediction. ∘ and ⋄—empirical estimates for networks of size N = 16 and 19 based on 104 and 103 random network realizations respectively.Probability of number of fixed points in Boolean networks.

For intermediate values of σh the picture is more complicated. Chaos is not suppressed but the symmetry between ±s is broken by the drive (in the open loop setting only, since once the loop is closed symmetry is restored). For σh ≪ 1 one might repeat our reasoning for σh = 0, this time with single fixed points rather than pairs, and therefore with Eq (25) being replaced by:

Numerical observations on small networks reported in Table 1 tell us that while a sharp transition in the distribution of fixed points is observed once a drive is introduced, the distribution of M does not immediately match Eq (26). Here again we tend to attribute such a discrepancy with theory to finite size effects compromising the fixed point orthogonality.

Table 1
Random Gaussian drive is scaled by σh. Poisson distribution with parameter λ=1,12 are provided for reference.
Numerical estimates for the distribution of M are shown as a function of external drive strength for an ensemble of dense random networks of size N = 16 with Gaussian i.i.d elements.
10.6070.3040.0760.0130.0022 × 10−4

To conclude, we present strong evidence, that the event of having a fixed point in the dynamical system (1) has probability of order one and it is the fixed point stabilities rather than their existence that determines the convergence probabilities of dynamical numeric simulations like in Fig 1.

Remarkably, this result is opposed to a case of continuous activation function (e.g. tanh) where an exponentially large amount of fixed points appear in the chaotic regime [47]. An in-depth analysis of these differences is beyond the scope of the current work.

Estimation of hub flipping time in perturbed closed loop network

To obtain a, fairly coarse, lower bound on Tflip one may argue as follows: the hub has approximately kh = mk〉 = O(1) incoming connections whose weighted sum determines its sign. Let Hin denote the set of these nodes. Given a random perturbation of magnitude δ ≪ 1 to the fixed point s* the probability for a single node in Hin to be flipped is

and the probability of q flips is pqp1q=O(δqN-q) which can be neglected for small δ. Finally, by the arguments developed above, the probability that the hub will flip due to a single node flip is 2πarctan1kh and consequently, the flip rate is given by:
This corresponds to a typical flip time of:
Hence, the stability condition Eq (10) clearly holds for −log ∇d of order one.

Even in the critical regime, where log ∇ → 0 is approached, condition Eq (10) translates into

which for our typical setting of N = 1500, < kh >≈ 20, and a reasonable small perturbed fraction of δ = 0.1 of the nodes translates into a requirement of:
This is very close to criticality compared to other effects that may dominate inaccuracy of mean field calculations even in open loop setting. (Compare to Fig 5 that depicts open loop and closed loop setting, and Fig 6 which includes plot of convergence probability vs. ∇d).

Beyond convergence: Dynamical properties of lumped hub approximation

Most of the analysis in the main text focused on the probability of convergence to quasi fixed points. Here we demonstrate that even when the networks do not converge, some of their dynamical properties are still captured by the lumped hub approximation. First, we computed the participation ratio as a measure of the dimensionality of dynamics in the non-converging case. The correspondence between SFO and lumped-hub approximation is seen in Fig 9, similar to the match of convergence fraction shown in Fig 4.

Trajectory dimensionality captured by lumped hub approximation.
Fig 9
The participation ratio of trajectories, defined as (trC)2tr(C2) with C the covariance matrix of the trajectories s(t ), is compared between scale free network and Lumped-hub approximations of increasingly coarser graining. Notations and methods are similar to Fig 4.Trajectory dimensionality captured by lumped hub approximation.

Second, we computed the autocorrelation functions of dynamic trajectories. Like other statistics, they greatly vary from one realization to the next in SFO ensembles. Nevertheless, average correlations functions—also, like other statistics—show a favorable comparison to the corresponding lumped-hub ensemble (Fig 10). This is quantified by the exponential fit to the autocorrelation function, after selecting only trajectories with low asymptotic autocorrelation.

Autocorrelation of SFO ensembles vs. their lumped hub approximation.
Fig 10
Top: Log of mean autocorrelation is shown for scale free (left) and corresponding lumped-hub (right) networks. The mean is obtained by averaging only over trajectories with asymptotic autocorrelation below 0.3. Insets show 100 individual trajectories obtained for γ = 2.4. Bottom: Exponential fits for the autocorrelation functions in the range 0 ≤ t ≤ 5 (γ color coded accordingly to the top panel.Autocorrelation of SFO ensembles vs. their lumped hub approximation.

Constructing samples of heterogeneous network ensembles

To construct network samples from an ensemble with a SFO dgree distribution, the adjacency matrix is constructed by first randomly sampling a sequence of N degrees from a scale-free distribution, and assigning a degree, k, from that sequence randomly to each node in the network. For each such node a set of k random outgoing connection is chosen. This procedure results in a scale-free outgoing distribution and Binomial incoming distribution as there are no constraint on the incoming distribution and the choice of incoming degrees is purely random. SFI ensembles where constructed by transposing SFO networks created as described above.

The scale-free sequences are obtained by a discretization to the nearest integer of the continuous Pareto distribution P(k)=(γ-1)kmin(γ-1)kγ. This is implemented by the continuous MATLAB Generalized Pareto random number generators with Generalized Pareto parameters k = 1/(γ − 1), σ = kmin/(γ − 1) and θ = kmin.

Sign function

Networks with sparse connectivity often have entire rows that are zero; namely, nodes are created that receive zero input. With the definition Eq (2), such nodes retain their value of zero and therefore effectively do not participate in the dynamics. The effective network thus has slightly modified size and statistical parameters. We have verified by simulation that all our results are weakly affected by this modification. In particular, the changes induced in ∇d are smaller than the finite-size effects.

Other definitions of the sign function, for example

do not effectively prune these source nodes. Instead, these source nodes become effectively equivalent to an external drive and thus increase the stability of the dynamics.

Probability of convergence to fixed points and quasi fixed points (QFPs)

The networks’ frozen cores were determined by running the dynamics for an initial convergence period of 4000 time steps, with 10% of networks nodes being updated at each step, and then measuring the fraction of nodes which were frozen within an additional time interval of 1000 steps (altogether 5000 time steps). Some simulations were run for double the time steps (a of total 10000 time steps, 2000 of which used for determining frozen core). No detectable difference was observed with such longer times. QFPs were defined as states with a frozen core larger than 90%, while fixed points are networks with a frozen core of 100%. The probability of convergence to a fixed point or QFP was calculated by simulating the dynamics of a an ensemble of 500 networks, unless stated otherwise, and measuring the fraction of networks in the ensemble which reached the relevant frozen core criterion.



    LH Hartwell, JJ Hopfield, S Leibler, AW Murray. . From molecular to modular cell biology. Nature. 1999;402(6761supp):, pp.C47.


    AL Barabasi, ZN Oltvai. . Network biology: understanding the cell’s functional organization. Nature Reviews Genetics. 2004;5(2):, pp.101, doi: 10.1038/nrg1272


    U Alon. . Network motifs: theory and experimental approaches. Nature Reviews Genetics. 2007;8(6):, pp.450, doi: 10.1038/nrg2102


    SH Strogatz. . Exploring complex networks. nature. 2001;410(6825):, pp.268, doi: 10.1038/35065725


    R Albert, AL Barabási. . Statistical mechanics of complex networks. Reviews of Modern Physics. 2002;74(1):, pp.47, doi: 10.1103/RevModPhys.74.47


    ME Newman. . The structure and function of complex networks. SIAM Review. 2003;45(2):, pp.167–256. , doi: 10.1137/S003614450342480


    S Boccaletti, V Latora, Y Moreno, M Chavez, DU Hwang. . Complex networks: Structure and dynamics. Physics Reports. 2006;424(4-5):, pp.175–308. , doi: 10.1016/j.physrep.2005.10.009


    G Caldarelli. Scale-free networks: complex webs in nature and technology. Oxford University Press; 2007.


    RM D’Souza. . Complex networks: Structure comes to random graphs. Nature Physics. 2009;5(9):, pp.627, doi: 10.1038/nphys1390


    D Eytan, S Marom. . Dynamics and effective topology underlying synchronization in networks of cortical neurons. Journal of Neuroscience. 2006;26(33):, pp.8465–8476. , doi: 10.1523/JNEUROSCI.1627-06.2006


    G Buzsáki, K Mizuseki. . The log-dynamic brain: how skewed distributions affect network operations. Nature Reviews Neuroscience. 2014;15(4):, pp.264, doi: 10.1038/nrn3687


    A Statman, M Kaufman, A Minerbi, NE Ziv, N Brenner. . Synaptic Size Dynamics as an Effectively Stochastic Process. PLoS Comput Biol. 2014.


    B Barbour, N Brunel, V Hakim, JP Nadal. . What can we learn from synaptic weight distributions?TRENDS in Neurosciences. 2007;30(12):, pp.622–629. , doi: 10.1016/j.tins.2007.09.005


    Jn Teramae, T Fukai. . Computational implications of lognormally distributed synaptic weights. Proceedings of the IEEE. 2014;102(4):, pp.500–512. , doi: 10.1109/JPROC.2014.2306254


    G Lima-Mendez, J van Helden. . The powerful law of the power law and other myths in network biology. Molecular BioSystems. 2009;5(12):, pp.1482–1493. , doi: 10.1039/b908681a


    B Barzel, AL Barabási. . Universality in network dynamics. Nature physics. 2013;9(10):, pp.673–681. , doi: 10.1038/nphys2741


    C Castellano, R Pastor-Satorras. . Relating topological determinants of complex networks to their spectral properties: Structural and dynamical effects. Physical Review X. 2017;7(4):, pp.041024.


    E Laurence, N Doyon, LJ Dubé, P Desrosiers. . Spectral dimension reduction of complex dynamical networks. Physical Review X. 2019;9(1):, pp.011042.


    J Gao, B Barzel, AL Barabási. . Universal resilience patterns in complex networks. Nature2016;530(7590):, pp.307, doi: 10.1038/nature16948


    F Farkhooi, W Stannat. . Complete mean-field theory for dynamics of binary recurrent networks. Physical Review Letters. 2017;119(20):, pp.208301, doi: 10.1103/PhysRevLett.119.208301


    CY Lin, CH Chin, HH Wu, SH Chen, CW Ho, MT Ko. . Hubba: hub objects analyzer—a framework of interactome hubs identification for network biology. Nucleic Acids Research. 2008;36(suppl_2):, pp.W438–W443. , doi: 10.1093/nar/gkn257


    E Zotenko, J Mestre, DP O’Leary, TM Przytycka. . Why do hubs in the yeast protein interaction network tend to be essential: reexamining the connection between the network topology and essentiality. PLoS Computational Biology. 2008;4(8):, pp.e1000140, doi: 10.1371/journal.pcbi.1000140


    U Harush, B Barzel. . Dynamic patterns of information flow in complex networks. Nature Communications. 2017;8(1):, pp.2181, doi: 10.1038/s41467-017-01916-3


    SA Kauffman. The origins of order: Self-organization and selection in evolution. Oxford University Press, USA; 1993.


    S Huang, G Eichler, Y Bar-Yam, DE Ingber. . Cell fates as high-dimensional attractor states of a complex gene regulatory network. Physical Review Letters. 2005;94(12):, pp.128701, doi: 10.1103/PhysRevLett.94.128701


    DJ Amit. Modeling brain function: The world of attractor neural networks. Cambridge university press; 1992.


    HI Schreier, Y Soen, N Brenner. . Exploratory adaptation in large random networks. Nature Communications. 2017;8.


    JJ Hopfield. . Neural networks and physical systems with emergent collective computational abilities. Proceedings of the national academy of sciences. 1982;79(8):, pp.2554–2558. , doi: 10.1073/pnas.79.8.2554


    M Aldana, S Coppersmith, LP Kadanoff. Boolean dynamics with random couplings In: Perspectives and Problems in Nolinear Science. Springer; 2003 p. , pp.23–89.


    B Drossel. . Random boolean networks. Reviews of Nnonlinear Dynamics and Complexity. 2008;1:, pp.69–110.


    N Guelzim, S Bottani, P Bourgine, F Képès. . Topological and causal structure of the yeast transcriptional regulatory network. Nature Genetics. 2002;31(1):, pp.60, doi: 10.1038/ng873


    TL Davis, I Rebay. . Master regulators in development: Views from the Drosophila retinal determination and mammalian pluripotency gene networks. Developmental Biology. 2017;421(2):, pp.93–107. , doi: 10.1016/j.ydbio.2016.12.005


    JP Eckmann, S Jacobi, S Marom, E Moses, C Zbinden. . Leader neurons in population bursts of 2D living neural networks. New Journal of Physics. 2008;10(1):, pp.015011, doi: 10.1088/1367-2630/10/1/015011


    K Rajan, LF Abbott, H Sompolinsky. . Stimulus-dependent suppression of chaos in recurrent neural networks. Phys Rev E. 2010;82:, pp.011903, doi: 10.1103/PhysRevE.82.011903


    A Rivkind, O Barak. . Local dynamics in trained recurrent neural networks. Physical Review Letters. 2017;118(25):, pp.258101, doi: 10.1103/PhysRevLett.118.258101


    E Wallace, HR Maei, PE Latham. . Randomly Connected Networks Have Short Temporal Memory. Neural Computation. 2013;25(6):, pp.1408–1439. , doi: 10.1162/NECO_a_00449


    H Sompolinsky, A Crisanti, HJ Sommers. . Chaos in Random Neural Networks. Phys Rev Lett. 1988;61:, pp.259–262. , doi: 10.1103/PhysRevLett.61.259


    F Mori, A Mochizuki. . Expected Number of Fixed Points in Boolean Networks with Arbitrary Topology. Phys Rev Lett. 2017;119:, pp.028301, doi: 10.1103/PhysRevLett.119.028301


    M Stern, H Sompolinsky, LF Abbott. . Dynamics of random neural networks with bistable units. Phys Rev E. 2014;90:, pp.062710, doi: 10.1103/PhysRevE.90.062710


    J Aljadeff, M Stern, T Sharpee. . Transition to Chaos in Random Networks with Cell-Type-Specific Connectivity. Phys Rev Lett. 2015;114:, pp.088101, doi: 10.1103/PhysRevLett.114.088101


    SN Dorogovtsev, AV Goltsev, JF Mendes. . Critical phenomena in complex networks. Reviews of Modern Physics. 2008;80(4):, pp.1275, doi: 10.1103/RevModPhys.80.1275


    R Pastor-Satorras, C Castellano, P Van Mieghem, A Vespignani. . Epidemic processes in complex networks. Reviews of Modern Physics. 2015;87(3):, pp.925, doi: 10.1103/RevModPhys.87.925


    M Aldana. . Boolean dynamics of networks with scale-free topology. Physica D: Nonlinear Phenomena. 2003;185(1):, pp.45–66. , doi: 10.1016/S0167-2789(03)00174-X


    R Baumstark, S Hänzelmann, S Tsuru, Y Schaerli, M Francesconi, FM Mancuso, et al. The propagation of perturbations in rewired bacterial gene networks. Nature Communications. 2015;6:, pp.10105, doi: 10.1038/ncomms10105


    Y Murakami, Y Matsumoto, S Tsuru, BW Ying, T Yomo. . Global coordination in adaptation to gene rewiring. Nucleic Acids Research. 2015;43(2):, pp.1304–1316. , doi: 10.1093/nar/gku1366


    N Berntenis, M Ebeling. . Detection of attractors of large boolean networks via exhaustive enumeration of appropriate subspaces of the state space. BMC Bioinformatics. 2013;14(1):, pp.361, doi: 10.1186/1471-2105-14-361


    G Wainrib, J Touboul. . Topological and dynamical complexity of random neural networks. Phys Rev Lett. 2013;110(11):, pp.118101, doi: 10.1103/PhysRevLett.110.118101

21 Nov 2019

Dear Dr Brenner,

Thank you very much for submitting your manuscript 'Scale free topology as an effective feedback system' for review by PLOS Computational Biology. Your manuscript has been fully evaluated by the PLOS Computational Biology editorial team and in this case also by independent peer reviewers. The reviewers appreciated the attention to an important problem, but raised some substantial concerns about the manuscript as it currently stands. While your manuscript cannot be accepted in its present form, we are willing to consider a revised version in which the issues raised by the reviewers have been adequately addressed. We cannot, of course, promise publication at that time.

Please note while forming your response, if your article is accepted, you may have the opportunity to make the peer review history publicly available. The record will include editor decision letters (with reviews) and your responses to reviewer comments. If eligible, we will contact you to opt in or out.

Your revisions should address the specific points made by each reviewer. Please return the revised version within the next 60 days. If you anticipate any delay in its return, we ask that you let us know the expected resubmission date by email at Revised manuscripts received beyond 60 days may require evaluation and peer review similar to that applied to newly submitted manuscripts.

In addition, when you are ready to resubmit, please be prepared to provide the following:

(1) A detailed list of your responses to the review comments and the changes you have made in the manuscript. We require a file of this nature before your manuscript is passed back to the editors.

(2) A copy of your manuscript with the changes highlighted (encouraged). We encourage authors, if possible to show clearly where changes have been made to their manuscript e.g. by highlighting text.

(3) A striking still image to accompany your article (optional). If the image is judged to be suitable by the editors, it may be featured on our website and might be chosen as the issue image for that month. These square, high-quality images should be accompanied by a short caption. Please note as well that there should be no copyright restrictions on the use of the image, so that it can be published under the Open-Access license and be subject only to appropriate attribution.

Before you resubmit your manuscript, please consult our Submission Checklist to ensure your manuscript is formatted correctly for PLOS Computational Biology: Some key points to remember are:

- Figures uploaded separately as TIFF or EPS files (if you wish, your figures may remain in your main manuscript file in addition).

- Supporting Information uploaded as separate files, titled Dataset, Figure, Table, Text, Protocol, Audio, or Video.

- Funding information in the 'Financial Disclosure' box in the online system.

While revising your submission, please upload your figure files to the Preflight Analysis and Conversion Engine (PACE) digital diagnostic tool, PACE helps ensure that figures meet PLOS requirements. To use PACE, you must first register as a user. Then, login and navigate to the UPLOAD tab, where you will find detailed instructions on how to use the tool. If you encounter any issues or have any questions when using PACE, please email us at

To enhance the reproducibility of your results, we recommend that you deposit your laboratory protocols in, where a protocol can be assigned its own identifier (DOI) such that it can be cited independently in the future. For instructions see here

We are sorry that we cannot be more positive about your manuscript at this stage, but if you have any concerns or questions, please do not hesitate to contact us.


Alexandre V. Morozov, Ph.D.

Associate Editor

PLOS Computational Biology

Mark Alber

Deputy Editor

PLOS Computational Biology

A link appears below if there are any accompanying review attachments. If you believe any reviews to be missing, please contact immediately:


Reviewer's Responses to Questions

Comments to the Authors:

Please note here if the review is uploaded as an attachment.

Reviewer #1: Rivkind et al. submitted an interesting numerical study of the scale-free networks of interacting binary units (Boolean). The authors first show the convergence probability (over finite realizations) to a fixed point is much lower for Scale-Free-Out (SFO) compared with Scale-Free-In (SFI) degree networks. The second result of their paper is the illustration of "lumped hub approximation" for the SFO network. Remarkably, lumped hub approximation captures the distribution of frozen states and convergence to probability to a fixed point. After that, the authors analyze a Mean-Field Theory (MFT) for the lumped hub approximation to give more insight for the effect of hubs on the recurrent dynamics. Here, the authors study two closed and open-loop interactions between the hub and bulk part of the approximate network. Furthermore, they investigate the applicability of the MFT for the open/closed loop approximation.

I believe this is undoubtedly nice work, and it is relevant for the various subfields of computational biology and theoretical neuroscience. Therefore, I support its publication in PLoS CB.

However, in my opinion, the authors have to address three major points in the current manuscript to secure the desired impact:

First, it is unclear to me whether SFO networks converge to point-wise deterministic dynamics. The variance of the column sum in these networks is diverging, and as a result, it breaks the conditions for the existence of a deterministic MFT, as described in Farkhooi and Stannat (PRL, 2017). This has to be discussed in the revision and the notion of a fixed point in this context has to be appropriately defined.

Second, the comparison between the SFO and lumped hub approximation can be improved. I'd suggest two different comparisons in temporal statistics of the networks, and time-averaged properties should be given. This comparison can be used to motivate further analysis in the paper.

Third, the open/closed-loop analysis using MFT can be improved. It is not so clear for me (in the way that it is written in the main text) if suppression of chaos and the phase transition depends on the approximation itself. I wonder if authors could provide auto-correlations of the system before phase transition for both the approximation and full SFO networks.

Reviewer #2: The authors investigate how the topology determines the controllability of a dynamic system, and investigate specifically the impact of the high-degree nodes in a scale-free topology. Scale-free topologies are characterized by many nodes with only one or a few connections, some nodes with intermediate number of connections, and a few nodes with very many connections. They bundle the highly connected nodes into a single one, and replace the remaining part of the scale-free network by a random network. This simplification facilitates analytical mean-field treatment. They then show under which conditions the network converges – or does not converge to a fixed point.

I am divided with respect to the impact and generality of the results. On the one hand, it is a very elegant and creative framework, and the results are convincing. On the other hand, the approximations for the mean-field are rough, and the dynamics on the network is very simple, including the fixed-point as an absorbing state. Thus, I am not yet convinced that the framework is of general biological relevance.

Specific comments:


It remains a question whether the “lumped-hub” approximation is a good one. Taking e.g. figure 4A, there is some correlation, but I guess almost any coarse-graining of the topology would result in some correlation here. Can you argue that your choice of coarse-graining is a good one in general?

Specifically, why lumping together four, not three or five of the most connected nodes? Why randomizing the other topology, not just keeping it? How much does it impact the results?


Abstract: “Based on the observation that in finite networks a small number of hubs have a

disproportionate effect on the entire system, we construct an approximation by lumping

these nodes into a single effective hub, which acts as a feedback loop with the rest of the


Only fairly late it becomes clear, why you combine the set of ‘most important nodes’ to one node, instead of averaging exactly in the opposite manner, namely over all the non-important hubs:

It seems to be much more intuitive, to reduce e.g. all the nodes with degree 1, which is the most frequent degree in a scale-free network, and then potentially also remove the other low-degree nodes, until a treatable topology is found.

It would be good to make the aim more clear at an earlier point.

The terms “source” and “sink” might help the description for the outgoing/ incoming hubs.

Please discuss whether scale-free distributions are really a reasonable description for neural topology: Is it really the case that the most frequent projection pattern of a neuron is to connect to just one other neuron? (p(k=1)>=p(k=N) for all N. – I doubt so.

Equation 1: I’d suggest to write out the sign function.

Figure 4: I’d suggest to use CDFs instead, and plot both PDFs in a single panel.

Fig. 5: “Mean field theory predicts a phase transition at \\sigma_crit”.

Could you elaborate about the position of the \\sigma_crit? I would have expected it at about \\sigma_h=1.2 for the red, and close to zero for the blue condition, i.e. at the point where QFPs start to emerge?


Have all data underlying the figures and results presented in the manuscript been provided?

Large-scale datasets should be made available via a public repository as described in the PLOS Computational Biologydata availability policy, and numerical data that underlies graphs or summary statistics should be provided in spreadsheet form as supporting information.

Reviewer #1: Yes

Reviewer #2: No:


PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: No

Reviewer #2: No

27 Jan 2020

Submitted filename: Replies_Jan 23_v2.pdf

26 Mar 2020

Dear Professor Brenner,

We are pleased to inform you that your manuscript 'Scale free topology as an effective feedback system' has been provisionally accepted for publication in PLOS Computational Biology.

Before your manuscript can be formally accepted you will need to complete some formatting changes, which you will receive in a follow up email. A member of our team will be in touch with a set of requests.

Please note that your manuscript will not be scheduled for publication until you have made the required changes, so a swift response is appreciated.

IMPORTANT: The editorial review process is now complete. PLOS will only permit corrections to spelling, formatting or significant scientific errors from this point onwards. Requests for major changes, or any which affect the scientific understanding of your work, will cause delays to the publication date of your manuscript.

Should you, your institution's press office or the journal office choose to press release your paper, you will automatically be opted out of early publication. We ask that you notify us now if you or your institution is planning to press release the article. All press must be co-ordinated with PLOS.

Thank you again for supporting Open Access publishing; we are looking forward to publishing your work in PLOS Computational Biology. 

Best regards,

Alexandre V. Morozov, Ph.D.

Associate Editor

PLOS Computational Biology

Mark Alber

Deputy Editor

PLOS Computational Biology


Reviewer's Responses to Questions

Comments to the Authors:

Please note here if the review is uploaded as an attachment.

Reviewer #1: At first, I apologize for my delay in responding to the revision of the submitted manuscript. I was rather ill in the last three weeks.

I read the revision and the authors response completely and now I am confident that their work is interesting and very relevant for the readership of PLoS CB. In their version, they provide convening additional results on the temporal dynamics of the network activities. Also, they revised the text to be more accessible to readers.

Therefore, I recommend the revision for publication in PLOS CB.


Have all data underlying the figures and results presented in the manuscript been provided?

Large-scale datasets should be made available via a public repository as described in the PLOS Computational Biologydata availability policy, and numerical data that underlies graphs or summary statistics should be provided in spreadsheet form as supporting information.

Reviewer #1: Yes


PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: No

4 May 2020


Scale free topology as an effective feedback system

Dear Dr Brenner,

I am pleased to inform you that your manuscript has been formally accepted for publication in PLOS Computational Biology. Your manuscript is now with our production department and you will be notified of the publication date in due course.

The corresponding author will soon be receiving a typeset proof for review, to ensure errors have not been introduced during production. Please review the PDF proof of your manuscript carefully, as this is the last chance to correct any errors. Please note that major changes, or those which affect the scientific understanding of the work, will likely cause delays to the publication date of your manuscript.

Soon after your final files are uploaded, unless you have opted out, the early version of your manuscript will be published online. The date of the early version will be your article's publication date. The final article will be published to the same URL, and all versions of the paper will be accessible to readers.

Thank you again for supporting PLOS Computational Biology and open-access publishing. We are looking forward to publishing your work!

With kind regards,

Bailey Hanna

PLOS Computational Biology | Carlyle House, Carlyle Road, Cambridge CB4 3DN | United Kingdom | Phone +44 (0) 1223-442824 | | @PLOSCompBiol free topology as an effective feedback system&author=Alexander Rivkind,Hallel Schreier,Naama Brenner,Omri Barak,Alexandre V. Morozov,Mark Alber,Alexandre V. Morozov,Mark Alber,Alexandre V. Morozov,Mark Alber,Alexandre V. Morozov,&keyword=&subject=Research Article,Computer and Information Sciences,Network Analysis,Scale-Free Networks,Computer and Information Sciences,Network Analysis,Computer and Information Sciences,Network Analysis,Protein Interaction Networks,Biology and Life Sciences,Biochemistry,Proteomics,Protein Interaction Networks,Engineering and Technology,Control Engineering,Control Theory,Computer and Information Sciences,Systems Science,Control Theory,Physical Sciences,Mathematics,Systems Science,Control Theory,Physical Sciences,Mathematics,Probability Theory,Probability Distribution,Biology and Life Sciences,Evolutionary Biology,Evolutionary Processes,Convergent Evolution,Computer and Information Sciences,Systems Science,Dynamical Systems,Physical Sciences,Mathematics,Systems Science,Dynamical Systems,Biology and Life Sciences,Molecular Biology,Macromolecular Structure Analysis,Protein Structure,Protein Structure Networks,Biology and Life Sciences,Biochemistry,Proteins,Protein Structure,Protein Structure Networks,