Network link prediction by global silencing of indirect correlations.

Predictions of physical and functional links between cellular components are often based on correlations between experimental measurements, such as gene expression. However, correlations are affected by both direct and indirect paths, confounding our ability to identify true pairwise interactions. Here we exploit the fundamental properties of dynamical correlations in networks to develop a method to silence indirect effects. The method receives as input the observed correlations between node pairs and uses a matrix transformation to turn the correlation matrix into a highly discriminative silenced matrix, which enhances only the terms associated with direct causal links. Against empirical data for Escherichia coli regulatory interactions, the method enhanced the discriminative power of the correlations by twofold, yielding >50% predictive improvement over traditional correlation measures and 6% over mutual information. Overall this silencing method will help translate the abundant correlation data into insights about a system's interactions, with applications ranging from link prediction to inferring the dynamical mechanisms governing biological networks.

The currently incomplete maps of molecular interactions between cellular components limit our understanding of the molecular mechanisms behind human disease[1-6]. Ultimately, high-throughput mapping projects[7-10] are expected to provide the accurate maps of interactomes necessary to systematically unlock disease mechanisms. Yet, as a complete interaction map is at least a decade away, we need to develop tools that allow us to infer the structure of cellular networks from empirically obtained biological data[11,12]. Many current tools designed to infer functional and physical interactions in the cell rely on the global response matrix, which captures the change in node i's activity in response to changes in node j's[13]. This matrix can be measured directly from gene knockout or overexpression experiments, or inferred indirectly using related measures such as Pearson or Spearman correlations[14], mutual information[15,16] or Granger causality[17]. Traditional methods for predicting links[15,16,18,19] assume that the magnitude of G correlates with the likelihood of a direct functional or physical link between nodes i and j. Yet G cannot distinguish between direct and indirect relationships: a path i → k → j can result in a significant response measured between i and j, falsely suggesting the existence of a direct link between them ().

Several methods correcting for such effects have been proposed: information theoretic approaches evaluate the association between nodes by measuring the entropy of their mutual activities, where a low entropy indicates a statistical dependence between the node activities[16,18,20]; probabilistic models, such as the graphical Gaussian model, allow one to evaluate the correlation between i and j, while controlling for the state of node k, and thereby provide a more indicative measure of direct linkage[21-25]; other models rely on assumptions pertaining to the network topology, such as the tendency of real networks to exhibit strong degree correlations[26]. The ultimate solution, however, should enable us to fully unwind the direct from the indirect effects, providing a measure which distinctly indicates the existence of direct links. Consequently, here we focus on the local response matrix in which the contribution of indirect effects is eliminated. In contrast with (1), which allows for global changes in i and j's environment, here the “∂” indicates that S is defined to capture only local effects, namely the response of i to changes in j when all surrounding nodes except i and j remain unchanged. Hence S > 0 implies a direct link between i and j.

Here we derive a method for calculating the local response matrix (2) from experimentally accessible correlation measures, allowing us to mathematically discriminate direct from indirect links. We show that the resulting S matrix, in which the contribution of indirect paths is silenced, is more discriminative than the empirically obtained G matrix, enhancing our ability to extract direct links from experimentally collected correlation data.

Results

The silencing method

To extract S from the experimentally accessible G, we formally link (1) and (2) via Equation (3) is exact and the sum accounts for all network paths connecting i and j (). It is of limited use, however, as it requires us to solve N2 coupled algebraic equations. In we show that (3) can be reformulated as where I is the identity matrix and sets the off-diagonal terms of M to zero. To obtain an approximate solution for S we use that fact 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 allows us to approximate with (), obtaining Equation (5), our main result, provides S from the experimentally accessible G. It achieves this through a 'silencing effect’, in which direct response terms are preserved, while indirect responses are silenced. To understand this consider a specific term in G, documenting the response of node i to j's perturbation. As indicated by Eq. (3), this response is a consequence of all direct and indirect paths leading from j to i. As we document below, the transformation (5) detects the indirect paths and silences them, maintaining only the contribution of the direct paths ().

Silencing in model systems

To demonstrate the predictive power of (5), we implemented Michaelis-Menten dynamics on a model network (), as commonly used to model generegulation[27,28]. We obtained G by perturbing the activity of each node and then calculated S using (5). Figure 2a shows the G and S terms associated with interacting (green) and non-interacting (orange) node pairs. Although G is higher for direct interactions, the overlap between the orange and the green symbols indicates a lack of a clear threshold q that separates direct and indirect interactions. In contrast, S displays a clear separation between direct and indirect interactions, accurately predicting each direct link. Indeed, the ROC curve derived from G () has an area of AUROC = 0.91, reflecting inherent limitations in separating direct from indirect interactions based on G only. In contrast for S we obtain AUROC = 0.997 (blue), where the true positive rate (TPR) reaches 100% with a false positive rate (FPR) of less than 10–3. Also, although for G precision increases gradually with the threshold q (), for S precision jumps to one for q > 10–4. Hence, in our well controlled model system effectively any non-zero S corresponds to a direct link.

