Abstract

Dissecting the underlying relations is crucial to predict the impact of targeted perturbations. However, a major challenge in identifying cell-context specific signaling networks is the enormous number of potentially possible interactions. Here, we report a novel hybrid mathematical modeling strategy to systematically unravel hepatocyte growth factor (HGF) stimulated phosphoinosi-tide-3-kinase (PI3K) and mitogen activated protein kinase (MAPK) signaling, which critically contribute to liver regeneration. By combining time-resolved quantitative experimental data generated in primary mouse hepatocytes with interaction graph and ordinary differential equation modeling, we identify and experimentally validate a network structure that represents the experimental data best and indicates specific crosstalk mechanisms. Whereas the identified network is robust against single perturbations, combinatorial inhibition strategies are predicted that result in strong reduction of Akt and ERK activation. Thus, by capitalizing on the advantages of the two modeling approaches, we reduce the high combinatorial complexity and identify cell-context specific signaling networks.

Author Summary

The interconnections between signaling pathways contribute to the high complexity of signaling networks, therefore playing an important role in response to treatment in pathological conditions. Thus, unraveling the network structure in a cell-context specific manner is key to predict cellular responses to perturbations. Here, we present a novel hybrid mathematical modeling strategy taking advantage of qualitative and quantitative modeling approaches. We combine interaction graph and dynamic modeling with quantitative experimental data to study the hepatocyte growth factor induced signaling network in primary mouse hepatocytes. Specifically, we analyze the interconnections within and between PI3K and MAPK signaling pathways involved in hepatocytes proliferation. Based on literature knowledge, more than 100000 potential network structures are possible. By applying our approach, we reduce this combinatorial compleXity and select 16 minimal model structures. Subsequently, by performing a systematic model selection we select the model structure representing the eXperimental data best. We eXperimentally validate the resulting best model structure, and, based on model simulation, we are able to predict the outcome of combinatorial treatments. Our hybrid approach is applied to unravel cell-context specific network structures and to predict the outcome of intervention strategies.

Introduction

Traditionally, signaling cascades were interpreted as linear chains of events. However, signaling pathways involve extensive crosstalk and feedforward as well as feedback loops resulting in complex, nonlinear intracellular signaling networks, whose topologies are often context-specific and altered in diseases [1].

HGF is the key growth and survival factor for hepatocytes [2, 3] and in response to liver damage facilitates restoration of the tissue mass by promoting proliferation of hepatocytes. Upon binding to its receptor Met, HGF activates the phosphoinositide-(PI)-3-kinase (PI3K) and the mitogen activated protein kinase (MAPK) signaling pathways. Conditional knock-in mice that harbor an inactivating mutation of Met show reduced activation of PI3K signaling and complete abrogation of the activation of the MAPK pathway in response to partial hepatectomy [2]. As a consequence, damage repair is impaired in these mice suggesting an important role for the crosstalk of the signaling pathways. Therefore, to gain insights into the mechanisms controlling hepatocyte proliferation during liver regeneration, it is important to unravel mechanisms of feedback and crosstalk regulation that are relevant in hepatocytes.

Akt is activated by phosphorylation on serine 473 and threonine 308 and subsequently phosphorylates multiple substrates with important functions in key biological responses. While PI3K can be activated by direct binding to the receptor, MAPK signaling requires the activation of SOS and Ras that in turn activates Raf initiating the MAPK signaling cascade. Activated Raf leads to phosphorylation of a dual specific kinase, the mitogen-activated protein kinase kinase (MEKI and 2), that phosphorylates the extracellular-signal regulated kinase (ERK1 and 2). Dual phosphorylated ERK regulates cytoplasmic and nuclear factors and thereby modulates numerous biological responses such as proliferation, differentiation and survival. Both signaling pathways are tightly interlinked and several mechanisms have been proposed for feedback loops within and crosstalk between PI3K and MAPK signaling (S2—S3 Tables). For example, it was shown that within the MAPK signaling pathway a negative feedback loop is triggered by ERK or p9ORSK targeting SOS [4], and a positive feedback loop operates from ERK to

Furthermore, a positive feedback loop enhancing Gab1 activation via PI3,4,5-P3 generation was identified within the HGF induced PI3K signaling pathway [6]. Upon IGF induced stimulation, a negative interaction between PI3K and MAPK signaling pathways was reported, as Akt mediated phosphorylation of Raf on serine 259 leads to inactivation of Raf [7]. The majority of studies identifying these mechanisms were performed with tumor cell lines that harbor key alterations in signaling pathways controlling, for example, cell proliferation. Therefore, studies with primary cells are essential to identify physiologically relevant mechanisms. Furthermore, if we assume the 17 most likely crosstalk and feedback mechanisms between PI3K and MAPK signaling pathways (S3 Table), 131072 (217) possible network structures are conceivable. Thus, due to the high combinatorial complexity, a systematic method is required to facilitate unbiased identification of cell-context specific structure of signaling networks.

Several compounds targeting individual components in PI3K and MAPK signaling have been developed, and multiple clinical trials were initiated [8—10]. Inhibitors targeting the PI3K signaling pathway include the reversible PI3K inhibitor LY29004, the more potent irreversible inhibitor Wortmannin and the allosteric Akt inhibitor VIII. Derivatives of Wortmannin, such as PX-866, other PI3K inhibitors such as XL-147 [11] and CAL 101 [12] as well as the allosteric pan-Akt inhibitor MK-2206 [13, 14] are currently used in clinical trials. To analyze the MAPK signaling pathway, the compounds PD 98059 and U0126 inhibiting MEK have been widely applied. Several MEK inhibitors entered the clinical trials including, for example, CI-1040 and its analogue PD 0325901 [15—18]. Despite the specificity of these inhibitors, their efficacy was reported to be very limited. Combinatorial treatment with the MEK inhibitor ARRAY-438162 in combination with the PI3K inhibitors BYL719 and BKM120 or the Raf inhibitors LGX818 and RAF265 are ongoing [19, 20]. Furthermore, BX-912 is the PDK1 inhibitor that is primarily used for the analysis of this signaling pathway, whereas the PDK1 inhibitor AR 12 was tested in a clinical trial and showed limited absorption as well as pharmacokinetic variability [21]. Although multiple compounds are available, strategies to improve their application and to select most promising combinations remain to be developed.

Mathematical models of the MAPK signaling pathway have been developed that only consider negative feedback [22] , negative and positive feedback loops [5] or that analyze the signal-to-response relation [23]. Mathematical models describing both PI3K and MAPK signaling pathways upon single or combinatorial stimuli reveal the presence of crosstalk mechanisms between MAPK and PI3K pathways [24—26] or differences in the stimulus specific network topology [27, 28]. As indicated, most of the studies considered only single feedbacks or a limited number of crosstalk mechanisms. Therefore, to unravel a more complex network structure, a systematic unbiased approach is required.

Many of them are based on qualitative modeling formalisms [29—31] that can deal with large networks, but are limited in their capacity to describe dynamic properties such as signal duration and amplitude. Interaction graphs as one example for qualitative models have been shown to be valuable tools for analyzing the structure of signaling pathways as they are implicitly contained as the underlying network structure in more complex modeling formalisms [32]. They can be used to make predictions on the possible qualitative behavior of a signaling network, and these predictions can be compared with experimental data. Resulting inconsistencies between data and network structure provide then a basis to identify missing or inactive interactions in the network structure. One possibility to derive such predictions is to use the concept of the dependency matrix [32, 33]. In contrast to related methods, which rely on the concept of sign consistency and require a steady state assumption [30, 34], exploitation of the dependency matrix is well-suited for the analysis of transient effects. The general idea is that some characteristics of the possible qualitative system response are determined by the paths and feedback loops that represent the interaction graph. Interaction graphs are closely linked to logical models. They can be derived from an interaction graph by adding rules that define how the discrete state of a node is governed by the states of other nodes [33]. While logical models are well-suit-ed for studying the input-output behavior of large signaling pathways, an appropriate description of crosstalk mechanisms within the logical formalism is often difficult. This is due to the fact that, in contrast to the main activation routes of a signaling pathway, crosstalk mechanisms typically enhance or downgrade certain effects, rather than being necessary for or completely blocking them. Therefore, interaction graphs that utilize continuous states are preferable to describe crosstalk mechanisms.

