Gene regulatory networks control the cellular phenotype by changing the RNA and protein composition. Despite its importance, the gene regulatory network in higher organisms is only partly mapped out. Here, we investigate the potential of reverse engineering methods to unravel the structure of these networks. Particularly, we focus on modular response analysis (MRA), a method that can disentangle networks from perturbation data. We benchmark a version of MRA that was previously successfully applied to reconstruct a signalling-driven genetic network, termed MLMSMRA, to test cases mimicking various aspects of gene regulatory networks. We then investigate the performance in comparison with other MRA realisations and related methods. The benchmark shows that MRA has the potential to predict functional interactions, but also shows that successful application of MRA is restricted to small sparse networks and to data with a low signal-to-noise ratio.

Introduction

Inferring regulatory relationships between genes or signalling molecules from perturbation data is one of the big challenges in computational biology. The major problem in disentangling the network structure is that a perturbation of a single node results in direct and indirect effects on other nodes. For example a transcription factor knockdown results in expression changes in direct targets among which other transcriptional regulators are encoded that propagate the signal further.

To discriminate cause from correlation it is vital to generate data that allow the re-construction of directed regulatory networks. Two data collection strategies have shown this potential: (i) time series data and (ii) pre- and post-perturbation steady state measurements.

Reverse engineering from time series data has been successfully applied by the use of, e.g. ordinary differential equations (ODE) [1,2] or dynamic Bayesian networks [3]. Here, causality of regulation is inferred by the time delay between the expression of regulator and targets which requires an appropriate sampling rate. Since the optimal sampling rate is often not known beforehand, a substantial number of time points have to be measured for each perturbation, rendering this approach experimentally costly to reconstruct larger networks.

In contrast, steady state measurement-based methods need a larger number of perturbations to infer causality. However, since for each perturbation only the steady state is measured, smaller overall datasets are required. Common modelling approaches for this data type are utilising Boolean or Bayesian models [4,5], or models that characterise the system using dynamical systems theory [6,7], which is particularly suited to incorporate feedbacks. Most often, the dynamics of a network is modelled using ODEs, where the change in the nodes xi is given by:

 
formula
(1)

Most approaches then apply a linearisation of the underlying dynamical system around the (unperturbed) steady state value [2,8]:

 
formula
(2)

with Δx representing small deviations in the system variables around the steady state. Jij are the elements of the Jacobian matrix J and represent the direct influence from variable j to i. These entries are zero for no direct connection and, if non-zero, the sign indicates the type of influence, where negative and positive values represent inhibition and activation respectively.

One way to model a small, persisting perturbation is by linear expansion of the effect of a parameter change Δp in (eqn 1), which can then be solved for steady state:

 
formula
(3)

This is the central equation that many reverse engineering methods aim to solve. If the perturbation Δp affects only node i, the entries of vector ∂fj/∂p are zero except for the entry of the perturbed node. Under most settings, the identity of the perturbed node is known but often the initial perturbation strength cannot be quantified. The entries of the vector (Δxj)ss correspond to the measured changes in steady state after perturbation. For multiple perturbation experiments, the equations can be arranged in matrix form. Then, (Δx)ss represents a matrix where all measurements for one perturbation are stored column-wise.

Two distinct approaches apply (eqn 3) to infer networks from steady state perturbation data [9]. The first, called Network Identification by multiple Regression (NIR) [6] formulates this task as multiple regression problems for rows of J under the constraint that these rows are sparse. The second approach is based on estimating the Jacobian matrix in its entirety. This method has been independently proposed in two flavours: one flavour, called Regulatory Strength Analysis (RSA), aims at estimating the Jacobian itself [10]. The other, called Modular Response Analysis (MRA), estimates the dimensionless, scaled and thus more robust version of the Jacobian matrix [11–13]. A disadvantage of NIR is that it is not parameter free and requires knowledge of the perturbation strength, requiring additional experiments.

MRA is an approach that uses concepts of metabolic control analysis (MCA) [14,15] and defines two types of response coefficients [12]: (i) global response coefficients that represent normalised steady-state changes and (ii) local response coefficients representing entries of the normalised version of the Jacobian matrix [12,16].

Entries of matrix R can be experimentally measured, whereas the entries of r represent the ‘strength’ of direct interaction between nodes in the network. The relationship of the two response coefficients is described in the central equation of MRA. Assuming the diagonal matrix represents the scaled quantification of direct single node perturbations (with ) (eqn 3) can be reformulated as:

 
formula
(4)

If each node has been perturbed, the quantifications of can be calculated from R [13], allowing to estimate the direct interactions r from the measurable global response matrix alone.

The workflow of MRA to derive r is exemplified in Figure 1. In the illustrated 4-node example it is shown how a sparsely connected network with feedback can still lead to a full global response matrix, rendering an intuitive derivation impossible.

The basic principle of Modular Response Analysis

Figure 1
The basic principle of Modular Response Analysis

Workflow of MRA-mediated network reconstruction from systematic perturbation data (P=-50%). As data input MRA requires the steady state data of all nodes before (x(0)) and after single node perturbation (x(1)). This information is then combined to derive the global response matrix R. In here R is a full matrix, stating that every perturbation reaches all other nodes. By applying the central equation of MRA to R, the local response matrix r can be determined, which contains the network structure.

Figure 1
The basic principle of Modular Response Analysis