Figure 2

Network inference in model systems

We numerically simulated Michaelis-Menten dynamics on a scale-free network [40-42], extracting the correlations G between all pairs of nodes (see Sec. S.III for details). (a) G and S associated with interacting (green) and non-interacting (orange) node pairs. S silences the correlations associated with indirect interactions, resulting in a clear separation between direct and indirect interactions, a phenomenon absent from G. (b) ROC curve obtained from G (red, area 0.91) and S (blue, area 0.997). The S network reaches 100% accuracy with a negligible amount of false positives. (c) Precision obtained for threshold q for G (red) and S (blue). The gradual rise of the G-based precision indicates that for a broad range of thresholds only a small fraction of the links will be identified. In contrast, the steep rise in precision for S indicates its enhanced discriminative power between direct and indirect links: virtually any non-zero S corresponds to a directly interacting pair. (d) The discrimination ratio, Δ, is much higher in S (blue) compared to G (red). This indicates that S is a much better predictor of direct vs. indirect interactions. The silencing (5), which captures the increase in the discrimination ratio is κ = 15.0. (e) Silencing increases with the path length d between i and j, so that the more indirect is the link the more dramatic is the silencing. (f) The source of S's success is the silencing effect, here illustrated on correlations measured for a linear cascade. The reconstruction of the cascade from G is confounded by numerous non-vanishing indirect correlations. In S the indirect correlations are silenced, providing a perfect reconstruction.

The performance of (5) is due to the silencing effect: it leaves G unchanged if i and j are linked, while it systematically lowers all G not rooted in a direct interaction. To quantify this effect we measured the discrimination ratio Δ = 〈G〉Dir/〈G〉Indir (Δ = 〈S〉Dir/〈S〉Indir) which captures the ratio between G (S) terms associated with direct links and those associated with indirect links (). We find that S is much more discriminative than G owing to its silencing of indirect responses. To quantify this effect we measure the silencing which captures the increased power of S to discriminate between direct and indirect links compared to G. In our model system we find that κ = 15, a silencing of more than an order of magnitude (). Furthermore, the longer is the distance d between two nodes, the larger is the silencing (). As an illustration, consider a linear cascade in which changes in any node result in a finite response G by all other nodes (). Equation (5) silences all indirect responses, while leaving the response of direct links effectively unchanged, offering a discriminative measure that enables a perfect reconstruction of the original network.

Predicting molecular interactions in E. coli

To test the predictive power of (5) on real data we used the E. coli datasets distributed by the DREAM5 network inference challenge[19]. The input data include a compendium of microarray experiments measuring the expression levels of 4,511 E. coli genes (141 of which are known transcription factors) under 805 different experimental conditions (). We constructed three separate global response matrices G between the 141 transcription factors and their 4,511 potential target genes, based on (i) Pearson correlations; (ii) Spearman rank correlations; and (iii) mutual information, which are three commonly used methods for link detection (). From each of the three G matrices we obtained S via (5), and compared the performance of G with the pertinent S. To validate our predictions we relied on the gold standard used in the DREAM5 challenge, consisting of 2,066 established gene regulatory interactions. Measuring AUROC from G and S, we find an improvement of 56% for Pearson correlations (), 67% for Spearman rank correlations () and a smaller improvement of 6% for mutual information (), e.g. allowing us to improve upon the top performing inference methods[19].

We further tested the discrimination ratio, Δ, and the silencing, κ, for each of these methods, finding that indirect correlations are subject to an average of two-fold silencing in the transition from G to S (). Silencing is especially crucial in the presence of the cascade and co-regulation motifs shown in Figures 3e-f, where most inference methods indicate a spurious link between X and Y owing to the indirect correlation mediated by node I. Indeed, the transformation (5) silences these indirect correlations by a factor of three or more for Pearson and Spearman correlations and by a smaller factor for mutual information, overcoming one of the most common hurdles of inference methods, which tend to over-represent triadic motifs[19].

Figure 3

Inferring regulatory interactions in E. coli