However, consideration and systematic analysis of a large number of potential mechanisms results in a high combinatorial complexity with many degrees of freedom and is therefore often not feasible with ODEs. Therefore, it is desirable to exploit the advantages of both qualitative and quantitative modeling and to develop strategies to combine both approaches.

We started with an interaction graph master model containing previously reported crosstalk, feedback and feedforward mechanisms and selected then minimal model structures of the interaction graph master model that can explain the observed qualitative characteristics of the experimental data. In this way, the search space of model candidates was vastly reduced. With the subsequent analysis using ODE models, we identified the model structure representing the experimental system best. We demonstrate that the inferred HGF model shows robust behavior against single perturbations, but enables predictions of combinatorial inhibition leading to strongly reduced Akt and ERK activation. This application demonstrates the potential of our network inference approach to guide the development of context-specific therapeutic intervention strategies.

Results

HGF induced signaling pathways

We distinguish between interactions that form the main activation routes (“core model”; black edges in Fig 2) and interactions describing feedback and crosstalk mechanisms (“candidate mechanisms”; turquoise edges) that have been reported for special cell types or under certain experimental conditions (SI—$3 Tables). The full graph including the core model and all candidate mechanisms is considered as a non-cell-type specific “HGF master model”.

Additionally, we employed siRNA targeting Akt and ERK1/2. By time-resolved quantitative immunoblotting and by protein array we analyzed the phosphorylation of Akt, MEK1/2, ERK1/2 and p9ORSK (SI and S2 Fig). p9ORSK

Quantitative Interaction experimental Graph data

Master model Discretized data Selected minimal model structures Translation into ODE models

ODE MOde'ing Iterative subtraction of building blocks from

complete model Combinations of minimal models

Best model

Experimental Prediction of

Iterative addition of candidate mechanisms to core model validation inhibitor combinations

Workflow of model selection strategy. Quantitative time resolved data is discretized for the selection of submodels of the interaction graph master model. The interaction graph master model is based on literature knowledge and consists of the “core model” and reported interactions between the signaling pathways of interest, the “candidate mechanisms”. Ordinary differential equation (ODE) modeling utilizes the entire information of the time resolved data. Model selection based on parameter estimation permits the selection of the best model structure.

We also measured the dynamics of 8081 activation upon MEK inhibitor treatment (SlA Fig). The fold change of protein phosphorylation for each treatment condition was calculated in comparison to the respective control and between treatment conditions (Fig 3A). For 8081 activation, the band shift was quantified as an indicator for its activation status (81A and 83A Figs). As eXpected, upon Met inhibitor treatment a strong decrease in phosphorylation of all measured proteins was detected. Upon inhibition of MEK, MEK phosphorylation was initially decreased, followed by an increased and sustained signal, whereas ERK, RSK_s and RSK_d phosphorylation was strongly reduced. MEK inhibitor also influenced the dynamics of 8081 activation resulting in a more sustained activation (81A and 83A Figs). During the entire observation period the effect of MEK inhibitor treatment on Akt phosphorylation was highly variable between the eXperiments. The siRNA targeting ERK1/2 showed a comparable effect on ERK and Akt as the MEK inhibitor. Treatment with PDKl inhibitor reduced Akt phosphoryla-tion on serine 473 and threonine 308 (Fig 3A and 83A Fig). Interestingly, a decreased phosphorylation of MEK, ERK, RSK_s and RSK_d upon PDKl inhibitor treatment was observed at early time points and a subsequent increase at later time points. When applying MEK and PDKl inhibitors simultaneously, we observed decreased phosphorylation of all measured proteins, but MEK phosphorylation showed an increased signal at later time points. A comparable phosphorylation response was elicited by combinatorial treatment with MEK and PDKl inhibitors or with PDKl inhibitor alone except for Akt phosphorylation, which was additionally decreased upon the combinatorial inhibitor treatment. The comparison between the combinatorial treatment with MEK and PDKl inhibitor and MEK inhibitor treatment alone showed decreased phosphorylation on Akt, MEK and RSK_d. Treatment with PI3K inhibitor and with siRNA targeting Akt were not considered in our study due to the high variability observed in our results (83B Fig).

Selection of minimal model structures

To relate the data to the interaction graph model (Fig 2), we discretized the measured responses to “increase”, “decrease” and “not measured or not conclusive” (Fig 3B): If a certain effect was consistently observed for the replicates within the same time points, it was included as an observed effect in the scheme of the discretized data. Otherwise, the effect was considered “not conclusive” and not taken into account. The fourth possible value "no change" was not observed in our data set. Additionally, we included the timing of the response in the discretized data, classifying an observed effect as an early or late response. response"Early" refers to the initial qualitative re-sponse of a node within 30 minutes after HGF stimulation. A response was termed “late” if the qualitative behavior at any successive time point is different from the initial (early) response (Materials and Methods). Thus, the late response characterizes the first effect that is opposite to the early response, or indicates that the qualitative response is similar for all time points.

If the model predictions are in accordance with the discretized data, the given structure is able to reflect the experimental results. Fig 3B shows the possible qualitative responses predicted from the interaction graph core model. Comparing model predictions with the discretized data, the majority of observed behaviors were not represented by the core model. Hence, we conclude that some of the candidate mechanisms must be active in primary mouse hepatocytes. To identify minimal substructures of the HGF master model consistent with the data, we added single candidate mechanisms and combinations thereof to the core model and tested the resulting model structures for their ability to represent the observed effects. In order to select all minimal submodels of the interaction graph master model that can explain the observed effects from the experimental data (Fig 3B, left panel) and that contain the core model, we proceeded as follows: starting from the core model, we added one candidate mechanism at a time (that is, all edges making up the mechanism, S3 Table) and derived the model predictions for the resulting interaction graph structure as explained above. If all experimental observations were in accordance with these predictions, the structure was added to the list of selected minimal model structures. If not, we combined the respective structure with other candidate mechanisms whose sole addition to the core model was not able to explain the data. Again, we checked the ability of each resulting model structure to explain the data and either considered the structure as selected minimal model or added further candidate mechanisms. We ensured that only minimal model structures were selected, that is, no submodel of a selected model can explain the data. Furthermore, the procedure assures that all minimal model structures that are able to explain the data were found.

Thus, among the possible candidate edges that may be relevant for the other cell types and conditions, we selected candidate mechanisms that are specific for the HGF induced MAPK and PI3K signaling pathways in primary mouse hepatocyte. The candidate edges that have not been selected by our approach might play a role in other cell types and for other conditions. In total, we identified 16 minimal model structures that can equally well explain our experimental data (Fig 3B, right panel and S4 Fig). Notably, in some cases, the qualitative response is not restricted by the model structure: If a path between the inhibited and measured species includes a negative feedback loop, or if paths of both signs exist between the two nodes, the actual response depends on the strength of the different mechanism and, thus, cannot be predicted by a purely qualitative model.

Each candidate mechanism is represented by a single edge. Hence, each selected minimal model structure is composed of the compressed core model and a combination of candidate edges, each corresponding to one candidate mechanism. We define as “building block” the characteristic combination of candidate edges of a selected minimal model structure. Fig 4 shows the compressed selected 16 minimal model structures, the compressed core model and a model containing all building blocks (“complete model”). In contrast to the master model, the complete model contains only those candidate mechanisms that are included in at least one selected minimal model structure. By comparing the 16 minimal model structures, we observed that all candidate models include (i) the edge from ERK to PI3K, (ii) a negative feedback from ERK to $081, either directly or indirectly via RSK_d, (iii) a positive route from ERK to MEK through various mechanisms and (iv) a positive route from PDKI to MEK. The latter route is included as a direct edge from PDKI to MEK in all models but model 8, Whereas model 8 contains a longer path via RSK_d, Ras and Akt.

Ordinary differential equation model selection