Workflow of MRA-mediated network reconstruction from systematic perturbation data (P=-50%). As data input MRA requires the steady state data of all nodes before (x(0)) and after single node perturbation (x(1)). This information is then combined to derive the global response matrix R. In here R is a full matrix, stating that every perturbation reaches all other nodes. By applying the central equation of MRA to R, the local response matrix r can be determined, which contains the network structure.

The resulting local response coefficients can be interpreted such that they represent the relative influence that the incoming node has on the outgoing node. A local response coefficient therefore means that a change in node j leads to the same (relative) change in node i. In this manner local response coefficients larger than one represent amplification and smaller than one attenuation [16]. The notion of nodes is not confined to represent single reaction molecules but modules that are defined by incoming and outgoing signal as long as linear approximations can be justified [16]. Hence the name Modular Response Analysis.

Due to the specifics of MRA it is a potent methodology to assist in reverse engineering of biological networks in general and genetic networks in particular for two reasons. First, the currently most widespread genetic perturbation techniques, siRNA and RNA overexpression (and its CRISPR counterparts) are fulfilling the requirement of steady state measurements as they are typically measured at suitable times to approach (quasi) steady states. And second, MRA is able to resolve the expected frequent occurrences of feedback and feedforward structures.

It should be noted, however, that the MRA approach in its original definition is not applicable to experimental data for various reasons:

  1. Experimental data are noisy and thus the global response matrix has to be derived from noisy data. A noisy R will propagate the noise to the local response matrix, where all entries will be different from zero pretending a fully connected network regardless of the real situation. In addition, through inversion of R noise might be amplified unequally.

  2. MRA assumes a linearised system, however, many biological reactions are not strictly linear, which is aggravated when modularisation is introduced. Although the response to small perturbations will be linear, often larger perturbations are required to obtain a sufficient signal-to-noise ratio.

  3. The steady state assumption itself proves to be another critical point in biological systems [13].

To enable the application of the promising approach of MRA to real experimental datasets, adjustments are required that enable MRA to reliably distinguish between real regulatory interactions and noise. The easiest solution would be to set an arbitrary threshold for the ’s to distinguish between noise and regulation, as has been done in [9] for the application of the MRA variant RSA. However, a general threshold will fail to encompass the varying specificity of the detection molecules, e.g. probe specificity. Furthermore, linear noise thresholding can not account for the non-linear propagation of noise terms that occur via inversion of R. Therefore an alternative procedure needs to be developed to distinguish between real and chance connections in r.

To tackle this problem, we developed an algorithm that estimates local response coefficients using a maximum likelihood estimation in combination with a greedy hill climbing model selection approach, termed MLMSMRA, and successfully applied it to reverse engineer a network from a combination of genetic and signalling perturbations [17]. In this study we benchmark the reverse engineering capacity for sensitivity to the obstacles of perturbation strength, noise, non-linearity, network size and sparsity and compare it with alternative approaches.

Derivation of MLMSMRA

By regarding the theoretical global response matrix as a function of r and the perturbation , the main (eqn 4) can be reformulated as:

 
formula
(5)

The goal is to minimise the difference between the theoretical () and the experimentally derived global response matrix (R). Therefore an objective function was developed that gives the negative logarithm of the likelihood:

 
formula
(6)

Here, is the error of the corresponding global response coefficient that can be estimated from replicate experiments by an appropriate error propagation model (see Supplementary Section A8). The perturbation strength can be estimated by setting the derivative of (eqn 6) to zero.

 
formula
(7)

Solving the thusly found global minimum for results in:

 
formula
(8)

Thus for any given r, the corresponding perturbation strength vector can be calculated by weighted least squares regression.

To find coefficients that match the observed global response matrix best in a maximum likelihood sense, r (and indirectly ) have to be selected such that they minimise the weighted sum of squared deviations of (eqn 6), which can be interpreted as χ2-values.

Since it is of interest to determine the optimal number of non-zero entries of r, one would ideally have to test all possible combinations, which is computationally feasible only for very small networks. To circumvent the curse of dimensionality a greedy hill-climbing model selection strategy was devised that is illustrated in Figure 2. At the start we assume a completely unconnected network, i.e. r = −I. Subsequently for each link it is systematically tested whether addition significantly (α = 0.05) improves the model fit, as judged by a likelihood ratio test (cf. [18]). The most significant link is added and this procedure is repeated until no further link significantly improves the fit. In the course of the extension procedure, redundant paths might be added which would not have been significant when added in a different order. Thus after completion of the extension phase the network is iteratively pruned of redundant links (i.e. insignificant reduction in log-likelihood upon removal) to arrive at the final network structure.

Model selection by MLMSMRA

Figure 2
Model selection by MLMSMRA

The procedure initialises the local response matrix r by a negative unity matrix. During the extension phase every possible link is added separately and parameters are estimated by a maximum likelihood approach. For the link that improves the likelihood most, it is tested whether this improvement is significant or just due to increased model complexity. If significant, the link is kept and the prediction will be taken as the new standard to compare with when trying to add another link. This procedure is continued until no significant improvement can be found by any added link. Next, in the reduction phase one link at a time is removed and those links are kept where the prediction would significantly decrease. After no further redundancies are present the final network is reached.

Figure 2
Model selection by MLMSMRA

