Impact of measurement noise, experimental design, and estimation methods on Modular Response Analysis based network reconstruction

Modular Response Analysis (MRA) is a method to reconstruct signalling networks from steady-state perturbation data which has frequently been used in different settings. Since these data are usually noisy due to multi-step measurement procedures and biological variability, it is important to investigate the effect of this noise onto network reconstruction. Here we present a systematic study to investigate propagation of noise from concentration measurements to network structures. Therefore, we design an in silico study of the MAPK and the p53 signalling pathways with realistic noise settings. We make use of statistical concepts and measures to evaluate accuracy and precision of individual inferred interactions and resulting network structures. Our results allow to derive clear recommendations to optimize the performance of MRA based network reconstruction: First, large perturbations are favorable in terms of accuracy even for models with non-linear steady-state response curves. Second, a single control measurement for different perturbation experiments seems to be sufficient for network reconstruction, and third, we recommend to execute the MRA workflow with the mean of different replicates for concentration measurements rather than using computationally more involved regression strategies.


S6 -Estimation methods to handle multiple replicates and solve the linear regression problem 22
S7 -Experiment design: number of replicates 27

S1 -ODE test-bed models
Suppl. Fig. S1 shows the two test-bed models which are used in this study. Both models are based on previously published models for the MAPK [REFKholodenko2006] and p53 [REFFey2016] systems. Parameters have been chosen in order to obtain an EGF-induced MAPK model that behaves relatively linear and a p53 DNA-damage response model that behaves strongly nonlinear. This reflects the observed behavior of these systems [REFKholodenko2010], [REFPurvis2012].
A model of signal transduction of the MAPK pathway upon EGF stimulation is illustrated in Suppl. Fig. S1a. It consists of a three-tiered cascade of phosphorylation-dephosphorylation cycles in which pRaf phosphorylates and thereby activates MEK, which then activates ERK, which negatively feeds back to Raf. Both MEK and ERK require phosphorylation at two sites to become fully activated, which is for simplicity assumed to happen in a single reaction step for both proteins. This system is described by a dynamical model in the form of Ordinary Differential Equationsẋ = f (x, θ ), x ∈ R 3 + . The state variables x(t) = [x 1 (t), x 2 (t), x 3 (t)] refer to the active states of the three proteins pRaf, ppMEK and ppERK. The input u(t) = EGF is set to a constant value u(t) = 1. The initial conditions of the state variables are set to zero, i.e. there are no active states at time t = 0. The total concentrations of the three states are Raf TOT , MEK TOT , ERK TOT = 20. The kinetic rates and other parameter values are: k 1 = 5, K + m1 = 20,V m1 = 10, K − m1 = 20, k 2 = 3, K + m2 = 20,V m2 = 10, K − m2 = 20, k 3 = 1, K + m3 = 20,V m3 = 10, K − m3 = 20, K m f = 5. Perturbation parameter p 1 , p 2 and p 3 correspond to fold changes in total protein amounts. One after the other they are set to some defined value smaller or greater than 1, the same for all perturbation experiments. To obtain the simulation results presented in the manuscript, we considered the values p j ∈ {0.2, 0.5, 0.75, 1.5}, ∀ j = 1, 2, 3, that represent 80%, 50% and 25% knockdown or 50% overexpression of the total protein concentrations.
Suppl. Fig. S1b illustrates activation of p53 by pATM. Active p53 triggers expression of MDM2, which is in turn involved in the degradation of p53, resulting in a negative feedback loop. Both p53 and MDM2 are subject to synthesis and degradation, and model variables x 1 , x 2 and x 3 correspond to pATM, p53 total amount and MDM2, respectively. As before, perturbation parameters p 1 , p 2 and p 3 describe fold changes in ATM total amount, and in p53 and MDM2 synthesis rates. The initial conditions of the state variables are set to zero, i.e. there are no active states at time t = 0. The input is set to a constant value u(t) = 1. The kinetic rates and other parameter values are: k 1 = 3, K 1 = 0.5, k 2 = 5, K2 = 0.5, k 3 = 1, n 5 = 5, K 5 = 0.1, k 4 = 1, KD = 0.01, k 6 = 1, n 6 = 5, K 6 = 0.5, k 7 = 1. The total concentration of ATM is ATM TOT = 1.
In both subfigures, graphs indicate the steady states of the system variables as functions of the perturbation parameters. While these curves can well be approximated by linear functions in the first case, they show a more pronounced non-linear behaviour in the second case. Figure S1. Two test-bed models. Shown are the reaction kinetic schemes (left), the ODE system (right top) and the dependencies of the logarithm of the state variables on the perturbation parameters (right bottom) for (a) the MAPK system; (b) the p53 system.