To test Which of the identified structures can quantitatively represent the transient and sustained effects observed in the experimental data, we translated each of the 16 compressed selected minimal model structures as well as the core and the complete model into an ODE model. While interaction graph modeling is based on discretized data, we utilized the full quantitative information of the data sets for ODE modeling (S1—S4 Datasets). In addition to the datasets used for the selection of the 16 model structures, we used datasets of time resolved measurements of HGF induced activation of the above mentioned proteins without inhibitor treatment. Furthermore, Met phosphorylation and active Ras measured by quantitative immu-noblotting and protein array, respectively, were considered. In addition, the degree of phos-phorylation of ERK and Akt after HGF stimulation was measured by mass spectrometry and included in the model. Due to its merely qualitative nature, the dataset on 8081 activation was not used for ODE modeling. Overall, the ODE models were calibrated on 2200 data points and 25 experimental conditions including four targeted perturbations. A complete list of kinetic reactions considered for our modeling approach is shown in S4 Table. For each model structure, parameter estimation was performed to determine the model performance in relation to the experimental data. We applied an adaptation of a likelihood ratio test (LRT) with a threshold of 95%, which takes into account the different degrees of freedom for each model structure (Materials and Methods). This “forward selection” facilitated an LRT-based ranking of model structures (Fig 5A). As expected, the complete model structure performed best, whereas the core model structure performed significantly worse than all selected minimal model structures. Interestingly, the model ranking showed several selected minimal models with similar likelihood (16, 4 and 10). Therefore, a clear distinction of a best performing model structure was not possible. Furthermore, a significant gap was observed in the likelihood values between the complete model and the best performing selected minimal model structure, model 16. Thus, none of the minimal models is a valid simplification of the complete model. This result suggested that models containing combinations of building blocks might perform better than minimal models. To reduce the complexity of combining all model structures, we applied a “backward selection” based on the removal of single building blocks from the complete model structure followed by parameter estimation to obtain a new ranking (Fig 5B). Strikingly, the forward and backward selections revealed very different results. As an example, model 16 is the best performing selected minimal model structure in the forward selection, whereas a reduction of model 16 in the backward selection does not lead to a significant loss in likelihood. The removal of the selected minimal model structures 4, 6, 8 and 12 showed a significant loss in performance according to the likelihood, suggesting the importance of their building blocks. Based on the backward selection result, we generated combinations of model structures focusing on these four minimal models and performed parameter estimation for all eleven possible combinations (Fig 5C). As expected, the new ranking showed that these model structures filled the gap between the complete model and the best performing minimal model 16. Interestingly, combinatorial model structures 4_8_12 and 4_6_8_12 displayed a similar performance as the complete model structure, indicating that they are valid simplifications of the complete model. Between these models, model 4_8_12 performed best as it contained fewer parameters than model 4_6_8_12. The list of estimated parameter values of model 4_8_12 is given in S5 Table. A comparable ranking of all models was obtained utilizing the Akaike Information Criterion (AIC) (S6 Fig). Based on the ranking, we proposed a model structure consisting only of a subset of feedback and crosstalk mechanisms of the complete model, namely eight mechanisms out of thirteen. As shown in Fig 6B—6F, model 4_8_12 can reproduce the dynamic behavior of the measured species under diverse experimental conditions. A comparison of the performance of the core model, the complete model and the model 4_8_12 is given in S8 Fig. We hypothesized that the advantage of a reduced model resides in an improved predictive power. To compare the predictive power of the complete model and the model 4_8_12, we analyzed with the model the dynamic behaviour of a protein that has not been measured experimentally. We selected “active PI3K” that has not been used in the parameter estimation approach. We calculated the prediction profiles of “active PI3K” as described in [39] for the entire observation time frame. The prediction profiles were translated into confidence intervals along the trajectory, giving a measure of the accuracy of the prediction of the PI3K dynamic for both models. As shown in Fig 5D, model 4_8_12 outperforms the accuracy of the complete model by having an approximately 10-fold smaller confidence interval. In addition, the model trajectories of active PI3K differ significantly between the complete model and model 4_8_12. This indicates that prediction of the complete model is not only uncertain, but also incorrect as the confidence interval of the complete model does not contain the well-defined trajectory of model 4_8_12. This discrepancy is due to the parameter non-identifiabilites in the over-parameterized complete model. These results show that the selected reduced model 4_8_12 has a better predictive power than the complete model.

For the random model structures, we performed parameter estimation and derived a ranking as in the forward selection step (Fig 5E). This showed that the majority of random models performed worse than our selected models. Random models that performed well had a similar structure as the best performing selected models (S6 Table and Fig 5E). As an example, random model structure 5, which possesses three edges that are contained in model 12 or model 16, performed similarly as the selected minimal model structure 16. The two additional edges are present in model structure 4 and 8, respectively. Therefore, random model structure 5 is similar to the best performing model structure 4_8_12.

Within the MAPK pathway, a positive and a negative feedback emerge from the candidate mechanisms ERK to Raf1 and RSK_d to 8081, respectively. Within the PI3K pathway, the edge PI3K to Gabl closes a positive feedback loop. Furthermore, model 4_8_12 is characterized by the presence of several crosstalk mechanisms between PI3K and MAPK. Notably, these mechanisms give rise to two positive and two negative feedback loops, each containing species from both the PI3K and MAPK pathways. As shown in Fig 6B—6F, model 4_8_12 can reproduce the dynamic behavior of the measured species under diverse experimental conditions. A comparison of the performance of the core model, the complete model and the model 4_8_12 is given in 88 Fig.

Negative crosstalk: experimental validation

The model predictions indicate that 3, 6 and 100 fold Akt inhibition results in a 77%, 83% and 99% reduction of the inhibitory effect of Akt on Raf1 (Fig 7A and 7B). The model includes active Rafl, which cannot be directly compared to Rafl phosphorylated at a specific phosphorylation site as phosphorylation on serine 338 contributes to Rafl activation, whereas phosphorylation on serine 259, the target of Akt, represents an inactivation signal [40]. Therefore, we used the effect of Akt inhibition on serine 259 phosphorylation of Rafl as a proxy to experimentally address the inhibitory impact of Akt on Rafl. Experimentally, we stimulated primary mouse hepatocytes with HGF in the absence and presence of the specific Akt inhibitor VIII and monitored the impact on Rafl phosphorylation on serine 259. As shown in Fig 7C and 89—811 Figs, we achieved an inhibition of Akt phosphorylation between 90% and 100% and concomitantly observed a decrease of Rafl phosphorylation on serine 259 (Fig 7D and 89—811 Figs). Additionally, a moderate increase of MEK phosphorylation is also observed upon Akt inhibitor treatment (Fig 7D). The results confirmed the presence of this interaction in our cellular model system upon HGF stimulation. This interaction is also present in models 7 and 8, but these models have not been selected as best performing models. We emphasize that, since model 4_8_12 is derived from the combination of the single models 4, 8 and 12, it does not contain unique edges. However, the validation of model 4_8_12 as the optimal network structure in primary mouse hepatocytes is given by our model selection process together with the experimental validation of the interaction between Akt and Raf1.

These characteristics may limit the efficacy of targeted therapies or create undesired effects. Therefore, to further validate our model structure, we analyzed how the system responds to different perturbation conditions.

Inhibitor combination: model predictions and experimental validation

Additionally, we performed simulations for every possible combination of double inhibition. As readout of PI3K and MAPK pathway activation, we calculated the area under the curve of Akt and ERK phosphorylation and the sum thereof to evaluate the effect of the targeted inhibition (Fig 8A). Single inhibition of PI3K or PDK1 drastically reduced pAkt and, surprisingly, exerted the opposite effect on pERK. This result indicated that the network structure was robust against single inhibitor treatment suggesting that the two signaling pathways compensated each other. We calculated the sum of the integrals of Akt and ERK phosphorylation to predict the double inhibition resulting in an effective reduction of both readouts (Fig 8A). Interestingly, while the sum of the integral of pAkt and pERK is sensitive to the combined inhibition of PI3K and MEK resulting in a 71% reduction of the sum, the combined inhibition of Met and RSK results in an increase of the sum to 118%. In some cases, we observed synergistic effects as the combination of inhibitors had a stronger effect on the readout than the sum of impact of the respective individual inhibitions. For example, single inhibition of Met or SOS1 led to the reduction of the sum of the integral of pAkt and pERK by 18%, while their combination led to 50% reduction. To address the synergistic effect of the inhibitor combination treatment, we calculated the effect of the inhibitor combination compared to the single inhibitor treatment (Fig 8B). Interestingly, we observed that the Ras inhibitor has a synergistic effect on pAkt and pERK with several other inhibitors such as Raf1, ERK and Rac inhibitors. We furthermore detected synergy between the PDK1 and Akt inhibitor on pAkt and between the Met inhibitor and the SOS1 inhibitor on pAkt and pERK. Surprisingly, the combination of PI3K and PDK1 inhibitors shows the lowest synergistic effect on pAkt and the highest on pERK, resulting in a mild synergy in the sum of pAkt and pERK.

