Signaling pathways are characterized by crosstalk, feedback and feedforward mechanisms giving rise to highly complex and cell-context specific signaling networks. 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 phosphoinositide-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.
Cells receive extracellular signals and process them through intracellular signaling pathways to regulate cellular responses. 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, non-linear intracellular signaling networks, whose topologies are often context-specific and altered in diseases [1]. An important factor that contributes to liver regeneration and has been implicated in the context of resistance to targeted tumor therapy is hepatocyte growth factor (HGF). 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. In general, activation of PI3K leads to the generation of phosphatidylinositol 3,4, 5–triphosphate (PI3,4, 5-P3) that serves as docking site for the serine/threonine protein kinase Akt at the plasma membrane. 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 (MEK1 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 p90RSK targeting SOS [4], and a positive feedback loop operates from ERK to Raf via RKIP, a protein playing a dual role as positive and negative regulator of MAPK signaling [5]. 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. Knowledge of cell-context specific feedback and crosstalk mechanisms is central to gain insights into mechanisms that cause undesired effects of targeted therapy. 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. To this aim, mathematical models provide unique tools to disentangle complexity and to predict the impact of perturbations. 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. Several computational methods have been developed to infer and analyze signaling networks. 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-suited 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. To analyze the impact of crosstalk and feedback regulation, dynamic modeling approaches using coupled ordinary differential equations (ODEs) are most suited and allow quantitative insights [24,35–38]. 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. Here, we present a novel hybrid approach (Fig 1), which combines qualitative and quantitative modeling techniques to unravel the HGF induced activation of MAPK and PI3K signaling in primary mouse hepatocytes based on time-resolved experimental data. 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.
To unravel the crosstalk between the MAPK and PI3K signaling pathways upon HGF stimulation in hepatocytes, we built an interaction graph model of HGF induced activation of these pathways based on literature knowledge (Fig 2). 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 (S1–S3 Tables). The full graph including the core model and all candidate mechanisms is considered as a non-cell-type specific “HGF master model”. To analyze which of the candidate mechanisms are active in our cellular system, we treated primary mouse hepatocytes with HGF in combination with specific Met inhibitor (PHA 665752), MEK inhibitor (U0126), PI3K inhibitors (LY294002, Wortmannin and PI-103) and PDK1 inhibitor (BX-912) alone or in combination. 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 p90RSK (S1 and S2 Fig). p90RSK activation occurs in two-steps: RSK_s denotes p90RSK phosphorylated at a single residue, RSK_d refers to the double phosphorylated, fully active form of p90RSK (S1 and S2 Tables). We also measured the dynamics of SOS1 activation upon MEK inhibitor treatment (S1A 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 SOS1 activation, the band shift was quantified as an indicator for its activation status (S1A and S3A 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 SOS1 activation resulting in a more sustained activation (S1A and S3A 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 PDK1 inhibitor reduced Akt phosphorylation on serine 473 and threonine 308 (Fig 3A and S3A Fig). Interestingly, a decreased phosphorylation of MEK, ERK, RSK_s and RSK_d upon PDK1 inhibitor treatment was observed at early time points and a subsequent increase at later time points. When applying MEK and PDK1 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 PDK1 inhibitors or with PDK1 inhibitor alone except for Akt phosphorylation, which was additionally decreased upon the combinatorial inhibitor treatment. The comparison between the combinatorial treatment with MEK and PDK1 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 (S3B Fig). 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." Early" refers to the initial qualitative response 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. Given an interaction graph model, we can predict the possible qualitative responses of the considered proteins for the given experimental conditions using the concept of the dependency matrix [33]. 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. Thereby, all possible minimal model structures were selected, each being composed of the core model and a subset of the candidate mechanisms. 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. To develop corresponding ODE models with an appropriate number of parameters, we compressed all selected minimal model structures by removing parallel mechanisms with the same sign and by compressing linear chains [30]. 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 SOS1, either directly or indirectly via RSK_d, (iii) a positive route from ERK to MEK through various mechanisms and (iv) a positive route from PDK1 to MEK. The latter route is included as a direct edge from PDK1 to MEK in all models but model 8, whereas model 8 contains a longer path via RSK_d, Ras and Akt. 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 immunoblotting and protein array, respectively, were considered. In addition, the degree of phosphorylation 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 SOS1 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. To challenge our approach, we generated 50 model structures consisting of randomly selected combinations of candidate edges. 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. At first glance, model 4_8_12 includes three feedback loops (Fig 6A). Within the MAPK pathway, a positive and a negative feedback emerge from the candidate mechanisms ERK to Raf1 and RSK_d to SOS1, respectively. Within the PI3K pathway, the edge PI3K to Gab1 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 S8 Fig. To experimentally validate model 4_8_12, we focused on the identified interaction from Akt to Raf1 and predicted the impact of different degrees of Akt inhibition on the inhibitory impact of Akt on Raf1. 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 Raf1, which cannot be directly compared to Raf1 phosphorylated at a specific phosphorylation site as phosphorylation on serine 338 contributes to Raf1 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 Raf1 as a proxy to experimentally address the inhibitory impact of Akt on Raf1. Experimentally, we stimulated primary mouse hepatocytes with HGF in the absence and presence of the specific Akt inhibitor VIII and monitored the impact on Raf1 phosphorylation on serine 259. As shown in Fig 7C and S9–S11 Figs, we achieved an inhibition of Akt phosphorylation between 90% and 100% and concomitantly observed a decrease of Raf1 phosphorylation on serine 259 (Fig 7D and S9–S11 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. The candidate mechanisms within model 4_8_12 generate crosstalk as well as feedforward and feedback loops within the network structure leading to robust network behavior. 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. To identify strategies to efficiently inhibit Akt and ERK signaling, we used the inferred 4_8_12 model to determine the effect of reducing the rate of outgoing reactions of each protein by 50%. 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. To experimentally validate the model predictions on combinatorial treatments, we first performed single inhibitor treatments targeting PDK1, PI3K, Met and MEK prior to HGF stimulation and estimated the inhibition strength parameters for each individual inhibitor for model 4_8_12 (Fig 8C). 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 S12 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.
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]. Our network inference approach combines qualitative (interaction graph) and quantitative (ODE) modeling and, hence, benefits from the strengths of both methods. 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 pre-selecting 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 non-linear effects such as saturation. Applying our approach to HGF signaling in primary mouse hepatocytes, the ranking of the resulting minimal ODE model structures by forward selection indicated that (i) model structures harboring the negative feedback edge from RSK_d to SOS1 in most cases outperformed their partner model containing the negative feedback edge from ERK to SOS1 and (ii) random model structures, which do not include any of these negative feedback edges, are strongly outperformed by the candidate models. 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 SOS1 could also be mediated by p90RSK [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 Raf1 from RKIP-Raf1 complex resulting in MEK activation. This represents an additional mechanism to fine-tune signaling responses within the negative feedback from downstream proteins to SOS1. Furthermore, our approach unraveled a negative interaction from Akt to Raf1, which we confirmed experimentally in our cellular model system. Raf1 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 Raf1; the positive feedback loop is mediated by PI3K activated Gab1 that strengths the MAPK activation via Rac and PAK. Thus, our proposed network structure shows positive and negative feedback loops that are highly interlinked. Feedback loops are the base of non-trivial 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]. Components of the MAPK and PI3K signaling pathways are frequently altered in pathological conditions and, thus, there is much interest in developing targeted therapies [56,57]. 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. Our modeling strategy delivers as output a dynamic ODE model validated with respect to both parameters and network structure, as an invaluable tool to not only elucidate cellular signaling pathways, but also to design cell-type or even patient-specific intervention strategies. 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 PDK1 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.
All animal experiments were approved by the governmental review committee on animal care of the state Baden Württemberg, Germany (reference number A24/10). Anesthesia was carried out by intraperitoneal injection of 5 mg ketamine hydrochloride 10% (w/v) (Bayer Health Care, Leverkusen, Germany) per 100 mg body weight and 1 mg xylazine hydrochloride 2% (w/v) (Pfizer, Berlin, Germany) per 100 mg body weight. Mice of 8–12 week-old C57BL/6N (Charles River, Sulzfeld, Germany) were used for primary hepatocyte isolation. 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 μM dexamethasone, 10 μg/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% CO2 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. For time course and dose response experiments after isolation, hepatocytes were seeded at confluence (2x106 cells/6 cm dish) and cultivated as described before. For mass spectrometry analysis, after isolation, hepatocytes were seeded at confluence (5x106 cells/10cm dish) and treated as described [59]. After 6 hours in serum-free cultivation medium depleted of dexamethasone, hepatocytes were stimulated with indicated concentration of recombinant mouse HGF (HGF) (R&D Systems) at indicated time points (Supporting Information Datasets). For inhibitor experiments, cells were treated for 30 to 45 minutes with 20 microM MEK inhibitor U0126,5–25 microM PI3K inhibitors LY294002, Wortmannin (all Cell Signaling Technologies), or PI-103 (Calbiochem), 20 microM PDK1 inhibitor BX-912 (Axon Medchem), 20 microM Akt inhibitor VIII (Millipore), or 5–20 microM Met inhibitor PHA-665752 (TOCRIS Bioscience) prior to stimulation with 40 ng/ml of HGF for the designated time points. The control condition was carried out by treating the cells with dimethyl sulfoxide (DMSO, Sigma-Aldrich). Model validation experiments were carried out by stimulating the primary mouse hepatocytes with the following inhibitor combinations: 1 microM MEK inhibitor U0126 (Cell Signaling Technologies) with 0. 1 microM PI3K inhibitor PI-103 (Calbiochem); 5 microM MEK inhibitor U0126 (Cell Signaling Technologies) with 0. 5 microM PI3K inhibitor PI-103 (Calbiochem); 10 microM PDK1 inhibitor BX-912 (Axon Medchem) with 5 microM Met inhibitor PHA-665752 (TOCRIS Bioscience); 10 microM PDK1 inhibitor BX-912 (Axon Medchem) with 1 microM Met inhibitor PHA-665752 (TOCRIS Bioscience); 10 microM PDK1 inhibitor BX-912 (Axon Medchem) with 0. 1 microM Met inhibitor PHA-665752 (TOCRIS Bioscience). 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. For siRNA treatment after isolation, hepatocytes were seeded at subconfluence (900000 cells/well in 6 wells/plate). 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 siRNAs 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. Cells were lysed on ice at the indicated time points using total cell lysis buffer (20 mM Tris pH7. 4,150 mM NaCl, 1mM EDTA pH8. 0 (AppliChem), 1% (v/v) NP40 (Roche Applied Sciences), 1mM ZnCl2,1mM MgCl2,1 mM Na3VO4,10 mM NaF, 10% glycerol, 0. 1 mg/ml AEBSF, 1μg/ml aprotinin). Total cellular lysate was used to assay MEK, ERK, Akt, Raf1, or they were subjected to immunoprecipitation (IP) using 500 micrograms of total protein and subsequently assayed on SDS-PAGE. For IP anti-Met, aanti-Raf (both Santa Cruz Biotechnology) and anti-SOS (Merk Millipore) antibodies were used. Samples were loaded in a random order to reduce the correlated error [60]. Primary antibodies anti-p44/42 threonine 202/tyrosine 204 phosphorylation (anti-pERK1/2), anti-p44/42 total protein (anti-ERK), anti-MEK serine 217/221 phosphorylation, anti-MEK total protein, anti-Akt serine 473 phosphorylation, anti-Akt threonine 308 phosphorylation, anti-Akt total protein, anti-p90RSK serine 227 phosphorylation, anti-p90RSK 380 phosphorylation, anti-p90RSK total protein (Cell Signaling Technology), anti-Raf serine 259 phosphorylation (Cell Signaling Technology), anti-Raf serine 338 phosphorylation (Upstate), anti-Raf total protein (Santa Cruz Biotechnology), anti-SOS (Merk Millipore), anti-phosphotyrosine 4G10 (Merck Millipore), anti-PDI (Upstate), anti-actin (Sigma) were used. For the quantification of active Ras, cells were lysed with total cell lysis buffer (125 mM HEPES (pH 7. 5), 750 mM NaCl, 1% NP40,50 mM MgCl2,10% glycerol, 0. 1 mg/ml AEBSF, 1 microgram/ml aprotinin and protease inhibitor cocktail (Roche) ). The whole lysate was subjected to IP supplemented with MgCl2 (50 mM) and Raf-RBD-GST-beads or processed on protein array (as described in the protein array section). After the IP, the lysates were assayed on SDS-PAGE and primary antibody anti-panRas (OP21) (Calbiochem) was used for detection. Secondary horseradish peroxidase-coupled antibodies (anti-rabbit HRP and anti-mouse HRP) were used (Amersham Biosciences). Target proteins were visualized using enhanced chemiluminescence (GE Healthcare), and signals were acquired using the CCD camera-based device LumiImager (Roche Diagnostic) or ImageQuant LAS 4000 biomolecular imager (GE healthcare). Immunoblot data were quantified using LumiAnalysis software (Roche Diagnostic) or ImageQuant TL version 7. 0 software (GE Healthcare). Primary mouse hepatocytes were isolated and cultivated as described above. For Akt and ERK phosphorylation analysis, after 6 hours of serum-free cultivation medium depleted of dexamethasone, 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 ERK1 and Akt (Santa Cruz Biotechnology). For IP, anti-ERK1 (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 manufacturer’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 analysis was performed for the simultaneous detection of Akt and ERK phosphorylation. 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 phosphorylated 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. Recombinant proteins were used as calibrators. 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-inhibitor 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 ddH2O and subsequently dried. Slides were scanned with LiCor Odyssey scanner at 700 and 800 nm (Resolution 21 μm, High Quality, Intensity: 5/4; Offset: 0, size 7x5 Scanning time approx. 35min/2slides). The quantification of the signal was performed as described [61]. We generated protein array analysis to quantify active Ras by using a mouse antibody (Pan-Ras, OP40, Calbiochem) and GST-tagged Ras-binding-domain (RBD) of Raf (GST-Raf-RBD 1–149) [56] for Ras detection. To generate the pGex plasmid harboring the GST-tagged Raf-RBD we cloned the Raf-RBD into the pGex vector via BamHI and EcoRI sites. The pGex-GST-Raf-RBD was expressed in BL21 bacteria and the fusion protein was purified using Gluthatione Sepharose beads (GE Healthcare). 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 gammaS-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 are an abstract representation of signaling pathways [32]. The nodes in these graphs represent the signaling molecules, in our application usually the activated forms of the proteins (S1 Table). Directed positive and negative edges represent activating or inhibiting influences between the proteins.
Cellular responses to extracellular stimuli are driven by activation of intracellular signaling pathways. 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.
Disentangling the Complexity of HGF Signaling by Combining Qualitative and Quantitative Modeling