S2 -Linear regression model for the LRCs
For N = 3 we can rewrite equation (5) of the manuscript as the following set of 3 linear systems: Each of the three systems with two equations in two independent variables can be rewritten in the following matrix form: By stacking all unknownsr i j , i, j = 1, 2, 3, in one single vector x, we obtain the following 6 dimensional linear system y = A · x: S3 -Variation of the steady-state variables, GRCs and LRCs for changing perturbation strengths Figure S2. Nonlinear variation of the steady-state variables for changing perturbation strengths for the MAPK test-bed model. The variability of the measured steady statesz i , i = 1, 2, 3 is shown in the control experiment (p j = 1, j = 1, 2, 3) and in the 50% knockdown experiments (p j = 0.5, j = 1, 2, 3). We generated n = 10, 000 realizations via Monte Carlo simulations from the noise model (9) with parameters σ η = 0.1 and σ ε = 0.2. The continuous solid lines correspond to the noise free trend of the logarithm of the steady states for changing perturbation parameters (lnx i (p j )).

Figure S3
. Nonlinear variation of the GRCs for changing perturbation strengths for the MAPK test-bed model. The variability of the GRCs R i j , i, j = 1, 2, 3 is obtained from the sampled steady-state realizations from Fig. S2 in the case of 50% knockdown experiments (p j = 0.5, j = 1, 2, 3). We generated n = 10, 000 realizations via Monte Carlo simulations from the noise model (9) with parameters σ η = 0.1 and σ ε = 0.2. The continuous solid lines correspond to the noise free trend of the approximation of the GRCs for changing perturbation parameters (R i j (p j )). Figure S4. Nonlinear variation of the LRCs for changing perturbation strengths for the MAPK test-bed model. The variability of the LRCs r i j , i, j = 1, 2, 3 is obtained from the sampled steady-state realizations from Fig. S2 in the case of 50% knockdown experiments (p j = 0.5, j = 1, 2, 3). We generated n = 10, 000 realizations via Monte Carlo simulations from the noise model (9) with parameters σ η = 0.1 and σ ε = 0.2. The continuous solid lines correspond to the noise free trend of the approximation of the LRCs for changing perturbation parameters (r i j (p j )).   Figure S7. Nonlinear variation of the steady-state variables for changing perturbation strengths for the p53 test-bed model. The variability of the measured steady statesz i , i = 1, 2, 3 is shown in the control experiment (p j = 1, j = 1, 2, 3) and in the 50% knockdown experiments (p j = 0.5, j = 1, 2, 3). We generated n = 10, 000 realizations via Monte Carlo simulations from the noise model (9) with parameters σ η = 0.1 and σ ε = 0.02. The continuous solid lines correspond to the noise free trend of the logarithm of the steady states for changing perturbation parameters (lnx i (p j )). Figure S8. Nonlinear variation of the GRCs for changing perturbation strengths for the p53 test-bed model. The variability of the GRCs R i j , i, j = 1, 2, 3 is obtained from the sampled steady-state realizations from Fig. S7 in the case of 50% knockdown experiments (p j = 0.5, j = 1, 2, 3). We generated n = 10, 000 realizations via Monte Carlo simulations from the noise model (9) with parameters σ η = 0.1 and σ ε = 0.02. The continuous solid lines correspond to the noise free trend of the approximation of the GRCs for changing perturbation parameters (R i j (p j )). Figure S9. Nonlinear variation of the LRCs for changing perturbation strengths for the p53 test-bed model. The variability of the LRCs r i j , i, j = 1, 2, 3 is obtained from the sampled steady-state realizations from Fig. S7 in the case of 50% knockdown experiments (p j = 0.5, j = 1, 2, 3). We generated n = 10, 000 realizations via Monte Carlo simulations from the noise model (9) with parameters σ η = 0.1 and σ ε = 0.02. The continuous solid lines correspond to the noise free trend of the approximation of the LRCs for changing perturbation parameters (r i j (p j )).

11/35
S4 -Experiment design: choice of the perturbation strength Figure S10. Effects of different perturbation strengths on MRA based network reconstruction for the MAPK test-bed model. Comparison of the boxplots of the estimated coefficients r i j , i, j = 1, 2, 3, i = j, for different perturbation strengths: 80%, 50%, 25% knockdowns (KD) and 150% overexpression (OE) of the total protein concentrations. These distributions were obtained from n = 10, 000 realizations via Monte Carlo simulations for plausible input noise levels of both multiplicative and additive components: σ η = 0.1 and σ ε = 0.2.