We estimated 10% strength for the PI3K inhibitor, 83% for the MEK inhibitor, 90% strength for the Met inhibitor, and 53% for the PDK inhibitor. Based on the estimated inhibitor parameters, we performed predictions of the dynamic behavior of pAkt and pERK upon combinatorial inhibitor treatments. In detail, we simulated the activation kinetics of pAkt and pERK upon combining the PI3K and the MEK inhibitor and upon combining the Met and the PDK1 inhibitor. We experimentally validated the model predictions by treating primary mouse hepatocytes with the indicated inhibitor doses and combinations and analyzed the activation kinetic of pAkt and pERK (Fig 8D and 812 Fig). The experimental results indicate that the combination of low dose of the PI3K and the MEK inhibitor slightly increases pAkt and reduces pERK and, therefore, are in good agreement with the model predictions. Interestingly, the combined application of the Met and the PDK1 inhibitor resulted in a reduction of pAkt, While it had only a moderate effect on pERK. Additionally, the calculated area under the curve of pAkt, pERK and their sum for the model trajectories and the experimental data are in agreement. In conclusion, model simulations and experimental verifications suggested that the considered signaling network is less sensitive to single interventions, but can be efficiently targeted by combinatorial treatments.

Discussion

The identification and quantitative description of relevant feedback, feedforward and crosstalk regulation of signaling pathways is an important step towards understanding cellular signaling networks and a key prerequisite for the development of successful drug targeting strategies [41—43]

The key advantage of our introduced hybrid approach is that interaction graph modeling enables pre-selection of minimal model structures from a vast search space of potential model candidates, which are then translated into ODE models for integration of quantitative details. Several other mathematical modeling approaches also deal with a family of candidate models aiming at the identification of the correct wiring. These approaches employ, for example, ensemble modeling [37] or Bayesian inference [44] and can usually only deal with a limited number (tens to hundreds) of competing ODE models. However, if one considers all possible combinations of reported crosstalk, feedback and feedforward mechanisms, these possibilities result in an enormous number of candidate models. As shown herein, by preselecting minimal model structures using interaction graph modeling, one can massively reduce the search space (here from 105 to 16) before continuing with a smaller library of ODE models. Other modeling approaches using perturbation data to unravel the network structure rely on modular response analysis (MRA), which requires steady state assumptions and linear equation based modeling. Therefore, these approaches are more suited to identify feedback and crosstalk mechanisms that determine the medium to long term behavior [24, 35]. Furthermore, a complete set of perturbation eXperiments that alter the state of each individual module is required. In contrast, our proposed modeling strategy also enables the elucidation of mechanisms influencing the immediate system response as well as taking into account nonlinear effects such as saturation.

The presence of a negative feedback within the MAPK signaling is well known, and ERK is usually considered as central player of this feedback [45]. However, in line with our findings, it was shown that the negative feedback to 8081 could also be mediated by p9ORSK [4]. Although forward selection did not allow us to identify a single selected minimal model structure as the one explaining the eXperimental data best, these structures served as important basis for further analysis. By backward selection, we identified minimal model structures containing feedback and crosstalk mechanisms relevant for HGF mediated responses in primary mouse hepatocytes resulting in a reasonable number of combinations of minimal model structures that had to be tested. The best performing combinatorial model 4_8_12 showed a positive feedback edge of ERK to Raf1. This feedback is known to involve RKIP [5], and it is the result of two negative effects: RKIP negatively regulates MEK activation by binding to phosphorylated Raf1 and, therefore, inhibiting the signaling cascade activation. A negative feedback loop on RKIP is triggered by ERK allowing the release of Rafl from RKIP-Rafl complex resulting in MEK activation. This represents an additional mechanism to fine-tune signaling responses within the negative feedback from downstream proteins to 8081. Furthermore, our approach unraveled a negative interaction from Akt to Raf1, which we confirmed experimentally in our cellular model system. Rafl phosphorylation on serine 259 by Akt has also been shown in the HGF stimulated human hepatoma cell line Hep3B [46]. Raf1 phosphorylation on serine 259 contributes to control the mitogenic response induced by MAPK signaling [47], and it has been observed in the context of genetic disorders [48]. Finally, a crosstalk mechanism from ERK to PI3K is indicated in the model 4_8_12, which gives rise to two additional feedback loops controlling both the MAPK and the PI3K branch. The negative feedback loop involves the negative effect of Akt on Rafl; the positive feedback loop is mediated by PI3K activated Gabl that strengths the MAPK activation via Rac and PAK.

Feedback loops are the base of nontrivial dynamic behavior. Positive feedback loops have been shown to cause bistable behavior [49] and must be contained in the system's structure to enable more than one steady state [50]. Negative feedback loops are known to stabilize the system's response [51] and are a structural prerequisite for oscillations [52]. Within the last decades, it has become more and more evident that complex and robust dynamics of cellular signaling networks arise from the interplay of positive and negative feedback loops [53]. Examples include that of periodic calcium spikes after growth factor or hormone stimulation [54] and the MAPK system, where positive and negative feedback loops allow the system to switch between monostable and bistable regimes and, therefore, to flexibly adapt the response to different stimuli [55].

A number of potential interactions between these two pathways have been studied in different cell types (S3 Table). By employing mathematical modeling, a study of the MAPK and PI3K pathways crosstalk showed that both compensate for each other [28]. Although these two pathways are well studied, a systematic strategy to unravel their crosstalk has been missing.

In our case study, various interconnections between the PI3K and MAPK pathway in the best performing model structure indicate redundancy within the HGF stimulated signaling network. Our model predictions suggest that the effects of intervention in one signaling pathway can be compensated by the impact of the other pathway, as in case of treatment with single PI3K or PDKl inhibitor. To maximize the response of inhibition, combinatorial treatments are required as suggested previously [24, 36]. However, a positive synergistic effect is predicted only for some combinatorial treatments suggesting robustness. Therefore, analyses focusing on the identification of the most promising synergetic effects are essential to develop novel combinatorial treatments. Our proposed strategy to disentangle complex signaling networks provides a generic and systematic approach to elucidate cell and context-specific feedback and crosstalk mechanisms and promises to significantly improve the development of effective combinatorial intervention strategies.

Materials and Methods