(a) Starting from gene expression data, we used Pearson correlations in expression patterns to construct G for 4,511 E. coli genes, obtaining S via (4). We compared our predictions to a gold standard of experimentally verified genetic regulatory links [19]. The area under the ROC curve (AUROC) is increased from 0.59 to 0.64 in the transition from G to S, representing a 56% improvement (above the baseline of 0.5 for a random guess). (b) An improvement of 67% is observed for Spearman rank correlations. (c) A less dramatic improvement of 6% is shown when G is constructed using mutual information. (d) The discrimination ratio for all three methods compared with that obtained from the pertinent S matrix. The transition to S (4) increases the discrimination between direct and indirect interactions by a factor of two or more, so that indirect interactions have a significantly lower expression in S. (e) - (f) This observation becomes even more dramatic when focusing on two specific motifs: cascades and co-regulators. In G the indirect correlation between X and Y, which is induced by the intermediate node, I, may lead to the false prediction of the spurious X – Y link. Thanks to silencing, the discrimination between the direct and indirect links in these motifs is increased by a factor of three or more for Pearson and Spearman correlations, and by a factor of about two for mutual information.

The effects of noise and uncertainty

As all experimental data is subject to noise, the global response matrix, G, is characterized by some degree of uncertainty. To test the performance of our methodology in the presence of noise, we repeated the numerical experiment of Figure 2, this time adding Gaussian noise to G, which allows us to calculate silencing as a function of increasing the signal to noise ratio θ (). As expected, silencing is unaffected by small values of θ, so that κ features a plateau below θ ≲ 0.1. For large θ, silencing decays as κ ~ θ–1, demonstrating that the performance of the method decreases slowly with increasing the signal to noise ratio. Indeed, as opposed to a rapid exponential decay, the observed slower power-law dependence indicates that the method is rather tolerant against noise. Silencing is lost only when the noise reaches the critical level θC ≈ 0.75, when the signal is almost completely overridden by noise, leading to κ = 1 ().

Hidden nodes offer another source of uncertainty. They represent the fact that in most cases we are unable to read the states of all nodes in the system[29]. To illustrate the effect of the hidden nodes on the performance of the silencing method, we consider the case of a simple cascade i → k → j, where the intermediate node k is hidden. In this scenario, Eq. (5) will not be able to silence the indirect i → j link because in the observable system the G term cannot be attributed to any indirect path. Hence, absent any other information about the system, it is mathematically impossible to infer the indirectness of G, as the removal of k isolated i from j[30]. This touches upon the fundamental mechanism of silencing: as illustrated in Figure 1 (and ) the silencing transformation (5) exploits the flow of information through indirect paths. Consequently, if as a result of hidden nodes the network fragments into several components such that the node pair i and j become isolated from each other, then all indirect paths between them became hidden and the pertinent G term will not be silenced (). Hence silencing is expected to fail only when the network breaks into many isolated components so that most node pairs become isolated. Fortunately, a fundamental property of complex networks is that with average degree 〈k〉 >> 1, one needs to remove a large fraction of the nodes to fragment the underlying giant connected component[31-34]. Therefore we can build on percolation theory, which allows us to analytically predict how the size of the largest connected component changes with the random removal of a certain fraction of nodes[35,36]. The calculation shows that silencing is maintained as long as the fraction of hidden nodes is smaller than where (). This equation indicates that for large 〈k〉 the method will be reliable even if a large fraction of the nodes are hidden.

Figure 1

Silencing indirect links

(a) The experimentally observed global response matrix, G, accounts for direct as well as indirect correlations, with no clear separation between them. The source of G could be gene coexpression data, statistical correlations or genetic perturbation experiments. (b) In the absence of a clear separation in G assigned to direct and indirect correlations, our ability to infer direct physical links (solid lines) is limited. Simple thresholding, i.e. accepting all links for which G exceeds a predefined threshold, is known to predict spurious links (strong dashed lines) and overlook true links (light solid lines). (c) While the average G terms associated with direct links (dark blue) are higher than the average terms associated with indirect links (light blue), as captured by the discrimination ratio, Δ, the difference is not sufficient to identify direct and indirect links. (d) Silencing is achieved through Eq. (5), which exploits the flow of information in the network: the flow from the source (j) to the target (i) is carried through the indirect effect G (orange) coupled with the direct impact S of the target's nearest neighbor κ (blue). By silencing the indirect contributions, Eq. (5) provides the local response matrix, S , whose non-zero elements correspond to direct links. (e) – (f) In S the terms associated with indirect links are silenced, allowing us the detect only the direct links of the underlying network. (g) As indirect terms become much smaller in S, we obtain a greater discrimination ratio, Δ. The degree of silencing, κ, captures the increase observed in the discrimination ratio by the transition from G to S (5).

