To the Editor:

In the August 2013 issue of this journal, Barzel & Barabási reported a method for reconstructing network topologies¹. Here we show that the Barzel & Barabási method is a variant of a previously published method, modular response analysis (MRA)². We also demonstrate that the implementation of their algorithm using statistical similarity measures as a proxy for global network responses to perturbations is erroneous and its performance is overestimated.

The reconstruction of network connections from data remains a fundamental problem in biology. It is not immediately obvious how to capture direct links between individual network nodes from experimental data because a perturbation to a component propagates through a network, causing widespread (global) changes, thereby masking direct (local) connections between nodes. This question has been previously studied in >100 publications, collectively representing MRA (reviewed in refs. 3,4,5,6,7). MRA quantifies direct interactions between network nodes (i and j) using the local response coefficients (also known as connection coefficients), which describe direct effects of a small change in node j on node i, while keeping the remaining nodes unchanged to prevent the spread of the perturbation. The local responses cannot be directly assessed, whereas the global responses can be measured; when following a perturbation to node j, the entire network relaxes to a new steady state. In this new state, nodes that are directly affected by node j have changed, but other indirectly affected nodes also change because the initial perturbation has propagated through the network. MRA has solved the problem of finding local, direct links between components through the global responses for networks of any size and complexity². The development and application of methods that are conceptually similar to MRA (e.g., regulatory strength analysis⁸ and maximum likelihood-based MRA^9,10) has reinforced the validity of using MRA-type methods to reconstruct network connections^{11,12,13,14,15,16}.

The Barzel & Barabási study¹ uses the same concept and strikingly similar terminologies to reconstruct networks by deriving the local connection coefficients from the global response coefficients. Key equations (3) and (4) in their silencing method¹ express the local coefficients in terms of the global coefficients and are a subset of the published MRA equations^{2,9,10,17,18,19} with a formal replacement of the diagonal elements of the local response matrix by zeros instead of minus ones (Supplementary Note 1). Another formal difference is that the variant of the global response matrix used by Barzel & Barabási¹ considers the global change in each node that results from the change in every other node of the network, whereas MRA more generally considers the global changes in network nodes that result from changes in any parameter that might affect several nodes simultaneously^18,19. Barzel & Barabási claim as one of the main outcomes of their study an approximate solution to equations (3) and (4) (equation (5))¹, whereas MRA offers an exact solution^2,18,19. Both the approximate¹ and the exact MRA² solutions require the inversion of the global response matrix. Consequently, Barzel & Barabási's approximation¹ does not decrease the computational complexity of the exact solution^2,18,19. In fact, the proposed approximate iterative method (equations S12 and S13)¹ needs repetitive matrix calculations, making it slower than the existing MRA algorithms that provide the exact solution.

The Barzel & Barabási approximate solution relies on the assumption that “typically, perturbations decay rapidly as they propagate through the network, so that the response observed between two nodes is dominated by the shortest path between them”¹. This assumption disregards well-documented biological evidence, such as the sensitivity amplification in signaling cascades^20,21 and the existence of positive feedback loops in biological networks, which amplify initial signals as they propagate through a network²². The common occurrence of positive or double-negative feedback loops invalidates the Barzel & Barabási assumption for many known signaling pathways, including the restriction point pathway²³, cell cycle signaling^24,25,26, mitogen-activated protein kinase (MAPK) cascades^27,28 that are evolutionary conserved from yeast to mammals, as well as transcription regulation networks²⁹. In these regulatory networks, global responses of the neighboring nodes outside of positive feedback loops are typically smaller than the response of a node that lies inside a positive feedback loop to an upstream node outside the loop (Supplementary Note 1, Supplementary Note 2 and Supplementary Fig. 1). Sensing and processing of stimuli is the normal function of most, if not all, regulatory biological pathways, whose common feature is to increase the response of a 'target' node to a 'source' node with the distance between them³⁰. For example, in an experimental study, Bastiaens and colleagues³¹ used MRA to unravel the direct linkage topology and strengths of connections in the MAPK cascade in PC12 cells that were stimulated with epidermal growth factor (EGF) versus nerve growth factor (NGF). They calculated the local response coefficients from experimentally measured global responses and found that EGF elicits negative feedback, whereas NGF induces positive feedback, imposed on the backbone of the same MAPK pathway that propagates signals from both growth factors. The experimentally measured response between immediate neighbors, such as the global response of MEK to small interfering RNA (siRNA) against RAF, was much smaller than the response of the more distant neighbor ERK to RAF siRNA under both EGF and NGF stimulation, regardless of the growth factor–specific difference in network wiring³¹. In contrast to Barzel & Barabási's assumption, real biological networks do not feature a rapid decay of specific perturbations because these pathways have evolved to sense and respond to external cues by processing, amplifying and integrating the signals. Thus, reconstruction of these pathways using the 'average perturbation decay' hypothesis misses key functional features of biological pathways.