Ethics statement All animal experiments were approved by the governmental reView committee on animal care of the state Baden Wiirttemberg, Germany (reference number A24/ 10). Anesthesia was carried out by intraperitoneal injection of 5 mg ketamine hydrochloride 10% (w/V) (Bayer Health

Isolation of primary mouse hepatocytes

Mice, which were housed at the DKFZ animal facility under a constant light/ dark cycle, maintained on a standard mouse diet and allowed ad libitum access to food and water, were used for primary hepatocyte isolation. Primary mouse hepatocytes were isolated as described [58]. After isolation, hepatocytes were seeded in full medium (phenol red-free Williams E medium (Biochrom) supplemented with 10% (v/v) fetal bovine serum (Life Technologies), 0.1 uM dexamethasone, 10 ug/ml insulin, 2 mM L-glutamine (Life Technologies) and 1% (v/v) penicillin/streptomycin 100x (Life Technologies)) using collagen I-coated cell dishes (BD Biosciences). Hepatocytes were cultured at 37°C, 5% C02 and 95% rH. Following cell adhesion, hepatocytes were washed with PBS (PAN Biotech) and subsequently cultivated in serum-free cultivation medium (phenol red-free Williams E medium supplemented with 0.1 microM dexamethasone, 2 mM L-glutamine and 1% (v/v) penicillin/streptomycin 100x) for 14 hours. Subsequently, hepatocytes were washed with PBS (PAN Biotech) and cultivated in serum-free cultivation medium depleted of dexamethasone for 6 hours prior treatment.

Time course and dose response experiments

Single inhibitor treatments were also performed using the inhibitor concentrations listed before. All inhibitors were applied 30 minutes prior to stimulation with 40 ng/ml of HGF for the designated time points.

After adhesion, cells were washed as described above and serum-free cultivation medium was added. Subsequently, hepatocytes were transfected with siRNAs by applying RNAiMAX (Invitrogen) and incubated for 16 hours. ON-TARGETplus SMART Pool siRNA targeting mouse Akt, ERK1 and ERK2 (Dharmacon, Thermo Scientific) were applied to a final concentration of 10nM, (ERK1 and ERK2 ON-TARGETplus SMART Pool siR-NAs were pooled). The control condition was carried out by using ON-TARGETplus Non-targeting Pool. After 16 hours, cells were washed, and serum-free cultivation medium depleted of dexamethasone was added 6 hours prior stimulation.

Quantitative immunoblotting

Mass spectrometry

For Akt and ERK phosphorylation analysis, after 6 hours of serum-free cultivation medium depleted of dexametha-sone, the cells were stimulated With 40 ng/ml or 100 ng/ml of HGF for 10 minutes or left untreated. Total cellular lysate was denaturated in 10% SDS and subsequently subjected to IP of ERKI

For IP, anti-ERKl (Santa Cruz Biotechnology) and anti-Akt (Cell Signaling Technology) antibodies were used. Immunoprecipitated proteins were resolved on SDS-PAGE, and gels were stained with SimplyBlue SafeStain Life Technologies according to man-ufacturer’s instructions. Quantitative determination of the degree of phosphorylation of Akt and ERK was performed as previously described [59, 61] and was used for ODE model calibration.

Protein array: Akt and ERK quantification

The spotting was performed with a sciFLEX-Arrayer SP5 (Scienion, Berlin) piezoelectric non-contact spotter on 16-pad nitrocellulose slides (Oncyte, Grace). For Akt detection, the following non-rabbit-derived capture antibodies were used: sc-55523 (Santa Cruz Biotechnology) and CS2967 (Cell Signaling Technology) for total Akt, Up05669 (Upstate) for phos-phorylated Akt. For ERK detection, the following spotting antibodies were used: BD610030 (BD Bioscience), Up05157 (Upstate), MAB1576, CS9107 (Cell Signaling Technology) for total ERK and CS9106, CS9109 (Cell Signaling Technology), M9692, for phosphorylated ERK. Bovine serum albumin from Sigma (A9418) was labelled with the DyLight800 fluorophore using DyLight 800-NHS-Ester (46422) and Fluorescent Dye Removal columns (22858) from Thermo Scientific were used for the dilution of the capture antibodies for a spot to spot normalization. Capture antibodies (diluted the following: M9692 (1:8 dilution), BD, CST, Sigma, Upstate (1:3,6 dilution), Santa Cruz Biotechnology undiluted) were mixed with 2x Whatman arraying buffer (S00537), containing a 1:100 dilution of DyLight800-labelled BSA to account for spotting differences. Spotted slides were blocked with Odyssey blocking buffer (Licor 927—40000) for 6 hours shaking at 4°C.

Commercially available phosphorylated ERK (PV3313, Life Technologies); unphosphorylated ERK (PV3312, Life Technologies) and phosphorylated Akt (14—276, Millipore) and unphosphorylated Akt (14—279, Millipore) were diluted for different dilution series, supplemented with 1% SDS and denatured for 5 minutes at 95°C. Cell lysates generated as for quantitative immunoblotting were diluted 1:16 in Array Buffer Plus (Array Buffer composed by 1% BSA, 0,5% NP40, 0,02% SDS, 50mM Tris pH7,4, 150mM NaCl, 1mM EDTA, 5mM NaF, 1mM Na3VO4) containing complete protease-inhibi-tor cocktail (Roche), supplemented with SDS and denatured. Slides were incubated with calibrators and cell lysates shaking at 4°C overnight. After washes with the Assay Buffer at room temperature the slides were incubated with the detection antibodies for 2 hours shaking at 4°C. Detection antibodies were used for total ERK Up06-182 (Upstate), sc-94 (Santa Cruz Biotechnology) diluted 1:800 and total AKT sc-1619 (Santa Cruz Biotechnology) diluted 1:200, CS9272 (Cell Signaling Technology) diluted 1:400. After washes with Array Buffer and Wash Buffer (PBS (Pan, ready-made), NP-40 (0.2%), NaF (5 mM), Na-Vanadate (1 mM)) the slides were incubated with the secondary detection antibody (goat anti-rabbit coupled to Alexa Fluor 680) for 30 minutes shaking at 4°C. Slides were light-protected and washed with Wash Buffer and ddeO and subsequently dried. Slides were scanned with LiCor Odyssey scanner at 700 and 800 nm (Resolution 21 um, High Quality, Intensity: 5/4; Offset: 0, size 7x5 Scanning time approx. 35min/2slides). The quantification of the signal was performed as described [61].

Protein array: Ras quantification

cloned the Raf-RBD into the pGeX vector Via BamHI and EcoRI sites. The pGeX-GST-Raf-RBD

To spot the eluted fraction, it was diluted with PBS and glycerol. The spotting was performed as described in the previous paragraph. Cell lysates were diluted 1:10 with Array Buffer PLUS. For calibration positive and negative control samples with gam-maS-GTP and GDP loading were used. Prior to incubation, the slides were blocked with LiCor Blocking Buffer for 2—6 h. Samples and calibrator-solutions were incubated on the slides shaking overnight. All incubations were performed at 4°C. The slides were then washed with array buffer and incubated with specific rabbit-derived detection Ras antibody sc-14022 and sc-520 (Santa Cruz Biotechnology). Slides were subsequently treated as described in the previous paragraph.

Interaction graphs

The nodes in these graphs represent the signaling molecules, in our application usually the activated forms of the proteins (SI Table). Directed positive and negative edges represent activating or inhibiting influences between the proteins.

We identified all minimal submodels from the interaction graph master model that can explain all qualitative properties of our experimental data based on the concept of the dependency matrix [33]. In the dependency matrix, the effect of one species on another is classified as follows: If all paths from species A to species B are positive, A is an activator of B. If all paths from A to B are negative, A is an inhibitor of B. If positive and negative paths lead from A to B, A is said to be an ambivalent factor for B, and A is a neutral factor for B if there are no paths from A to B. Activators and inhibitors are further classified as weak or strong depending on whether any of the paths runs through a node that is involved in a negative feedback or not [33]. Weak and strong acti-vators/ inhibitors differ in how they restrict the possible qualitative behavior. Decreasing the activity of a strong activator of species B leads to a decrease of the activation of B at every time point. In contrast, decreasing the activity of a weak activator of species B also leads to an initial decrease of the activity of B, but the negative feedback may cause a subsequent increase. This concept allows deriving predictions of the possible qualitative behavior in response to inhibitors given a certain model structure. To select the minimal model structures, we used the discretized fold change (‘increase’, ‘decrease’, ‘no change’, ‘not conclusive’) of the phosphorylation state of two experimental conditions, C1 and C2. C1 represents the “control condition”, whereas in condition C2, one or two additional inhibitors compared to C1 are applied (Fig 3B, left panel). For each comparison, we checked in the dependency matrix of the model how the additional inhibitors in C2 act on the measured species leading to qualitative predictions for the early and late differential behavior of this species in condition C2 versus C1. If all additionally inhibited species in C2 are activators of the measured protein, no matter whether they are weak or strong, the model predicts "decrease" of the measured species (compared to its level in C1) as early response. Analogously, if all inhibited species act as inhibitors for a measured species, the model prediction for the early differential response of this node is "increase". For the late response, the model predicts "decrease" as the only possible behavior if all inhibited species are strong activators, and "increase" if all inhibited species are strong inhibitors. If all inhibited species are neutral factors, the model prediction is "no change", both for the early and late re-sponseresponse. If both positive and negative paths from the inhibited to the measured species exist (either because one of the inhibited species is an ambivalent factor or because one inhibited species is an activator, another one an inhibitor for the measured species) any qualitative early and late response is allowed. Finally, if there is a weak activator or weak inhibitor among the inhibited species, the late response is not restricted by the model.

The applied MEK inhibitor, for example, blocks MEK kinase activity, thus inhibiting the outgoing edges from MEK. Therefore, we introduce a “dummy” node that is activated by MEK and itself activates all downstream nodes of MEK (in our particular case only ERK). The effect of the MEK inhibitor is then reflected by this dummy node.

The pseudocodes of the algorithm are given in 81 Code. The HGF interaction graph model is provided in CellNetAnaly-zer format at http://www2.mpi-magdeburg.mpg.de/projects/cna/repositoryhtml.

Ordinary differential equation modeling

We can associate an interaction graph with a given system of ODEs: the interaction graph is defined based on the signs of the entries of the Iacobian of the system, which represents the signs of the partial derivatives of the state variables [63]. In this way, one can draw some conclusions on the possible dynamic behavior of the system based on structural properties of the interaction graph [64]. Here, we built for each identified interaction graph substructure a corresponding ODE model, that is, an ODE model whose underlying interaction graph reflects the respective minimal model structure (S5 Fig). A related approach enables the automatic conversion of logical models to ODE systems [65], but cannot be applied here as we utilize interaction graphs in our approach. Here, we built the ODE models from the interaction graphs in the following way: For the core model, each edge pA>pB from the interaction graph, where pA denotes the activated form of a species A and pB the activated form of a species B, gives rise to an ODE reaction describing the formation of pB. With mass action kinetics and assuming pA catalyzes the formation of pB but is not consumed, we get the reaction B->pB with kinetic rate law k1*pA*B, where k1 denotes a kinetic parameter for this reaction. In addition, we introduce the deactivation (for example, dephosphorylation) reaction pB>B with kinetic rate law k2*pB (S4 Table). An exception from this rule is the treatment of the species RSK. The interaction graph model describes activation of RSK as a two-step process: RSK->RSK_s->RSK_d. In this case, we assume that both activated forms (RSK_s and RSK_d) are deactivated to RSK (S4 Table, lines 26—29). For the Met receptor, an additional formation and degradation reaction was necessary (S4 Table, lines 3 and 4) to reflect the observed receptor dynamics. PDKl is assumed to be constitutively active [66, 67] and was thus included as a constant parameter in the ODE model. To incorporate the effect of inhibitor treatments into the model, inhibitor parameters for PDKl, PI3K, Met and MEK have been introduced. These parameters allow for a potential reduction of the kinetic parameters within a range of 0% to 100%. Furthermore, these parameters are coupled to binary switches, allowing them to be activated only in their respective condition.

We distinguish candidate mechanisms that provide an alternative, independent way of activation of the respective species and those that influence the core activation mechanism. An example for the first is the activation of PI3K by Ras, which is described as PI3K -> active PI3K with rate law k*Ras_active*P13K (S4 Table, line 41). An example for the latter is the effect of PAK on Raf, which is assumed to promote Raf activation through Ras, and which is thus included as Raf> pRaf with rate law k* Ras_active*pPAK* Raf in the ODE model (S4 Table, line 35). Translation of inhibiting edges of the interaction graph follows the same rules. As an example, deactivation of Raf by Akt is described as pRaf->Raf with rate law k* pAkt* pRaf (S4 Table, line 37).

Here, we built all model structures in a framework-specific format, which allows import and export to the common data format SBML [69]. The model 4_8_12 is available as a SBML in 81 Model.

To find the optimal parameter sets that describe the experimental data for each model structure, we performed parameter estimations. The framework is using a parallelised implementation of the CVODES ODE solver [70]. The procedure of parameter estimation is based on multiple local optimizations of different parameter starting values. For the optimizations, the LSQNONLIN algorithm (MATLAB, R2011a, The Mathworks Inc. (Natick, MA)) was used. Parameter values are limited between a range of 0.00001 and 1000, that is, eight orders of magnitude on a logarithmic scale. We assume that this parameter range is sufficient for our model selection approach. Parameter values close to upper or lower boundaries result from partial practical non-identif1ability of the model structures. This does not influence the resulting ranking of the model structures. Parameter estimation is based on multiple local optimizations of different parameter starting values. For the random sampling of the starting points, a latin hypercube method is utilized [71]. For each set of randomly selected initial parameter values, a local optimization procedure is performed. Here, parameter estimation of the software has been modified to cope with the high complexity of the parameter estimation problem: As an addition, parameter values that were estimated as being close to their respective upper or lower boundary were automatically fixed for a limited period of optimization cycles. As boundary-close parameters slow down the optimization procedure, temporal removal of these degrees of freedom for a short interval during the process led to an increase in fitting performance. To prove the validity of our optimization procedure, a Latin hypercube sampling approach with 1000 initial parameter sets has been performed (S7 Fig). As shown, the parameter estimation algorithm is able to converge into a global minimum within the parameter landscape in a reproducible manner. In addition to kinetic parameters, the ODE model is defined by scaling and noise parameters. These non-kinetic parameters were fitted in parallel to the kinetic parameters as described [68].

We assume that non-stimulated measurements of signalling components in our data set is sufficient for training the system to an initial steady state in an unstimulated setting. Additionally, our ODE model would be able to capture oscillatory behaviour if present.

This allows, in principle, for the occurrence of oscillatory behaviour of signalling components within the considered parameter space. However, as we did not observe oscillations in our experimental data, oscillations are not exhibited by the ODE model.

Performance of an estimated set of parameter values is characterized by its likelihood L. Here, we minimized the negative logarithm of the likelihood function (-Zlog(L)) during the optimization process. Therefore, the best performing model structure in combination with the estimated parameter set is described by the lowest -Zlog(L). To take into account different degrees of freedom of the models, two different methods are common in literature. First, the Akaike Information Criterion (AIC) is defined as Where k denotes the degrees of freedom in the respective model [72]. The second common method is known as the Likelihood-Ratio-Test (LRT), Where one model is defined as the null model and another nested model is compared against the null model [73]. Here, the comparison is performed by where Adf denotes the difference in degrees of freedom of the null model and the nested model. Hereby, one tests if a nested model is a valid simplification of the null model. The complete model is defined as the null model. Each of the tested model structures is therefore compared pairwise against the complete model. The result of each pairwise likelihood ratio test is then used to obtain the ranking of the corresponding model structures. In this work, we present all rankings with unprocessed -210g(L) values, AIC ranking and LRT for all model structures against the complete model. While the AIC allows the creation of a complete ranking and therefore a comparison of different candidate models against each other, the LRT provides us with more detailed information regarding the pairwise comparison of a nested model against the null model. In practice, the AIC slightly favours larger models due to the linear penalisation of the degrees of freedom of a model.

To compare the predictive power of the selected candidate model 4_8_12 and the complete model, we utilize prediction profiles as described [39]. For our analysis, prediction profiles have been calculated along the complete time course of the selected model 4_8_12 and the complete model for “active PI3K”. Through the calculation of prediction profiles, a range for the specified trajectories of the protein dynamic is given for each calculated time point, in which the likelihood value of the model stays within a 95% threshold.

The selected model structures contain either four or five candidate edges. We decided to consider random models with five candidate edges. Thus, we randomly selected subsets of size five from all candidate edges present in the complete model. The so-de-rived random models were compared with the minimal model structures: if a random model was identical to or comprised a minimal model structure, this model was rejected. In this way, we generated 50 random models and compared their performance with the minimal model structures (S6 Table and Fig 5D).

To study the effects of certain protein inhibitions, we performed in silico eXperiments. Here, we down-regulated all outgoing reactions for each protein, one protein at a time, by reducing the respective kinetic parameter to 50% of the original value and simulated the pathway dynamics under this perturbed condition. Furthermore, each possible combination of two protein inhibitions has been simulated. As readout, the area under curve of pAkt, pERK and the sum of pAkt and pERK have been evaluated for the measured time of 120 minutes and compared to the control condition. The corresponding heat maps are shown in Fig 8A).