The procedure initialises the local response matrix r by a negative unity matrix. During the extension phase every possible link is added separately and parameters are estimated by a maximum likelihood approach. For the link that improves the likelihood most, it is tested whether this improvement is significant or just due to increased model complexity. If significant, the link is kept and the prediction will be taken as the new standard to compare with when trying to add another link. This procedure is continued until no significant improvement can be found by any added link. Next, in the reduction phase one link at a time is removed and those links are kept where the prediction would significantly decrease. After no further redundancies are present the final network is reached.

Because the method combines MRA, maximum likelihood (ML) parameter estimation and model selection (MS), the resulting algorithm is termed as MLMSMRA (implementation details in Supplementary Section A4). Within this framework, introduction of previous knowledge is straightforward: if links are known, they can be introduced in the starting network, and if links are biologically not reasonable, the corresponding entry in r can be fixed to zero.

Benchmark on in silico networks

To investigate the potential of MLMSMRA to reconstruct genetic networks from perturbation data, artificial genetic networks were generated and evaluated (see Supplementary Section A1). The test sets were chosen such that they most closely resemble real experiments in structure, dynamics and parameter distributions. The networks consist of mRNA and their corresponding proteins which can interact as transcription factors or as post-translational modifiers. Systematic perturbations were only acting on mRNA level, mimicking classical siRNA or overexpression experiments. Thus not all nodes are directly perturbed but all signals are measured.

The artificial networks are now subjected to be reverse engineered from scratch by MLMSMRA. To investigate the optimal experimental settings, parameters were varied for external (perturbation strength, network size) as well as internal (noise, non-linearity, sparsity) characteristics. The respective default network consists of a ten-species network with 20% link density, mild non-linearity (hill coefficient = 2.5 ± 0.7), perturbation strength of −50% and multiplicative noise averaging 20%, marked by cyan arrows in the following figures.

The predictions were evaluated by interpreting three statistical measures: (i) Matthew’s correlation coefficient (mcc), (ii) sensitivity and (iii) precision (see Supplementary Section A2). Matthew’s correlation coefficient reflects the overall performance and was chosen due to its reported superiority among single-number classifiers [19]. In addition, the interpretation of mcc is straightforward. In analogy to correlation coefficients 1 indicates perfect reproduction, 0 no better than random and −1 stands for absolute disagreement. Thus by mcc it can be directly assessed whether the prediction performs better than random guessing. Sensitivity describes the depth of link recovery and precision defines the validity of all predicted links.

With protein–RNA (pR) and protein–protein (pp) regulation, two qualitatively different network interactions are included with the latter interaction type chosen to monitor MRA capability to predict indirect interactions. To determine differences in predictability, networks containing only one regulation type were tested in addition to the combined scenario. The results for all cases are shown in Figure 3. Due to the similar trend in the three network types, afterwards only transcription factor networks (pR) as the most important subgroup will be depicted. However, in Supplementary Section A3 boxplots for all scenarios are provided for all parameter benchmarks as well as performance on sole RNA data.

Absolute perturbation strength but not direction effects MLMSMRA performance

Figure 3
Absolute perturbation strength but not direction effects MLMSMRA performance

Evaluation of MLMSMRA prediction on 100 simulated ten-node networks for different perturbation strengths. The initial perturbation, indicated as percent change of the transcription rate, includes an approximated knockout (KO) and overexpression (+100) scheme as well as three knockdown scenarios. Method performance is depicted as medians with upper and lower quartiles for mcc, sensitivity and precision. Three interaction scenarios with the same number of links (18) were investigated: pR (protein–RNA), pp (protein–protein) and pR+pp. The arrow indicates default setting.

Figure 3
Absolute perturbation strength but not direction effects MLMSMRA performance

Evaluation of MLMSMRA prediction on 100 simulated ten-node networks for different perturbation strengths. The initial perturbation, indicated as percent change of the transcription rate, includes an approximated knockout (KO) and overexpression (+100) scheme as well as three knockdown scenarios. Method performance is depicted as medians with upper and lower quartiles for mcc, sensitivity and precision. Three interaction scenarios with the same number of links (18) were investigated: pR (protein–RNA), pp (protein–protein) and pR+pp. The arrow indicates default setting.

Perturbation strength

One of the underlying assumptions of MRA is that the ODE system can be well approximated by a linear system, which is most accurate if perturbations are small. However, if the perturbations are small, the magnitude of the response may be smaller than the noise fluctuations. Therefore, a trade-off between perturbation strength and noise is to be expected. To investigate this systematically, in silico datasets with different perturbation strengths were generated, i.e. changing the RNA production rate while leaving noise constant.

Figure 3 shows that stronger perturbations lead to better performance of the algorithm, such that knockout conditions yield the best predictions, whereas 10% knockdown resulted in predictions that were hardly distinguishable from random. It is noteworthy that not perturbation direction, i.e. knockdown or overexpression, but perturbation strength seems to determine the performance.

When we compared the benchmark for the different network types, we noted that the performance on sole pR and pp networks was largely comparable, indicating that MLMSMRA can predict interactions of nodes that were not directly perturbed (i.e. pp). The slightly lower performance on networks containing both interaction types (pR+pp) is possibly due to the increased number of possible links (cf. network size influence in Figure 5A). To conclude, although MRA relies on linearisation, strong perturbations enable the algorithm to detect more interactions with a higher precision.

Noise and non-linearity