Network reconstruction methods that exploit the global network responses require systematic perturbation measurements. Clearly, many high-throughput approaches do not provide such perturbation data, whereas statistical similarity measures can be calculated from omics data. In the absence of perturbation data, Barzel & Barabási reconstructed the network by substitution of the global response coefficients with statistical similarity measures, such as the Pearson and Spearman correlation coefficients, and mutual information¹. However, this substitution yields a symmetric correlation matrix with very different mathematical properties from the global response matrix. Inference based on these symmetric measures (using equation (5) and equations S11–S13 of ref. 1) results in networks where the local response matrix is always symmetrizable (Supplementary Note 3). In these inferred networks, the absence of a direct connection from node j to node i (S_ij = 0) inevitably implies that there is no direct connection in the reverse direction (node i to node j, S_ji = 0). Thus, the inferred networks cannot have a one-way connection between any two nodes. This approach violates the reality of cellular networks where one-way connections dominate. For example, ubiquitous post-translational protein modifications are typically one-way connections. When a kinase phosphorylates a substrate, the substrate usually does not phosphorylate the kinase. Equally important, the symmetrizable local response matrix implies that the overall signal amplification or attenuation along a circular path (for instance, formed by a feedback loop) is exactly the same as in the reverse direction along this path, which is also biologically unrealistic. Therefore, the local response coefficients inferred from the correlation or mutual information matrices¹ instead of the global response matrices do not represent real cellular networks and predict erroneous network connections (Supplementary Note 3).

There are many established methods that use statistical similarity measures to identify connections between network nodes, including lasso regression-based methods^32,33,34, the partial correlation method^35,36, Gaussian graphical models³⁷ and elastic nets³⁸, which rank the predicted edges on the basis of correlations between experimental observations or regression coefficients, rather than incorrectly using MRA equations to express the local connections through the correlation or mutual information matrices.

Barzel & Barabási claim that their method is robust against noise¹. However, neither the Barzel & Barabási algorithm nor the standard MRA² take explicit precautions against extrinsic or intrinsic noise in the data (Supplementary Note 4). Several statistical reformulations of MRA have been developed to allow robust inference of network topology in the presence of noise. For example, the statistical adaptations of MRA based on the Monte Carlo³¹ and maximum likelihood^9,10 methods were successfully applied to infer signaling pathway topologies from noisy perturbation data in mammalian cells. A recent Bayesian MRA reformulation is also capable of inferring networks from noisy and even incomplete data sets³⁹.