Positive values mean that the double inhibition condition is having a stronger effect than the two single inhibitions, that is, there is a synergy effect. Negative values mean that the single inhibitions are more effective than the combination. A neutral value of 0 means the two inhibitors work in an additive manner.

Supporting Information

Species in the interaction graph model. The table describes the species present in the interaction graph model. The species are considered in their active form, for example, Akt is considered as Akt phosphorylated on serine 473. Species 11: MEK1/2 phosphorylated at serine 217/221 Species 12: MEK1/2 phosphorylated at serine 298Species 12: MEK1/2 phosphory-lated at threonine 292 Species 19, 21, 22: considered in their active GTP-bound form. Species 25, 26: RSK_s refers to p9ORSK phosphorylated on a single serine residue. RSK_d refers to p9ORSK phosphorylated on two serine residues and considered as active p9ORSK. (DOCX)

Mechanisms in the interaction graph model (core model). The table describes the list of reactions of the HGF interaction graph (core model). (DOCX)

Candidate mechanisms in the interaction graph model. The table describes the list of reactions defined as candidate mechanisms in the HGF interaction graph master model. Each number corresponds to a candidate mechanism. The letters indicate the reactions for the candidate mechanisms composed by more than one reaction. The asterisks indicate the edges included in more than one candidate mechanism. (DOCX)