We reasoned that noise in experimental perturbation studies will play a critical role, with more noise leading to worse predictions. Indeed, when changing the noise level and keeping perturbation constant at -50% it can clearly be seen from Figure 4A that the overall performance is declining. This decrease in mcc is due to the decreased median sensitivity (from 77% for no noise to 6% for 50% noise), while the median precision is much less changed (from 47 to 29%). This indicates that the quality of the prediction is comparatively robust to noise but comes at the cost of less recovered links. Thus the best performance is achieved when noise is minimised and perturbations are maximised.

Noise but not co-operativity limit performance

Figure 4
Noise but not co-operativity limit performance

(A) Median statistical performance in response to different noise levels exemplified on transcription factor networks. (B) Boxplots comparing the predictability with 20% noise added to the measurements (Xn20) or added directly to the global response coefficients derived from unnoised measurements (Rn20). (C) Median statistical performance for different co-operativity coefficients h as a proxy for varying non-linearity. The x-axis displays the average h with all σh = 0.71. The asterisk denotes the linear case with fixed h = 1. The arrow indicates default setting.

Figure 4
Noise but not co-operativity limit performance

(A) Median statistical performance in response to different noise levels exemplified on transcription factor networks. (B) Boxplots comparing the predictability with 20% noise added to the measurements (Xn20) or added directly to the global response coefficients derived from unnoised measurements (Rn20). (C) Median statistical performance for different co-operativity coefficients h as a proxy for varying non-linearity. The x-axis displays the average h with all σh = 0.71. The asterisk denotes the linear case with fixed h = 1. The arrow indicates default setting.

Of concern is the low precision hinting to a high percentage of wrongly predicted links, which might be due to the noise propagation. To investigate this, two cases were compared: once with the simulation noised by 20% and once where the global response coefficients were calculated from the unnoised steady states and appropriate noise was subsequently added to the global response matrix (Figure 4B). It can be noted that the predictions on the noised R (Rn20) show a much better sensitivity and more importantly, a much improved precision (from 36 to 72%). We interpret this such that it is essential to have a highly reliable estimate of the steady state before perturbation, and therefore to have sufficient experimental replicates for the reference state.

The basis of MRA is a linearisation of the underlying ODE system around its unperturbed steady state. If the system is non-linear, this approximation might not apply. In order to investigate the effect of non-linearity, the co-operativity coefficient h used to model the protein interaction similar to a hill coefficient (see Supplementary Section A1) was varied (Figure 4C). Neither increasing μh to 4.5 nor assuming linear interactions (constant h = 1) could strongly affect the performance.

Network size and sparsity

While for a long time it was infeasible to generate sufficient data for networks with more than a dozen nodes, experimental advances such as pooled CRISPR screening with single-cell transcriptome read-out [20] allowed to perturb hundreds of transcription factors. Therefore it is of importance to investigate how network size influences performance. Figure 5A shows that the proposed method performs better on smaller networks for both quantity as well as quality. For networks with 20 nodes the sensitivity approaches zero (4%) and the precision drops to 23%. We attribute this strong size dependency to the local optimisation strategy of the link selection and therefore methodological advances are required for larger networks.

Small and sparse networks yield better predictions

Figure 5
Small and sparse networks yield better predictions

Median performance for (A) networks of different sizes and (B) networks exhibiting differing percentage of link densities. (C) Comparison of the median number of real non-zero links (black) and the median number of predicted links (red) for the networks in B. The arrow indicates default setting.

Figure 5
Small and sparse networks yield better predictions

Median performance for (A) networks of different sizes and (B) networks exhibiting differing percentage of link densities. (C) Comparison of the median number of real non-zero links (black) and the median number of predicted links (red) for the networks in B. The arrow indicates default setting.

Finally, we investigated how network sparsity influences the performance of MLMSMRA. Figure 5B shows that sparse networks lead to higher sensitivity, however precision remains largely invariant especially for networks with densities larger than 10%. By comparing the real average number of non-zero links with the average number of predicted links it can be noted that the latter saturates for higher link densities (Figure 5C).

Comparison with alternative methodologies

After the general characterisation of MLMSMRA, we compared the performance with other established reverse engineering methods that calculate quantitative networks from steady state data, namely: NIR [21,6], MRA with thresholding of the response coefficients (termed as ‘plain MRA’), an MRA variant connected to Monte-Carlo (MC) sampling [11,7], in here referred to as MCMRA, and a Bayesian Variable Selection Approach (BVSA) [22] combined with sign information of the original MRA (see methodology description and application details in Supplementary Section A5).

First, we investigated how the algorithms performed for different interaction types on the in silico networks with default parameter settings. For better comparability and lack of consistent cut-off measures the evaluated network connections were limited to the number detected by MLMSMRA, i.e. all networks for the same test set contained the same number of recovered links which were selected after the methods’ inherent scoring system. Since plain MRA, NIR and BVSA were incapable to incorporate measured but not perturbed nodes (in here proteins) all methods were initially tested on RNA only data (bars without outline in Figure 6). Afterwards MCMRA as well as MLMSMRA that are capable to model unperturbed data were provided with RNA and protein data (bars with black outline).

Comparison of steady state reverse engineering methods for different interaction types

Figure 6
Comparison of steady state reverse engineering methods for different interaction types