To test this prediction, we revisited the numerically obtained G analyzed in Figure 2 and measured the degree of silencing after randomly removing an increasing fraction of nodes. In each case we also measured the ratio between isolated and connected node pairs (ρ). We find that, as predicted, the degree of silencing is driven mainly by ρ, approaching κ ≈ 1 (no silencing) when ρ ≥ 1, namely when the isolated pairs begin to dominate the network (). Here as 〈k〉 = 4, Eq. (7) predicts ηC ≈ 0.57, i.e. the method will fail only when almost 60% of the nodes are hidden. Note that for biological networks 〈k〉 is expected to be in the range of[37] 〈k〉 ≲ 10, predicting ηC ≲ 0.8. Namely, one needs to lose access to 80% of the nodes for silencing to lose its effectiveness.

Discussion

With computational complexity , Eq. (5) is scalable and requires no assumptions about the network topology. By silencing indirect effects, it turns the raw correlation data into a predictive S matrix, dominated by direct interactions. It is especially suited to treat perturbation data, such as genetic perturbation experiments, in which case G describes the response of all genes (dx) as a consequence of the perturbation of the source gene (dx)[38]. In practice, however, G could be the result of a broader set of experimental realizations where other measures are used to evaluate the association between nodes, typically statistical measures such as Pearson or Spearman correlation coefficients. Still, our empirical results () clearly show that the transformation (5) successfully applies to these empirically accessible measures as well. Hence silencing is largely insensitive to the specific process by which G was constructed.

The method's broad applicability is rooted in the fact that it does not depend on the value of each specific term in G, but rather on the global relationships between them. Indeed, the global structure of G reflects the patterns of propagation of the perturbations along the network. Equation (5) helps uncover these paths from the raw data, disentangling the direct from the indirect effects. These patterns of information flow are inherent to the underlying network structure, and should not depend on the specific experimental realization of (1). For instance, a cascade i → j → k will be characterized by a decreasing correlation propagating along the arrows, a large correlation between i and j and a weaker one between i and k. Although the magnitude of these correlations might depend on the size or the form of i's perturbation as well as on the statistical measure we used to evaluate them, the decay pattern required to infer the structure of the cascade is an inherent property of the network flow and can be successfully detected by the silencing method ().

The silencing transformation is derived from fundamental mathematical principles of dynamical correlations in networks. Hence it is expected to apply under rather general conditions. However, as Equation (5) indicates it requires that the input matrix, G, is invertible. This imposes some limitations when constructed from statistical correlation measures. For instance in the empirical results of Figure 3a we constructed G from Pearson correlations, using the states of 4,511 nodes measured under 805 experimental conditions. In general, if the number of experimental conditions is smaller than the number of nodes the resulting Pearson correlation matrix may be singular. In this case additional processing will be required before (5) could be applied. Here, following the DREAM5 protocol, we only focused on the correlations between the 141 known transcription factors and the rest of the nodes, which lead to an invertible G (). Other means to ensure G's invertibility are discussed in .

Isolating indirect effects in correlation data, a fundamental challenge of network inference, is typically approached through local probabilistic tools[12,14-18]. In contrast, the success of the silencing method is rooted in its exploitation the global network topology[39]: it relies on the fundamental principles of network structure and dynamics to identify and silence the effects of indirect paths. The ability to extract S from G could also have implications for our understanding of network dynamics. Indeed, G is a global network measure, as its magnitude is determined by the numerous indirect paths connecting i and j. Hence, for a given dynamics, the G matrix will take a different form depending on the network topology, making it a poor predictor of the system's dynamics. By eliminating indirect effects, S measures the effect gene i would have on gene j had they been isolated from the rest of the network. It thus helps us quantify the dynamical mechanism that governs individual pairwise interactions, avoiding the convolution of dynamical and topological effects present in experimental data. For instance, consider a set of perturbation experiments providing G. The structure of G reflects the microscopic mechanisms that govern the pairwise interactions, e.g. genetic regulation and biochemical processes. It is difficult, however to extract this information from G since its terms are a convolution of many interactions, reflecting the many paths leading from i to j. The transition to S , via (5), allows us to treat each isolated interaction on its own, providing a direct observation into the microscopic interaction mechanism. Direct application of this fact could be the derivation of a rate equation that governs the system's dynamics from G, as well as predicting the universality class and the scaling laws governing the system's response to perturbations. Hence (5) helps translate the ever-growing amount of data on global correlations into valuable local information.