Reactions in the ODE models. The table describes the list of reactions of the ODE models. The first list describes the reactions of the core model (reaction 1—29); the second list describes the reactions of the candidate mechanisms (reaction 30—42). In the second column, each reaction is shown in a schematic representation; in the third column, the respective kinetic rate law is shown. Values of each kinetic parameter for model 4_8_12 are shown in SS Table. Parameters “Met_inh”, “PDK_inh” and “MEK_inh” represent binary values dependent on the

List of parameter names and values of the model 4_8_12. The table describes the list of parameter names and values of the ODE model 4_8_12. The model was calibrated on 2200 data points and 25 experimental conditions. The second column shows the loglo value of each kinetic parameter involved in the model reactions shown in S4 Table. The X indicates that the parameter is not present in the model 4_8_12. (DOCX)

List of models and candidate mechanisms. The table describes the list of candidate mechanisms for each selected minimal model structure (gray), model combinations (gray) and random models (rand): 1 indicates that the candidate mechanism is present, 0 indicates that it is not present. The models are sorted according to the Likelihood ratio test value as 81 Fig. Quantitative immunoblotting. Quantitative immunoblotting of primary mouse hepa-tocytes treated with MEK inhibitor (U0126), PDKl inhibitor (BX-912), their combination or DMSO (-) and subsequently stimulated with 40 ng/ml HGF for the indicated time points. A) Immunoprecipitation (IP) of total SOS protein. The band shift indicates the phosphorylated inactive form of SOS (higher band). B) Total cellular lysate was used to detect phosphorylation of MEK1/2, ERK1/2 and Akt on serine 473 and total actin used as normalizer. Quantification of the corresponding immunoblots is shown for pAkt, pMEK, pERK. C) Total cellular lysate was used to detect p9ORSK phosphorylation on serine 227 and 380, total p9ORSK and protein disulfide-isomerase (PDI) used as normalizer. Quantification of the corresponding immuno-blots is shown for p9ORSK Ser380 and p9ORSK Ser227. Samples were loaded in a random order to reduce the correlated error [60]. Experiments shown in panels B and C were performed at least in triplicates and one representative dataset is shown. (EPS)

A) Quantitative immunoblotting of primary mouse hepatocytes treated with siRNA targeting ERK1/2 and Akt and subsequently stimulated with 40 ng/ml HGF for the indicated time points. Total cellular lysate was used to detect MEK1/2 phosphorylation, total ERK, total MEK1/2 and total actin used as normalizer. Samples of hepatocytes treated with siRNA targeting Akt were processed with the samples of siRNA targeting ERK1/2. Data quantification is shown in S3B Fig. Data of Akt siRNA treatment has not been used for model calibration. Samples were loaded in a random order to reduce the correlated error [60]. B) Knockdown efficiency of siRNA targeting ERK1/2. Signal intensity of total ERK1/2 and actin was measured by LumiAnalysis software (Roche Diagnostic). Total ERK1/2 signal was divided by the respective actin signal and the average of 6 replicas is shown with the corresponding standard deviation. The knockdown efficiency for ERK1/2 siRNA is indicated as percentage of the control condition. (EPS)

For each protein, the fold change of the phosphorylation state measured by quantitative immunoblotting of two different experimental conditions is shown on a logarithmic scale at the indicated time points after 40 ng/ml of HGF stimulation. A) Fold change of Akt phosphorylation on threonine 308 and the active, not phosphorylated, SOS is shown (representative blot shown in SlA Fig) B) Fold change of the phosphorylation state of the indicated proteins upon PI3K inhibitors treatment. These data resulted to be not conclusive, and, therefore, it has not been discretized and it was not used for model calibration. Each row refers to one experiment; same experimental conditions are grouped with magenta lines. (EPS)

Shown are the predictions from all selected minimal model structures (Fig 4) and their comparison with the discretized data shown in Fig 3B. “Early” refers to the first response following HGF stimulation. For later time points, either the first effect that is opposite the initial response is given, or an arrow of same direction indicates that the qualitative response was equal for all time points. Arrow pointing down/ up: the inhibition can only cause a decreased/ increased activation of the measured protein. Bullet point: the inhibition does not affect the measured protein. Combined up/ down arrow and bullet point: the model does not restrict the response to the inhibition. (EPS)

A) Shown is the interaction graph underlying the complete ODE model, thus representing the sign-structure of its Jacobian matrix. Black arrows indicate positive, red arrows negative influences. Note that PDKl receives no input and is not consumed in the ODE model. For simplification, it was therefore treated like a parameter, and the edges indicate its influences (partial derivatives) on other components. The same interaction graph would follow if we included PDKl as a (constant) state variable in the model. Furthermore, in addition to the nodes that are contained in the interaction graph given in Fig 4, which represent the activated forms of the proteins, the underlying interaction graph of the ODE model contains for each protein one node representing its inactive form. With the chosen mass action kinetics (S4 Table), the interaction graph contains for each reaction where the activated form of protein A, pA, positively influences the activation of protein B (B —> pB), the positive edge pA —> pB and the negative edge pA —> B, the latter thus giving rise to an additional indirect negative influence of pA on pB (pA —>B —> pB). In general, these indirect effects can become visible in the transient dynamics. However, as in our ODE model formulation the active and inactive form of a protein are conserved moieties, the interaction graph of the Iacobian can be reduced to the graph given in panel (B). B) This graph contains the same set of nodes as the interaction graph given in Fig 4, except for an additional node representing the inactive form of the Met receptor: the moieties of Met and pMet are not conserved as additional uptake and degradation is considered (S4 Table). Still, the paths downstream the activated receptor are the same as in the interaction graph in Fig 4. Compared to panel (A), all indirect effects except one could be eliminated: as activation of RSK is modeled as a two-step process, the interaction graph of the ODE model contains with PDK —> pRSK_s and pRSK_d —> pRSK_d two negative edges that are not contained in the interaction graph in Fig 4. However, even if these edges are considered, the selected substructures of the complete interaction graph model are the minimal structures needed to eXplain all qualitative effects from the eXperimental data. (EPS)

A) Rankings represent the forward selection of selected minimal model structures; B) backward selection, where the building blocks are removed from the complete model; C) the model combinations selection. D) Rankings of model selection including minimal model structures, model combinations and random models. Random model structures consist of randomly selected combinations of candidate edges (S6 Table). Approximately 50% of the random models performed worse than the worst performing selected minimal model structure. Well performing random model structures are closely related to well performing selected minimal model structures and their combination (S6 Table). The Akaike Information Criterion (AIC) has been utilized to penalize the likelihood. All rankings of model selection present the negative logarithmic likelihood penalized by parameter difference on the y-aXis. Model identifiers are

A) Mul-tistart optimization has been performed With 1000 initial parameter sets. Fits have been sorted according to the resulting likelihood. B) Detailed View of panel (A) Where the trend of the likelihood of 780 fits is shown.