Average performance on 100 in silico networks with default parameters (error bars=SD). Three different scenarios are depicted: pR, transcriptional regulation, pp, post-translational regulation and pR+pp, a mix of both. Predictions are based on RNA data (bars without outline) or RNA+protein data (bars with black outline). Asterisks denote significant differences of MLMSMRA estimated by a two-sided U-test (P≤0.05) with bracket ends indicating the compared entities.

Figure 6
Comparison of steady state reverse engineering methods for different interaction types

Average performance on 100 in silico networks with default parameters (error bars=SD). Three different scenarios are depicted: pR, transcriptional regulation, pp, post-translational regulation and pR+pp, a mix of both. Predictions are based on RNA data (bars without outline) or RNA+protein data (bars with black outline). Asterisks denote significant differences of MLMSMRA estimated by a two-sided U-test (P≤0.05) with bracket ends indicating the compared entities.

When networks are predicted from RNA data, MLMSMRA significantly outcompetes plain MRA and NIR. Interestingly, MCMRA and BVSA are significantly worse in the mixed scenario but perform similarly well on transcription factor networks (pR). In the test sets comprising protein–protein interaction networks (pp) no method is capable to predict post-translational regulation from RNA data alone. Only when protein information is added a prediction better than random can be achieved (black outlined bars).

When RNA and protein data are provided, MLMSMRA prediction seems to be always slightly better than MCMRA, albeit only for post-translational networks this advantage is significant.

Overall it can be stated that MLMSMRA has improved MRA to be more versatile for different interaction types and seems to be the better choice when indirect interactions (in here post-translational regulations) have to be inferred.

For transcriptional regulation, however, MCMRA performs equally well in the default test set. As this interaction type is the most basal in genetic network modeling a closer inspection on parameter dependencies for those two methods was conducted (see Supplementary Section A6). In there it was evident that MLMSMRA can reproduce TF networks better on RNA+protein data whereas if only RNA data were given MLMSMRA produced more reliable networks from non-linear, small and sparse networks whereas MCMRA fared better on scenarios with low noise and strong perturbations.

External benchmark: DREAM4 network challenge

As the artificial networks tested and evaluated in here present only one of many attempts to simulate and assess genetic regulatory networks, a bias towards MLMSMRA might exist. Therefore another analysis on externally created genetic networks with a different statistical evaluation scheme was conducted.

The Dialogue for Reverse Engineering Assessments and Methods (DREAM) initiative, round 4, challenge 2 [23–25] represents an interesting framework to test predictions of transcription factor networks from RNA data. The challenge encompasses the prediction of five ten-node networks from five different datasets: steady state data for (i) wild type, (ii) knockdown, (iii) knockout and (iv) multifactorial perturbations (all genes perturbed) and (v) five time series of unspecified perturbations affecting a third of the nodes. The networks were modelled by stochastic differential equations (SDE) with structure and parameters selected from existing networks in Escherichia coli and Saccharomyces cerevisiae. As part of the challenge all networks were selected to incorporate feedback structures. The task was to predict the directed unsigned topology by providing a ranking for all 90 possible interactions. This ranking was used to create Receiver Operation Curves (ROC) and Precision-Recall (PR) curves which were scored by weighting the respective area under the curves.

Table 1 depicts the results for ranking the links by BVSA (posterior probability), MLMSMRA (order of link addition), MCMRA (confidence) and plain MRA (absolute strength of ). When predicting from the knockout dataset, all methods would be ranked as third among the participants with BVSA achieving a slightly higher score than MCMRA followed by plain MRA and MLMSMRA, suggesting that MRA-based predictions from steady states represent clearly a valid prediction method. This is particularly remarkable, as the best performing method integrated all five available datasets, including time series [26], whereas the MRA predictions were based on only a subset of those.

Table 1
DREAM 4 challenge 2 results for MRA variants
MethodScoreRank
Team 543 7.13 1 
Team 549 5.29 2 
BVSA KO 4.95 
MC KO 4.79 
plain KO 4.28 
MLMS KO 4.17 
Team 532 3.97 3 
… … … 
Team 459 2.63 15 
BVSA KO 2.62 16 
MC KD 2.16 20 
MLMS KD 1.24 25 
plain KD 0.85 27 
Team 540 0.5 29 
MethodScoreRank
Team 543 7.13 1 
Team 549 5.29 2 
BVSA KO 4.95 
MC KO 4.79 
plain KO 4.28 
MLMS KO 4.17 
Team 532 3.97 3 
… … … 
Team 459 2.63 15 
BVSA KO 2.62 16 
MC KD 2.16 20 
MLMS KD 1.24 25 
plain KD 0.85 27 
Team 540 0.5 29 

Overall score and ranking among the 29 DREAM challenge participants (italic entries) of the four MRA methods BVSA, MLMSMRA, MCMRA and plain MRA applied on knockout data (KO) or knockdown data (KD). The score generated by the evaluation software is based on the AUROC and AUPR curves for the prediction of five networks (cf. Supplementary Section A7). Note that the challenge provided KO, KD, time series and multifactorial perturbation data, of which only the first two could be used by MRA approaches.

The MRA approaches performed better with knockout data than with knockdown data, which re-emphasises our previous findings that small perturbations produce worse predictions than large perturbations.

As the overall score was a combination of scores from five networks, we also evaluated the performance on individual networks between MCMRA, plain MRA and MLMSMRA (Supplementary Section A7). Here, we found that in two of the five networks MCMRA and in another two MLMSMRA performed best, and for one network both predictions were equal. Thus it is likely that specific network characteristics have to be considered in order to choose between MCMRA and MLMSMRA.