We are also unconvinced by the claim that the Barzel & Barabási method “improves upon the top-performing inference methods”¹ in the DREAM5 network inference challenge (Network 3, http://wiki.c2b2.columbia.edu/dream/data/scripts/DREAM5/files/DREAM5_NetworkInference_Evaluation.zip). The DREAM5 data set contains the expression levels of 4,511 Escherichia coli genes, including 334 known transcription factors, and contestants were asked to rank the likelihood of interactions between transcription factors and target genes. Performances of the contenders were then estimated against a gold standard network, which involved only 141 out of 334 transcription factors. After we were unable to reproduce the results of Barzel & Barabási, the authors provided us with the source code for their algorithm. We found three issues in their analysis that rendered the comparison of their algorithm performance with the performances of the DREAM5 contestants invalid (Supplementary Note 5). First, although the identity of the 141 transcription factors of the gold standard network on which the performance was evaluated was unknown to the contestants, Barzel & Barabási exploited this information to zero out the correlations involving the other 193 transcription factors in their correlation matrix (used as a proxy for matrix G). Second, Barzel & Barabási incorrectly calculated the receiver operating characteristics (ROC) curves that estimated the performance of their algorithm (Supplementary Note 5). The DREAM5 challenge gold standard network contains only 141 transcription factors and 1,080 out of 4,511 potential target genes, whereas the interactions involving the remaining 3,431 genes are neither included in the gold standard network nor considered for the performance evaluation in the DREAM5 challenge (http://wiki.c2b2.columbia.edu/dream/data/scripts/DREAM5/files/DREAM5_NetworkInference_Evaluation.zip). Even so, Barzel & Barabási erroneously count these omitted interactions (between the 141 transcription factors and the remaining 3,431 genes) as true negatives, resulting in the inflated area under the ROC curve (AUROC) estimate (Supplementary Note 5). Finally, to evaluate performance of their algorithm, Barzel & Barabási disregarded the precision recall score (AUPR) and used only the AUROC score, which is known to be insufficient and can be misleading when the numbers of true positive (the interactions present in the gold standard network) and true negatives (the interactions absent in the gold standard network) differ significantly⁴⁰. In the E. coli network, true negatives are about 100-fold more abundant than true positives. Consequently in the DREAM challenge, the performances were estimated using scores that combined both AUROC and AUPR. After we corrected these errors (Supplementary Note 5) and properly recalculated AUROC and AUPR using the evaluation script of the DREAM5 challenge (http://wiki.c2b2.columbia.edu/dream/data/scripts/DREAM5/files/DREAM5_NetworkInference_Evaluation.zip), the Barzel & Barabási inference algorithm performed poorly compared with the best performers in the DREAM5 competition (using either AUROC or AUPR as criteria; Supplementary Note 5; supplementary materials are also available on GitHub and figshare: http://figshare.com/articles/NBT_correspondence/1356170 GitHub: https://github.com/SBIUCD/NBT_correspondence.git). Specifically, the performance of the Barzel & Barabási algorithm estimated by AUPR ranks between 20th and 28th, and using AUROC it ranks between 3rd and 14th of the 29 participants, depending on whether the algorithm was applied to the Spearman correlation, Pearson correlation or Mutual Information. Predictions, in which we used merely 'raw' statistical similarity measures as substitutes of the local connections, ranked better (AUPR ranked between 7th and 13th, and AUROC ranks between 2nd and 10th depending on the statistical similarity measure used) than their algorithm (Fig. 1).

As in the original DREAM5 challenge, the performance was estimated using two scores: (a–c) AUROC. (d–f) AUPR. TPR, true-positive rate; FPR, false-positive rate; 'Precision' indicates the fraction of correctly inferred true interactions; 'Recall' equals TPR. Insets show the AUROC and AUPR scores for the Barzel & Barabási algorithm, the Pearson and Spearman correlations, and mutual information. In all three cases, the performance of the Barzel & Barabási algorithm was lower than the performance of the raw similarity measures alone (note that in Fig. 3 of the original publication¹, the Barzel & Barabási algorithm was inappropriately applied and scored).

Full size image

In summary, the concerns raised in this Correspondence cast doubt both on the level of conceptual advance and the practical usefulness of the silencing method proposed in the Barzel & Barabási study¹. Ironically, many of the practical issues could have been remedied by consulting the extensive published literature describing MRA and MRA-based methods, which Barzel & Barabási unfortunately seem to have overlooked.