A) Structure of the core model, the complete model (D), model 4_8_12 (G) and corresponding plots showing representative model trajectories (solid lines) of the phosphorylation kinetic of MEK (B-E-H) and Akt (C-F-I) measured by quantitative immunoblotting in primary mouse hepatocytes pretreated with the indicated inhibitors and stimulated with 40 ng/ml HGF for the indicated time (stars). Values of the likelihood ratio test with the threshold of 95% are shown under each corresponding model. In the plots, y-axes show the concentration of the respective measured protein in arbitrary units on a logarithmic scale. The colored area surrounding the model trajectory represents the confidence interval delimited by the dashed line. Treatments are color-coded as indicated in the figure. The core model (A) does not capture the difference between the experimental conditions while the complete model (D) and model 4_8_12 (G) clearly distinguish them. For example, the complete model and model 4_8_12 capture the effect of MEK inhibitor and PDKl inhibitor on pMEK1/2 showing a delayed and sustained response. (EPS)

Quantitative immunoblot-ting of primary mouse hepatocytes treated With Akt inhibitor VIII or DMSO and subsequently stimulated with 40 ng/ml HGF for the indicated time points. A) Total cellular lysate or immu-noprecipitation (1P) was used to detect the indicated proteins. Rafl protein has a distinguished phosphorylation pattern: Phosphorylation at serine 338 contributes to Rafl activation, while phosphorylation at serine 259 is an inactivation signal and targeted by Akt (shown in Fig 7). B) Immunoblot data were quantified by ImageQuant TL version 7.0 software (GE Healthcare). (EPS)

Quantitative immunoblot-ting of primary mouse hepatocytes treated with Akt inhibitor VIII or DMSO and subsequently stimulated with 40 ng/ml HGF for the indicated time points. A) Total cellular lysate or immu-noprecipitation (1P) was used to detect the indicated proteins. Rafl protein has a distinguished phosphorylation pattern: Phosphorylation at serine 338 contributes to Rafl activation, while phosphorylation at serine 259 is an inactivation signal and targeted by Akt. B) Immunoblot data were quantified by ImageQuant TL version 7.0 software (GE Healthcare).

Quantitative immunoblot-ting of primary mouse hepatocytes treated with Akt inhibitor VIII or DMSO and subsequently stimulated with 40 ng/ml HGF for the indicated time points. A) Total cellular lysate or immu-noprecipitation (1P) was used to detect the indicated proteins. Rafl protein has a distinguished phosphorylation pattern: Phosphorylation at serine 338 contributes to Rafl activation, while phosphorylation at serine 259 is an inactivation signal and targeted by Akt. B) Immunoblot data were quantified by ImageQuant TL version 7.0 software (GE Healthcare).

Quantitative immunoblotting of primary mouse hepatocytes treated with the inhibitor combinations and 0.1uM Met inhibitor PHA665752 and 10 microM PDKl inhibitor (BX-912), 0.1 microM P3K inhibitor (PI-103) and 0.1 microM MEK inhibitor (U0126) or DMSO (-) and subsequently stimulated with 40 ng/ml HGF for the indicated time points. A) Total cellular lysate was used to detect phosphory-lation of ERK1/2 and Akt on serine 473 and total actin used as normalizer. The quantification plots are shown in Fig 8D. B) Replica of the experimental validation. Total cellular lysate was used to detect phosphorylation of ERK1/2 and Akt on serine 473 and total actin used as normalizer. The plots show the quantification of the corresponding immunoblots.

Pseudocodes. Pseudocodes for the selection of minimal model structures from the interaction graph master model. (PDF)

ODE model 4_8_12 formatted in SBML. The 4_8_12 model contains the core edges and all candidate edges deriving from the combination of models 4_8_12. (XML)

pERK_au: sum of ERK1 and ERK2 phosphorylation. Units: arbitrary units (au). pMEK_au: sum of MEKl and MEK2 phosphorylation. Units: arbitrary units (au). sin-gle_pRSK_au: p90RSK phosphorylation corresponding to the p90RSK_s of the model. Units: arbitrary units (au). double_pRSK_au: p90RSK phosphorylation corresponding to the p90RSK_d of the model. Units: arbitrary units (au). Actin_au: actin. Units: arbitrary units (au). PDI_au: Protein disulf1de isomerase. Units: arbitrary units (au). EXp_Gel_ori: experiment identifier is indicated in letters (i.e. “A_tc”); reference number of gel. Gel_ori: internal reference gel number. Exp: internal reference eXperiment identifier. NaN: no datapoint has been measured for the indicated condition. Columns defining perturbation treatment (inhibitor or siRNA): 0 indicates that the corresponding inhibitor or siRNA was not present; 1 indicates that the corre-82 Dataset. ExperimentalData_Ras_ProteinArray. The datasets contain the experimental data used for the selection of minimal model structures and model calibration. The datafile is characterized by a header row as follow: Time: duration of HGF stimulation. Units: minutes. Gel: unique number to identify the gel for each analyzed protein. In this special case, “gel” refers to protein array. HGF_input: concentration of HGF used for stimulation. Units: ng/ml Met_ihn: Met inhibitor Mek_inh: MEK inhibitor Pdk_inh: PDKl inhibitor ERK siRNA: siRNA targeting ERK1/2 RasneW_au: active Ras measured by protein array. Units: arbitrary units (au) NaN: no datapoint has been measured for the indicated condition. Columns defining perturbation treatment (inhibitor or siRNA): 0 indicates that the corresponding inhibitor or siRNA was not present; 1 indicates that the corresponding inhibitor or siRNA was present. (XLSX)

ExperimentalData_Ras_Immunoblot. The datasets contain the experimental data used for the selection of minimal model structures and model calibration. The datafile is characterized by a header row as follow: Time: duration of HGF stimulation. Units: minutes. Gel: unique number to identify the gel for each analyzed protein HGF_input: concentration of HGF used for stimulation. Units: ng/ml Met_ihn: Met inhibitor Mek_inh: MEK inhibitor Pdk_inh: PDKl inhibitor ERK siRNA: siRNA targeting ERK1/2 aRas_au: active Ras measured by immunoblot. Units: arbitrary units (au) EXp_Gel_ori: experiment identifier is indicated in letters (i.e. “A_tc”); reference number of gel. Gel_ori: internal reference gel number Exp: internal reference experiment identifier. (XLSX)

ExperimentalData_ERK phosphorylation_MassSpec. The datasets contain the experimental data used for the selection of minimal model structures and model calibration. The datafile is characterized by a header row as follow: Time: duration of HGF stimulation. Units: minutes. Gel: unique number to identify the gel for each analyzed protein HGF_input: concentration of HGF used for stimulation. Units: ng/ml Met_ihn: Met inhibitor Mek_inh: MEK inhibitor Pdk_inh: PDKI inhibitor ERK siRNA: siRNA targeting ERK1/2 pERK_au: sum of ERK1 and ERK2 phosphorylation measured by mass spectrometry. Units: arbitrary units (au) EXp_Gel_ori: eXperiment identifier is indicated in letters (i.e. “A_tc”); reference number of gel. Gel_ori: internal reference gel number EXp: internal reference eXperiment identifier. (XLSX)

ExperimentalData_Akt phosphorylation_MassSpec. The datasets contain the eXperimental data used for the selection of minimal model structures and model calibration. The datafile is characterized by a header row as follow: Time: duration of HGF stimulation. Units: minutes. Gel: unique number to identify the gel for each analyzed protein HGF_input: concentration of HGF used for stimulation. Units: ng/ml Met_ihn: Met inhibitor Mek_inh: MEK

Units: arbitrary units (au) Akt_rel_obs: unphosphorylated Akt measured by mass spectrometry. Units: arbitrary units (au) EXp_Gel_ori: experiment identifier is indicated in letters (i.e. “A_tc”); reference number of gel.Gel_ori: internal reference gel number EXp: internal reference eXperiment identifier.

Acknowledgments

We thank Dr. Stephan Feller for providing the Raf-RBD construct, Dr. Bettina Hahn for the mass spectrometry measurements and Dr. Clemens Kreutz for the support for the analysis of the predictive power of the models.

Author Contributions

Performed the experiments: LAD SB KW AK RM. Analyzed the data: LAD RS TM SHR NI MS. Contributed reagents/ma-terials/analysis tools: AR. Wrote the paper: LAD RS TM UK IT SK.