It should be noted that although BVSA, MCMRA and MLMSMRA performed similarly well on our transcription factor networks, in the DREAM Challenge BVSA performs best on knockout as well as on knockdown data.

Discussion

MRA represents a potentially powerful framework to recover an underlying network structure from systematic perturbation data of medium-sized networks. Several approaches have been extending MRA so that it can be efficiently applied to a wide range of experimental datasets. In this work, we review and benchmark a variant of MRA that uses a maximum-likelihood framework and a greedy hill-climbing model selection approach, called MLMSMRA. Systematic benchmarking on in silico datasets unveiled a number of features. It was reassuring for real world applications that only sensitivity but not precision of the recovered networks is affected by signal-to-noise ratios of up to 1. This means that links in the reconstructed network are not dominated by false positives, irrespective of the noise level. Our results suggest when reducing the effective noise (by e.g. using replicate measurements) especially for the unperturbed state, much higher sensitivities can likely be obtained. Furthermore, it was found that stronger perturbations result in better predictions. This is somewhat counter-intuitive, as MRA relies on the linearisation of the underlying dynamical system and is therefore assumed to be only valid for small perturbations. Variation of non-linearity in the system, represented by co-operativity, did not strongly alter performance. It should be noted, however, that other sources of non-linearity such as saturation effects for overexpression or network-inherent features that cause bistable or oscillatory behaviour were not systematically investigated and will most likely affect MRA performance. An obstacle of applying this framework to contemporary datasets is that network size negatively correlated with overall performance. This points to the necessity of developing alternative model selection strategies such as using sparsity constraints and linear programming [27].

While most formulations of MRA require that all nodes are perturbed and all are measured, MLMSMRA could deal well with data with incomplete perturbation settings, as it can incorporate prior knowledge on network structure. Such data may be important in the presence of protein–protein interactions, as protein levels themselves are difficult to perturb without changing the mRNA level.

The benchmark also shows that MLMSMRA can handle different interaction types and additional data, often better than the other benchmarked approaches (BVSA, MCMRA, plain MRA and NIR). This is interesting, as NIR has for a long time been considered the most efficient algorithm for reconstructing networks from perturbation data [8,9,2]. It further was shown that MLMSMRA performs reasonably well also on RNA data for sparse, small and non-linear networks. Unbiased benchmarking using the DREAM challenge confirms that MRA-based methods are competitive to reconstruct the networks from perturbation data, and it also verifies that larger perturbations give better results. The DREAM benchmark reproduces previous findings that BVSA performs best among MRA variants [28,22] which is surprising given that the benchmark on the in silico transcription factor networks shows equal performance for MLMSMRA, MCMRA and BVSA. The reasons might be a more beneficial network setting as elaborated for MCMRA in Supplementary Section A7 or explained by the fact that our benchmark was done on thresholded data and required sign information that was not provided by BVSA itself. In the light of these results when reverse engineering transcription factor networks from RNA measurements, BVSA can be recommended.

If proteins are measured, hidden layers are present or a more complex perturbation scheme is present, MLMSMRA could excel. MLMSMRA was initially developed to dissolve the regulatory structures that lie underneath the oncogenic transformation of ovarian cancer cells upon KRAS mutation [17]. In that study, a network, consisting of RAS proximal signalling upstream of seven transcription factors, was systematically perturbed and RNA and proteins were measured. Of the 25 interactions predicted by MLMSMRA, 15 were confirmed by independent literature findings whereas only one false positive interaction could be ascertained. Our in silico benchmarks suggests that this outstanding performance was due to the beneficial experimental setup: (i) reduced noise influence due to redundant measurements of mRNA and protein, (ii) strong perturbations with knockdown efficiencies of 80% and (iii) a rather conservative error estimate of 0.4 for the entries of R probably limiting false positive.

While the real-world example shows that the approach is flexible and fit to resolve biological networks, several challenges remain ahead. One relates to the rather low precision when the experimental setup is less beneficial than in the example of the RAS network. The other challenge will be upscaling of these approaches, which will be required to leverage the huge amount of data that is expected from single-cell CRISPR screens [20]. While this is also an algorithmic challenge, it will be essential to include prior knowledge from systematic compilation of available literature data or curated databases such as TRANSFAC [29], ORegAnno [30], HTRIdb [31] or FasTForward DNA [32]. Bayesian approaches, such as those previously used for signalling, may help to integrate such knowledge [28].

In the present study we have strived to create biologically valid test scenarios. This can further be improved by going in vitro studying synthetic networks in living cells. The current problem is that the realised networks are very small ranging from two to three nodes in human HEK293 cells [33,34] and up to five nodes in yeast [21]. However, even from those small networks interesting aspects of MRA could be deduced. For example Kang et al. [33] found that the quantification of local response coefficients was partly sensitive to perturbation strength in a three-node feed forward motif but insensitive in a derived cascade motif.

Next to more realistic settings the main challenge for a more comprehensive understanding of these networks will be to model and appreciate the dynamic nature of the gene regulatory network, where interactions are context- and time-dependent. Obviously, steady state approaches can not capture these features, but may still be informative if applied to cells in multiple states.

Summary

  • MRA is a method suitable to disentangle direct from indirect interactions.

  • Several extensions of MRA have been developed to apply it to biological data with noisy and/or incomplete measurements.

  • Here we describe one such extension (the MLMSMRA approach) in more detail on transcriptional networks.

  • We benchmark multiple MRA approaches on simulated data and highlight their respective strengths.

Funding

This work was supported partially by the BMBF (German Ministry for Education and Research) grant MapTorNet [grant number 031A426A (to N.B.)].

Competing interests

The authors declare that there are no competing interests associated with the manuscript.

Abbreviations

     
  • BVSA

    Bayesian variable selection approach

  •  
  • DREAM

    dialogue for reverse engineering assessment and method

  •  
  • mcc

    Matthew’s correlation coefficient

  •  
  • MRA

    modular response analysis

  •  
  • NIR

    network identification by multiple regression

  •  
  • ODE

    ordinary differential equation

  •  
  • pp

    protein–protein

  •  
  • pR

    protein–RNA

  •  
  • RSA

    regulatory strength analysis

  •  
  • TF

    transcription factor

References

References
1
Barenco
M.
,
Tomescu
D.
,
Brewer
D.
,
Callard
R.
,
Stark
J.
and
Hubank
M.
(
2006
)
Ranked prediction of P53 targets using hidden variable dynamic modeling
.
Genome Biol.
7
,
R25
[PubMed]
2
Stark
J.
,
Brewer
D.
,
Barenco
M.
,
Tomescu
D.
,
Callard
R.
and
Hubank
M.
(
2003
)
Reconstructing gene networks: what are the limits?
Biochem. Soc. Trans.
31
,
1519
1525
[PubMed]
3
Geier
F.
,
Timmer
J.
and
Fleck
C.
(
2007
)
Reconstructing gene-regulatory networks from time series, knock-out data, and prior knowledge
.
BMC Syst. Biol.
1
,
11
[PubMed]
4
Morris
M.K.
,
Saez-Rodriguez
J.
,
Clarke
D.C.
,
Sorger
P.K.
and
Lauffenburger
D.A.
(
2011
)
Training signaling pathway maps to biochemical data with Constrained Fuzzy Logic: quantitative analysis of liver cell responses to inflammatory stimuli
.
PLoS Comput. Biol.
7
,
e1001099
[PubMed]
5
Saez-Rodriguez
J.
,
Alexopoulos
L.G.
,
Epperlein
J.
,
Samaga
R.
,
Lauffenburger
D.A.
,
Klamt
S.
et al (
2009
)
Discrete logic modelling as a means to link protein signalling networks with functional analysis of mammalian signal transduction
.
Mol. Syst. Biol.
5
,
331
[PubMed]
6
Gardner
T.S.
,
di Bernardo
D.
,
Lorenz
D.
and
Collins
J J
(
2003
)
Inferring genetic networks and identifying compound mode of action via expression profiling
.
Science
301
,
102
105
[PubMed]
7
Santos
S.D.M.
,
Verveer
P.J.
and
Bastiaens
P.I.H.
(
2007
)
Growth factor-induced MAPK network topology shapes Erk response determining PC-12 cell fate
.
Nat. Cell Biol.
9
,
324
330
[PubMed]
8
Brazhnik
P.
(
2005
)
Inferring gene networks from steady-state response to single-gene perturbations
.
J. Theor. Biol.
237
,
427
440
[PubMed]
9
Camacho
D.
,
Licona
P.V.
,
Mendes
P.
and
Laubenbacher
R.
(
2007
)
Comparison of reverse-engineering methods using an in silico network
.
Ann. N.Y. Acad. Sci.
1115
,
73
89
[PubMed]
10
de la Fuente
A.
and
Mendes
P.
(
2002
)
Quantifying gene networks with regulatory strengths
.
Bioinformatics
29
,
73
77
[PubMed]
11
Andrec
M.
,
Kholodenko
B.N.
,
Levy
R.M.
and
Sontag
E.
(
2005
)
Inference of signaling and gene regulatory networks by steady-state perturbation experiments: structure and accuracy
.
J. Theor. Biol.
232
,
427
441
[PubMed]
12
Bruggeman
F.J.
,
Westerhoff
H.V.
,
Hoek
J.B.
and
Kholodenko
B.N.
(
2002
)
Modular response analysis of cellular regulatory networks
.
J. Theor. Biol.
218
,
507
520
[PubMed]
13
Kholodenko
B.N.
,
Kiyatkin
A.
,
Bruggeman
F.J.
,
Sontag
E.
,
Westerhoff
H.V.
and
Hoek
J.B.
(
2002
)
Untangling the wires: a strategy to trace functional interactions in signaling and gene networks
.
Proc. Natl. Acad. Sci. U.S.A.
99
,
12841
12846
14
Heinrich
R.
and
Rapoport
T.A.
(
1974
)
A linear steady-state treatment of enzymatic chains. General properties, control and effector strength
.
Eur. J. Biochem.
42
,
89
95
[PubMed]
15
Kacser
H.
and
Burns
J.A.
(
1973
)
The control of flux
.
Symp. Soc. Exp. Biol.
27
,
65
104
[PubMed]
16
Kholodenko
B.N.
,
Hoek
J.B.
,
Westerhoff
H.V.
and
Brown
G.C.
(
1997
)
Quantification of information transfer via cellular signal transduction pathways
.
FEBS Lett.
414
,
430
434
[PubMed]
17
Stelniec-Klotz
I.
,
Legewie
S.
,
Tchernitsa
O.
,
Witzel
F.
,
Klinger
B.
,
Sers
C.
et al (
2012
)
Reverse engineering a hierarchical regulatory network downstream of oncogenic Kras
.
Mol. Syst. Biol.
8
,
601
[PubMed]
18
Timmer
J.
,
Müller
T.
,
Swameye
I.
,
Sandra
O.
and
Klingmüller
U.
(
2004
)
Modeling the nonlinear dynamics of cellular signal transduction
.
Int. J. Bifurcat. Chaos
14
,
2069
2079
19
Baldi
P.
,
Brunak
S.
,
Chauvin
Y.
,
Andersen
C.A.
and
Nielsen
H.
(
2000
)
Assessing the accuracy of prediction algorithms for classification: an overview
.
Bioinformatics
16
,
412
424
[PubMed]
20
Datlinger
P.
,
Rendeiro
A.F.
,
Schmidl
C.
,
Krausgruber
T.
,
Traxler
P.
,
Klughammer
J.
et al (
2017
)
Pooled CRISPR screening with single-cell transcriptome readout
.
Nat. Methods
14
,
297
301
[PubMed]
21
Cantone
I.
,
Marucci
L.
,
Iorio
F.
,
Ricci
M.A.
,
Belcastro
V.
,
Bansal
M.
et al (
2009
)
A yeast synthetic network for in vivo assessment of reverse-engineering and modeling approaches
.
Cell
137
,
172
181
[PubMed]
22
Santra
T.
,
Kolch
W.
and
Kholodenko
B.N.
(
2013
)
Integrating Bayesian variable selection with modular response analysis to infer biochemical network topology
.
BMC Syst. Biol.
7
,
57
[PubMed]
23
Marbach
D.
,
Schaffter
T.
,
Mattiussi
C.
and
Floreano
D.
(
2009
)
Generating realistic in silico gene networks for performance assessment of reverse engineering methods
.
J. Comput. Biol.
16
,
229
239
[PubMed]
24
Marbach
D.
,
Prill
R.J.
and
Schaffter
T.
,
Mattiussi
C.
,
Floreano
D.
,
Stolovitzky
G.
(
2010
)
Revealing strengths and weaknesses of methods for gene network inference
.
Proc. Natl. Acad. Sci. U.S.A.
107
,
6286
6291
25
Prill
R.J.
,
Marbach
D.
,
Saez-Rodriguez
J.
and
Sorger
P.K.
,
Alexopoulos
L.G.
,
Xue
X.
et al (
2009
)
Towards a rigorous assessment of systems biology models: the DREAM3 challenges
.
PLoS ONE
5
,
e9202
26
Küffner
R.
,
Petri
T.
,
Windhager
L.
and
Zimmer
R.
(
2010
)
Petri Nets with Fuzzy Logic (Pnfl): reverse engineering and parametrization
.
PLoS ONE
5
,
27
Bosdriesz
E.
,
Prahallad
A.
,
Klinger
B.
,
Sieber
A.
,
Bosma
A.
,
Bernards
R.
et al (
2018
)
Comparative network reconstruction using Mixed Integer Programming
.
bioRxiv
,
28
Halasz
M.
,
Kholodenko
B.N.
,
Kolch
W.
and
Santra
T.
(
2016
)
Integrating network reconstruction with mechanistic modeling to predict cancer therapies
.
Sci. Signal.
9
,
ra114
[PubMed]
29
Matys
V.
,
Kel-Margoulis
O.V.
,
Fricke
E.
,
Liebich
I.
,
Land
S.
,
Barre-Dirrie
A.
et al (
2006
)
TRANSFAC and its module TRANSCompel: transcriptional gene regulation in eukaryotes
.
Nucleic Acids Res.
34
,
D108
D110
[PubMed]
30
Griffith
O.L.
,
Montgomery
S.B.
,
Bernier
B.
,
Chu
B.
,
Kasaian
K.
,
Aerts
S.
et al (
2008
)
ORegAnno: an open-access community-driven resource for regulatory annotation
.
Nucleic Acids Res.
36
,
D107
D113
[PubMed]
31
Bovolenta
L.A.
,
Acencio
M.L.
and
Lemke
N.
(
2012
)
HTRIdb: an open-access database for experimentally verified human transcriptional regulation interactions
.
BMC Genomics
13
,
405
[PubMed]
32
Thomas
P.
,
Durek
P.
,
Solt
I.
,
Klinger
B.
,
Witzel
F.
,
Schulthess
P.
et al (
2015
)
Computer-assisted curation of a human regulatory core network from the biological literature
.
Bioinformatics
31
,
1258
1266
[PubMed]
33
Kang
T.
,
Moore
R.
,
Li
Y.
,
Sontag
E.
and
Bleris
L.
(
2015
)
Discriminating direct and indirect connectivities in biological networks
.
Proc. Natl. Acad. Sci. U.S.A.
112
,
12893
12898
[PubMed]
34
Kang
T.
,
White
J.T.
,
Xie
Z.
,
Benenson
Y.
,
Sontag
E.
and
Bleris
L.
(
2013
)
Reverse engineering validation using a benchmark synthetic gene circuit in human cells
.
ACS Synthetic Biol.
2
,
255
262
,

Supplementary data