Advertisement

  • Loading metrics

Open Access

Peer-reviewed

Research Article

Context-specific network modeling identifies new crosstalk in β-adrenergic cardiac hypertrophy

Abstract

Cardiac hypertrophy is a context-dependent phenomenon wherein a myriad of biochemical and biomechanical factors regulate myocardial growth through a complex large-scale signaling network. Although numerous studies have investigated hypertrophic signaling pathways, less is known about hypertrophy signaling as a whole network and how this network acts in a context-dependent manner. Here, we developed a systematic approach, CLASSED (Context-specific Logic-bASed Signaling nEtwork Development), to revise a large-scale signaling model based on context-specific data and identify main reactions and new crosstalks regulating context-specific response. CLASSED involves four sequential stages with an automated validation module as a core which builds a logic-based ODE model from the interaction graph and outputs the model validation percent. The context-specific model is developed by estimation of default parameters, classified qualitative validation, hybrid Morris-Sobol global sensitivity analysis, and discovery of missing context-dependent crosstalks. Applying this pipeline to our prior-knowledge hypertrophy network with context-specific data revealed key signaling reactions which distinctly regulate cell response to isoproterenol, phenylephrine, angiotensin II and stretch. Furthermore, with CLASSED we developed a context-specific model of β-adrenergic cardiac hypertrophy. The model predicted new crosstalks between calcium/calmodulin-dependent pathways and upstream signaling of Ras in the ISO-specific context. Experiments in cardiomyocytes validated the model’s predictions on the role of CaMKII-Gβγ and CaN-Gβγ interactions in mediating hypertrophic signals in ISO-specific context and revealed a difference in the phosphorylation magnitude and translocation of ERK1/2 between cardiac myocytes and fibroblasts. CLASSED is a systematic approach for developing context-specific large-scale signaling networks, yielding insights into new-found crosstalks in β-adrenergic cardiac hypertrophy.

Author summary

Pathological cardiac hypertrophy is a disease in which the heart grows abnormally in response to different motivators such as high blood pressure or variations in hormones and growth factors. The shape of the heart after its growth depends on the context in which it grows. Since cell signaling in the cardiac cells plays a key role in the determination of heart shape, a thorough understanding of cardiac cells signaling in each context enlightens the mechanisms which control response of cardiac cells. However, cell signaling in cardiac hypertrophy comprises a complex web of pathways with numerous interactions, and predicting how these interactions control the hypertrophic signal in each context is not achievable by only experiments or general computational models. To address this need, we developed an approach to bring together the experimental data of each context with a signaling network curated from literature to identify the main players of cardiac cells response in each context and attain the context-specific models of cardiac hypertrophy. By utilizing our approach, we identified the main regulators of cardiac hypertrophy in four important contexts. We developed a network model of β-adrenergic cardiac hypertrophy, and predicted and validated new interactions that regulate cardiac cells response in this context.

Citation: Khalilimeybodi A, Paap AM, Christiansen SLM, Saucerman JJ (2020) Context-specific network modeling identifies new crosstalk in β-adrenergic cardiac hypertrophy. PLoS Comput Biol 16(12): e1008490. https://doi.org/10.1371/journal.pcbi.1008490

Editor: Martin Meier-Schellersheim, National Institutes of Health, UNITED STATES

Received: August 24, 2020; Accepted: November 5, 2020; Published: December 18, 2020

Copyright: © 2020 Khalilimeybodi et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

Data Availability: All relevant data are within the manuscript and its Supporting Information files.All images obtained from western blot and immunofluorescence methods and used for making graphs are accessible via our data repository https://doi.org/10.6084/m9.figshare.12824162.v1.

Funding: This study was funded by the awards received by J.J.S. from National Institutes of Health (HL137755, HL137100, HL127624) at https://www.nih.gov/ and the National Science Foundation (1252854) at https://www.nsf.gov/ (to J.J.S.). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Competing interests: The authors have declared that no competing interests exist.

Introduction

Cardiac hypertrophy, a putative risk factor of heart failure, is a context-dependent phenomenon [1–3] which is regulated by a broad range of biochemical and biomechanical factors interacting in a complex network of signaling pathways [4,5]. Homogenous growth of the heart occurs in the contexts of exercise and pregnancy [6,7]. But pathological contexts such as hypertension lead to an irreversible thickening of the heart wall [8–10], distinct from heart chamber dilation which is triggered in the contexts of myocardium infarction and valvular heart disease [11]. Identifying signaling pathways specific to each context has long been of interest [12], but most experimental studies have focused on single stimuli and isolated pathways for practicality and limited resources [13–15]. While computational models of cardiac hypertrophy signaling [16–21] address this issue to a certain extent by considering multiple pathways, predicting interactions between pathways in each context has remained a challenge due to the size and connectivity of cardiac signaling network.

Recent advances in experimental and computational approaches [22–27] facilitate the investigation of whole signaling networks and inference of new causal interactions. Among the established approaches to derive a network model [28–30], data-driven approaches that do not use prior knowledge has been widely applied for modeling gene regulatory and metabolic networks [31,32]. But modeling large and complex signaling networks via this method is not feasible owing to its low performance in model scale-up and interpretability [33]. In contrast, prior knowledge approaches have been implemented successfully by several groups to model large-scale signaling networks [34,35]. Yet prior knowledge models comprise numerous parameters to train, and only limited experimental data are available for cell signaling networks. Accordingly, approaches with less required data, e.g. Boolean and Logic-based ODE (LDE), are prevalent for modeling large-scale signaling networks [36–39]. The other challenge in implementing a prior knowledge network model is its inability to predict signaling outside the pre-defined scope of the network. Hence, recently, researchers prefer to utilize hybrid approaches for modeling large-scale networks [40,41].

We previously developed a large-scale signaling network model of cardiac hypertrophy using logic-based differential equations [4]. This model was built from a wealth of prior knowledge and was validated against in vitro and in vivo rat data [4,42]. However, this model does not account for the network specificity in each context and needs revision to accurately predict myocyte response to inhibitory perturbations [4,23]. Moreover, emerging context-specific and patient-specific high-throughput data increases the need for systematic approaches to develop such network models.

In this study, we developed a systematic method to derive a hybrid context-dependent model for a signaling network, starting from its interaction graph and context-specific data. Employing the interaction graph of our previous hypertrophy network and context-specific data curated from the literature, we identified the key signaling reactions that regulate the cardiomyocyte response in four different contexts (hypertrophic agonists). More importantly, in developing a context-specific model for β-adrenergic-induced hypertrophy, we predicted and then experimentally validated previously unknown crosstalks between Ca2+-dependent and β-adrenergic non-classical pathways, CaMKII-Gβγ and CaN-Gβγ, which regulate hypertrophy network response in the context of β-adrenergic stimulation. These findings indicate the strength of the CLASSED approach in systematically predicting unknown context-specific network connections to provide a better understanding of a large-scale network response in different contexts.

Materials and methods

Ethics statement

All experimental procedures were conducted in compliance with the Guide for the Care and Use of Laboratory Animals published by NIH and approved by the University of Virginia Institutional Animal Care and Use Committee.

Research pipeline

Fig 1 illustrates the pipeline of the CLASSED (Context-specific LArge-Scale Signaling nEtwork Development) approach to develop context-specific large-scale signaling network models. This approach employs 1) the interaction graph of a prior-knowledge network, and 2) the context-specific perturbation data, in predefined Microsoft Excel formats as inputs. Automated validation of the model is conducted by converting the interaction graph to a logic-based ODE (LDE) model through Netflux software (available at https://github.com/saucermanlab/Netflux) and simulating experimental conditions determined by context-specific data. We constructed a four-stage model revision process to develop a context-specific network model by implementing the automated validation module as a core (available at https://github.com/mkm1712/Automated_Validation). The default parameters of the model are estimated in the first stage from available qualitative data in all contexts. Next, the “Classified Qualitative Validation” module represents all validation data graphically in different contexts and data classifications, to demonstrate the heterogeneous performance of the model predictions. Performing a two-step Morris-Sobol global sensitivity analysis with the context-specific data in the third stage identifies the key signaling reactions, single and pairs, that affect model validity. This provides insight into context-dependency of the signaling network and reduces the model’s degrees of freedom to a practical range for subsequent calibration with semi-quantitative data. Furthermore, a single reaction deletion technique is presented as a parallel approach, with a lesser accuracy but a higher computational efficacy, to identifying key reactions in each context. Finally, all interactions between signaling pathways are explored via adding new reactions to signaling network nodes through the “OR” gates, and to network reactions through the “AND” gates. This analysis uncovers missing context-specific signaling interactions and crosstalks.

thumbnail
Download:
Fig 1. Pipeline of CLASSED approach.

(A) Four-stage revision process to derive a hybrid context-specific model with an automated validation module as a core, and context-specific data and signaling network interaction graph as inputs. (B) Schematic of context-specific global sensitivity analysis (GSA). Two-step Morris-Sobol global sensitivity analysis identifies key reactions of model validity in each context and reduces the model dimension for calibration of main parameters with semi-quantitative data.

https://doi.org/10.1371/journal.pcbi.1008490.g001

Model construction

Prior-knowledge network.

In this study, we made use of a hypertrophy network model (S1 Table) constructed by Ryall et al. [4] as the prior-knowledge network for our approach. The network topology was based on the 14 most established hypertrophic pathways from literature which links receptor inputs to transcription factors activities, cardiomyocyte size, and gene expression. The network contains 106 nodes, 191 reactions and 17 receptor inputs (S1 Fig).

Modified normalized-Hill logic-based differential equation approach.

We applied a modified version of the normalized-Hill LDE approach to convert the interaction graph to a system of ordinary differential equations. Although these modifications are minor, we found that they are necessary for the method to be accurate in various cases. In this regard, we revised the formula of inhibition reaction to address two limitations of the previous version. As shown in S2A and S2B Fig, we expect a decrease in the weight factor (w) of an inhibition reaction increases the activity of its downstream effectors. This is not valid for the previous formulation. Moreover, setting the weight factor to zero in the “single reaction deletion” technique leads to incorrect baseline activities for a few signaling components such as GSK3β and foxo. Both issues have been addressed by updating the inhibition formula. It is noteworthy that in the previous studies, the weight factor was fixed to one in all cases and, consequently, their results were not affected by the mentioned limitations.

Furthermore, we modified the “AND” gate formula. As illustrated in S2C and S2D Fig, we divided the original “AND” gate output to the mean of its inputs to the power of n-1. Where n is the number of “AND” gate inputs. This adjustment increases the “AND” gate output’s sensitivity to the inputs’ perturbations when we have a large number of inputs (>3). It equalizes the order of magnitude between the output and inputs. S2C and S2D Fig illustrate the variations of an “AND” gate output in terms of its two inputs for the original and the modified formula. S2E Fig illustrates the sensitivity of model validation percent (for all data) to the in silico threshold for defining an increase or decrease in the model. While we have usually seen more than 10% relative changes in the activity of stimulated or inhibited species in experiments, capturing these variations by the model requires a lower threshold due to the reduced sensitivity of the model because of utilizing the same default parameters for all reactions and the existence of “AND” and “OR” gates with a large number of inputs in the network. As shown in S2E Fig, the reduced sensitivity of the model can be easily compensated by choosing an in silico threshold equal to 1% without a considerable effect on computational cost (6.6% increase in computational time in comparison with 10% threshold).

Model parameters.

Eight default parameters regulate signaling nodes’ activity in the model. The first three are “EC50”, “n” (Hill-equation parameters) and “WR” (reaction weight) which are specific to intracellular signaling reactions. While “Wi” adjusts the initial baseline activity of all inputs in a resting state, “We” represents the elevated activity of a specifically stimulated input reaction, such as adding an agonist like isoproterenol. “Ymax” determines the maximum activity of each node and is set to zero for the node inhibition. Finally, “Y0” and “tau” control the initial value and time constant of each node activity, respectively. The overexpression of a signaling node is simulated by setting its “Y0” and “tau” equal to “Ymax” and 109 (a large value), respectively. The single reaction deletion was performed by setting the weight factor “WR” of all signaling reactions to zero, one at a time.

Parameter estimation.

We utilized the GlobalSearch optimization method of MATLAB Global Optimization Toolbox to find the optimized default parameters. We selected the GlobalSearch method due to its better performance in finding global minimum in terms of accuracy and repeatability in comparison with other methods such as pattern search and multistart. The GlobalSearch method [43] employs a scatter-search algorithm to create a set of trial points and heuristically chooses the time to execute a local optimization. The “fmincon” solver with the “sqp” algorithm has been selected for the GlobalSearch method. The repeatability of finding global minimum has been examined by repeating the estimation process with 100 random initial sets.

Context-specific data collection

Context-specific data for the cardiac hypertrophy were manually curated from the literature. We obtained qualitative and semi-quantitative data at steady-state or an appropriate time point from research articles with similar cell lines, assay, and experimental conditions with a preference for data from neonatal rat ventricular myocytes. In total, 450 qualitative (all contexts) and 100 semi-quantitative (ISO-specific context) experimental data points have been obtained from 230 articles. The qualitative and semi-quantitative data are displayed in S2 and S3 Tables, respectively.

Global sensitivity analysis

Two global sensitivity analysis methods were conducted sequentially to identify signaling reactions with the most influence on the model agreement with data in each context. In the first step, we utilized Morris’s elementary effects method which is considered a proper screening method for models with numerous parameters and/or when the computational cost for model simulation is relatively high. The “Sampling for Uniformity” strategy with oversampling size, input factor level, and trajectories’ number equal to 300, 8, and 16, respectively, has been applied to generate the Morris sampling data by EE sensitivity package developed by Khare et al. [44].

For the second step, the Sobol sensitivity analysis was employed. This variance-based method is considered a reference for evaluating other global sensitivity analysis approaches [45]. However, the Sobol method is not practical for models with many parameters owing to its extremely high computational cost. Hence, applying a screening approach to reduce the number of estimated parameters prior to Sobol analysis is indispensable for large-scale models. The Sobol-Jansen formula was utilized for estimating the Sobol indices [46]. Also, the 95% confidence intervals for Sobol indices comprising main, total and pair indices were computed by using the Bootstrap method with a resampling size equal to Sobol’s sample size and 100,000 replicates. Considering that no definite value is recommended for the Sobol’ sample size in literature, we performed the sensitivity analysis with different sample sizes from 1,000 to 100,000. The Sobol indices converged to their final value for sample size equal to and greater than 80,000 (4–5 million simulations) with all indexes greater than zero. Each run with sample size equal to 80,000 took eight hours on 20 parallel CPUs with 40 core for each CPU.

Experimental approaches

Materials.

The Janus Kinase 2 protein (JAK2) inhibitor AG 490, Mitogen-Activated Protein Kinase Kinase (MEK) inhibitor U0126, small molecule Gβγ inhibitor Gallein, Ca2+/calmodulin-dependent protein kinase II (CaMKII) inhibitor KN93, Calcineurin (CaN) inhibitor Cyclosporin A (CsA) and 4′,6-diamidino-2-phenylindole (DAPI) were all purchased from Sigma (USA). Isoproterenol hydrochloride (ISO) and Dimethyl sulfoxide (DMSO) were purchased from R&D Systems (USA). A monoclonal anti-sarcomeric actin antibody and Antibodies against p44/42 MAPK (ERK1/2), phospho-p44/42 MAPK (Thr202/Tyr204), α-Tubulin (DM1A) and secondary antibodies Alexa4 568 and 647 were purchased from Cell Signaling Technology (USA). ISO was dissolved in H2O and used at 10 μmol/L final concentration. Additional reagents were diluted in DMSO and added to the serum-free media with the following final concentrations: AG 490 (10 μmol/L), U0126 (1 μmol/L), Gallein (10 μmol/L), KN93 (1 μmol/L), and CsA (0.5 μmol/L). The final concentration of DMSO was lower than 0.3%.

Cell culture.

Cardiac myocytes were isolated from 1–2-day-old Sprague-Dawley rats using Neomyts isolation kit (Cellutron, USA). The cells were cultured in plating media (66% low-glucose Dulbecco's modified eagle media (DMEM), 17% M199, 10% horse serum, 5% fetal bovine serum (FBS), 1% L-Glutamine, 100 U/mL penicillin, and 50 mg/mL streptomycin) at a density of 40,000 cells per well of a 96-well Corning CellBIND plate and 750,000 cell per well of a 6-well Corning CellBIND plate coated with SureCoat (Cellutron, USA). All medium solutions were buffered with HEPES (pH 7.5) to a final concentration of 20 mM. After 48 hours of incubation with daily media change, the medium was replaced with 0.1% FBS medium and the cardiac myocytes were serum-starved for 18 hours before treatment.

Western blot.

Treated or control cardiac myocytes were washed with ice-cold PBS and digested with RIPA buffer [47]. After transferring the cell lysates to tubes, we pass lysates through the needles about 10 times and incubated it on ice for 20 minutes. The samples were centrifuged at 12,000 rpm for 20 minutes at 4°C, and the supernatant liquids were collected. Total protein concentration was quantified using the Pierce BCA protein assay kit (Fisher Scientific, USA). For each sample, 20 μg of total protein were electrophoresed on a 7.5% sodium dodecyl sulfate-polyacrylamide gel (Bio-Rad, USA) and then transferred to Polyvinylidene fluoride (PVDF) membranes (Sigma, USA). After blocking with 0.5X Odyssey (LICOR) blocking solution in TBST (Tris-buffered saline, 0.1% Tween 20) at room temperature for 1 hour, the membranes were incubated with primary antibodies (1:1000) overnight at 4°C. After extensive washing with TBST, proteins were detected with secondary antibodies (1:15000) IrDye 800CW and IrDye 680LT (Licor, USA) followed by quantification using a LICOR Odyssey system. The amounts of phospho-p44/42 MAPK and p44/42 MAPK were normalized to α-tubulin for all samples, and the significance of the difference between the treated and control samples was determined by two-tailed Student’s t-test.

Immunofluorescence and Imaging.

Cell fixation after treatment was conducted with 4% paraformaldehyde (PFA) and followed by permeabilization with 0.2% Triton-X/TBS for 20 min and blocking by 1% BSA/TBS for 1 hour. Then, we incubated the cells with primary antibodies overnight at 4°C and subsequently blocked them by 2% normal goat serum (NGS) for 1 hour before incubating with appropriate secondary antibodies. The cells were stained with DAPI for visualizing nuclei (blue), monoclonal anti-α-Actinin antibody (red) for cardiac myocytes, and phospho-p44/42 MAPK antibody (green) for ERK1/2 activity. The fluorescence images were captured using the Operetta High Content Screening System (Perkin Elmer, USA). The central images were analyzed via MATLAB.

Results

Estimation of default parameters for myocyte-specific context

Qualitative observations representing variations in signaling network nodes in response to hypertrophic agonists were obtained from the literature at the steady-state, or suitable time points for nodes with transient activity. While our logic-based differential equation approach predicts signaling dynamics, due to lack of sufficient kinetic experimental data we focused mostly on steady-state conditions to study causal interactions in cardiac hypertrophy signaling network. All observations were recorded as an increase, decrease, or no change according to their published statistical significance in comparison to the control group in each study. Data were categorized in four classes: Input-Output, Input-Intermediate, Intermediate-Inhibition, and Intermediate-Overexpression (Fig 2A). Validations in the original manuscript reporting this model were from only Input-Output and Input-Intermediate data types [4]. But given that many context-specific responses originate from mediators of signal transduction, we found it important to expand the validation dataset in terms of data classes and numbers (Fig 2B).

thumbnail
Download:
Fig 2. Estimating default parameters enhances model prediction accuracy in an expanded data compendium covering four validation categories.

(A) Classification of myocyte-specific data in four classes based on the function of measured and perturbed signaling nodes in the network. (B) Comparison of the original [4] and expanded datasets in terms of frequency of each data class. (C) Data distribution in the original and expanded datasets regarding the context (hypertrophic agonist). (D) Comparison of the model prediction accuracy in all data classes before and after estimation of default parameters. (E, F) Variations of model validation percent in terms of four default parameters of the model (n, WR, Wi, We). In each figure, the non-variant parameters were fixed to optimum values. The white star shows the optimum parameter set.

https://doi.org/10.1371/journal.pcbi.1008490.g002

In this study, we acquired data for 16 different agonists but largely focusing on four established hypertrophic agonists as the studied contexts involving isoproterenol (ISO), phenylephrine (PE), angiotensin II (Ang II) and biomechanical stretch. The number of data for each context is displayed for the original and expanded datasets in Fig 2C. Each of studied contexts includes 75 qualitative experimental data points. Next, by importing the model interaction graph and expanded dataset to the automated validation module, we computed the model validation percent. The validation module simulates each defined experimental condition; categorizes the simulation results based on their percent of changes from control with 1% change as the threshold; and computes the percent of agreement between model predictions and the matched biochemical measurements.

As shown in Fig 2D, the hypertrophy model with original default parameters from Ryall et al. study [4] predicts the data on Input-Output and Input-Intermediate categories quite well, with more than 80% correct prediction in both classes. Its performance was lower for Intermediate-Overexpression data, with ~60% accuracy. In the Intermediate-Inhibition class of data, we observed a low prediction accuracy (24%).

To test whether the poor accuracy could be explained by non-optimal default parameters, we performed a parameter estimation. The parameter’s values in the original model [4] had been set based on the default values from beta-adrenergic signaling model [39]. In this study, we estimated the four default parameters of the model including Hill coefficient “n”, reaction weight “WR”, initial input weight “Wi” and experimental weight “We” based on the hypertrophic experimental data. As shown in Fig 2D, estimating these default parameters significantly enhances the model prediction accuracy particularly in the Intermediate-Inhibition class of data from 24% to 63%. Also, we observed improvements in Intermediate-Overexpression class with more than 75% correct predictions, similar to the total model validation percent. To examine the extent to which uncertainty in parameters affects the model prediction accuracy, the variations of model validation percent in terms of estimated parameters have been illustrated in Fig 2E and 2F. While a gradual variation in model predictive power was seen for the pair “WR”, “n” in the diagonal direction, its variation for the pair “We”, “Wi” is horizontally and stepwise. This result indicates that while the change in the agonist activity “We” has the lowest effect on model validity, elevation of signaling nodes’ basal activity significantly deteriorates the model validity. Furthermore, variations of estimated parameters in the neighborhood of the optimum point with 77% agreement with data have no significant impact on the model validation percent.

Identifying candidates for context-dependent model revision

After calibrating the four default parameters based on qualitative data of all contexts, we employed validation module’s output to illustrate validation data in regard to the contexts and data classes. Fig 3 displays the classified qualitative validation of the model in ISO-specific context. As can be seen, prediction of network outputs activity (blue nodes in S1 Fig) in response to ISO had no disagreement with experimental data within the input-output class (7/7 correct prediction). Also, validations where intermediate nodes were measured as the output had similarly a high agreement with data (27/30 correct prediction). In ISO-specific context, the model had lower efficacy in Intermediate-Inhibition data class (22/38 correct prediction), and CaMKII- and CaN-related data (2/13 correct prediction). Thus, CaMKII and CaN seemed to be the main sources of inefficiency in the model structure, which we later targeted for model revision. Similar results can be seen in other contexts in regard to the model prediction accuracy for different data classes (S3 Fig). However, we observed different sources of inefficiency, in addition to CaMKII and CaN, in other studied contexts such as NFκβ and RhoA in Ang II and stretch contexts, respectively [48].

thumbnail
Download:
Fig 3. Classified qualitative validation reveals incomplete parts of the ISO-specific network.

Comparison of model predictions with 75 experimental data points in three classes of data including input-output (7/7 correct prediction), input-intermediate (27/30 correct prediction), and intermediate-inhibition (22/38 correct prediction) after isoproterenol stimulation. The red, blue, and gray boxes illustrate increase, decrease, and no change, respectively. In the model, changes in the measured node activity greater than +1% or smaller than -1% have been considered as an increase or decrease, respectively. In experimental data, statistically significant changes in comparison with the control have been considered as an increase or decrease based on the direction of changes.

https://doi.org/10.1371/journal.pcbi.1008490.g003

The main reason for considerable disagreement in the Intermediate-Inhibition class is the context-dependency of the measured data. This reflects its impact on model predictions in Intermediate-Inhibition class. This context-specific response could be due to multiple reasons including inaccurate reaction parameters, erroneous combinatory logic, and incomplete model structure. Inaccurate parameters such as reaction weight could be influencing the individual contributions of each pathway to the observed outputs. In the current state of the model, all reactions have equal reaction weights, when in reality some reactions may be exerting more or less of an influence than others [12,49]. The model could become more sensitive to a given inhibition if certain reactions were given higher weight than others, which is applied to the model in the third stage. Incorrect reaction logic could also lead to a greater insensitivity of model to nodes inhibition, as combinations with "AND" gates are generally more sensitive to a node inhibition than "OR" gates. The low predictive power may also be the outcome of an incomplete model structure, where new reactions are required to correctly predict the experimental outcomes. We address this issue in the final stage of CLASSED approach.

Reactions required for context-specific agreement with experimental data

One of the challenges in analyzing large-scale networks is to systematically determine how modification in individual reactions will affect the comprehensive set of validations. In the same vein, it is complicated to find out which reactions are responsible for the model's inconsistency with experimental data. To address this issue, firstly, we employed automated validation module to test the local effects of each reaction on model validation percent by revalidating the model with context-specific data following systematic, one-by-one removal of each 174 intermediate reactions in the network for each context.

Fig 4A illustrates the results of single reaction deletion technique for four studied contexts in the hypertrophy network. Based on the results, three categories of reactions can be recognized in each context. In the “non-sensitive reactions” category, deleting a reaction has zero or a small positive effect on model validity. Interestingly, the majority of hypertrophy network reactions have been categorized as non-sensitive reactions, which can be explained by several intrinsic properties of the hypertrophy network. As shown in S1 Fig, there are many redundant parallel pathways in the hypertrophy network, which lead to the insensitivity of model to the elimination of a single reaction. Applying the same default parameters for all reactions also intensifies this insensitivity. Moreover, there exist several hubs in the hypertrophy network like Ras, PI3K, p38, and ERK1/2 with many inputs or outputs which transduce the signal regardless of deleting a single reaction. Hence, applying more comprehensive approaches like global sensitivity analysis which consider the variations of all reaction parameters at the same time seems to be essential to scrutinize the role of each reaction in the overall response of the network.

thumbnail
Download:
Fig 4. Single-reaction deletion identifies and categorizes context-specific and shared key reactions.

(A) Changes in model validation percent after one-by-one removal of each 174 intermediate reactions in the hypertrophy network for ISO (cross), PE (circle), AngII (triangle), and stretch (square) contexts. (B) The key reactions in ISO-specific and (C) AngII-specific contexts have been visualized on the hypertrophy network map. The thicker and bolder red, or blue arrows illustrate more decrease, or increase in validation percent, respectively, after removing each reaction.

https://doi.org/10.1371/journal.pcbi.1008490.g004

The “direct reactions” category includes the reactions which are specific to each context and have little or no effect on other contexts. This category mostly comprises the direct downstream reactions of the main receptor for that context. On the other hand, the “shared reactions” category involves reactions that are crucial for multiple contexts and could play a switching role in different context-dependent responses. To determine key reactions, we ranked the reactions in terms of change in the validation percent after removing each reaction. Revision of a model should be done with caution when modifying high score reactions because these reactions are required for a large number of correct predictions. Conversely, non-sensitive reactions when removed from the model indicate a non-sensitive area of the network and could be targeted for model reduction. Under the context of ISO stimulation, single reaction deletion affected model accuracy only for 40 out of 174 model reactions.

To examine how key reactions differ across contexts, the results of the single reaction deletion were visualized on the hypertrophy network map for ISO-specific and AngII-specific contexts in Fig 4B and 4C, respectively. As expected, upstream reactions in the classical β-adrenergic signaling pathway [50] from βAR to protein kinase A (PKA) have the highest effect on model validity in ISO-specific context and would be categorized in the “direct reactions”. Similar results could be obtained for the AngII-specific context in which upstream reactions from AngII receptor (AT1R) to Inositol trisphosphate (IP3) have the highest score. In these two contexts, reactions such as PI3K = >PDK1 = >Akt and MEK12 = >ERK1/2 (ERK1/2 pathway in general) have significant impacts on both contexts and identified as “shared reactions”. These reactions might be the potential targets for further experiments to understand distinct myocyte responses after ISO and angiotensin II stimulations [51,52]. The results of single reaction deletion technique for other contexts, PE-specific and Stretch-specific, are depicted in S4A and S4B Fig. Comparing all contexts, S4C Fig, shows that more shared reactions are involved in the cardiac hypertrophy signaling in PE-specific and AngII-specific contexts than ISO-specific and Stretch-specific contexts.

Whereas the single reaction deletion method is easily interpretable, a limitation is that it does not account for nonlinear interactions between reactions in the hypertrophy network. Therefore, we applied a two-step global sensitivity analysis to detect the key reactions for hypertrophy model validity (for details see description in method's section). In this stage, we selected EC50 as a variable of sensitivity analysis due to the following reasons: 1) the weight factor “WR” was estimated before in the first stage to enhance the model accuracy in predicting myocyte-specific data (EC50 was fixed); 2) based on the LDE formulation, EC50 is the most influential parameter of each reaction [39] and more diverse reactions could be obtained by using EC50 than WR (S5 Fig) as a variant parameter of the Morris method which provides identification of transactivation in each context; and 3) unlike “WR”, using EC50 as a variant parameter predicts similar key reactions for each context regardless of choosing different model outputs (blue nodes in S1 Fig) or model validation percent as the objective function which provides robust results.

The results of Morris sensitivity analysis for ISO-specific context are displayed in Fig 5A and 5B. Since the level of significance for μ* is case-dependent, we assumed μ* equals to 1.5 as a threshold to obtain key reactions in ISO-specific context. While this selection could be changed based on the model structure and the number of parameters (e.g. 1.2 for other contexts), we selected this value to obtain a range of 25 to 40 effective reactions considering that these reactions would be utilized as the inputs of Sobol sensitivity analysis. Considering higher threshold results in the reduced number of inputs for Sobol sensitivity analysis and higher chance of missing an important reaction or nonlinearity effect between reactions. On the other hand, lowering the threshold increases the number of effective reactions which results in a much higher computational cost without gaining much information about key reactions in most cases. Therefore, choosing an appropriate value for the level of significance in Morris results is case-dependent. As shown in Fig 5B, 35 out of 174 reactions have a considerable effect (μ* > 1.5) on model validation percent in ISO-specific context. Interestingly, all reactions were above the diagonal line (μ* = σ) in Fig 5A, which means that they have a dominant interaction (nonlinear) effect on model accuracy. In other words, these reactions primarily influence cell signaling through combination with other network reactions. Thus, when regulating these reactions experimentally, for example to find a treatment, we should consider controlling their pair reactions, too. As the orange dots illustrate the reactions with monotonic effect, most of the key reactions in ISO-specific network have a non-monotonic effect on model agreement with data. This means that we could observe different responses in regulating a specific reaction in β-adrenergic-induced hypertrophy and there is an optimal range for the concentration of the drug. In contrast to the findings of single reaction deletion technique, Morris analysis indicated that the β-adrenergic non-classical pathways [50] play the main role in regulating the overall model response in ISO-specific context. This contradiction seems to be due to the importance of nonlinear interactions in β-adrenergic hypertrophy network. Morris sensitivity analysis’ results for other studied contexts (AngII, PE, stretch) are depicted in S6 Fig. As can be seen, there are more reactions under the diagonal line in other contexts which indicates reactions with higher direct effect on cell overall response. These reactions could be the potential targets for regulating cell response with less concern of interference from other signaling reactions.

thumbnail
Download:
Fig 5. Two-step Global sensitivity analysis identifies key single and pair reactions in an ISO-specific context.

(A) Morris sensitivity analysis result including Morris index (μ*) and standard deviation (σ) illustrates important and non-important reactions in the ISO-specific context. Greater Morris index (μ*) demonstrates more influence on model validation percent. Larger σ to μ* ratio (above μ* = σ diagonal line) demonstrates a more nonlinear effect on model prediction accuracy which implies on more interactions with other network reactions. The orange dots and blue stars illustrate monotonic and non-monotonic effects of each reaction on model validation percent, respectively. The vertical red line indicates the significance level for identifying non-important reactions. (B) Visualization of important reactions in ISO-specific context based on their ranks from Morris sensitivity analysis. The thicker red arrows illustrate the larger Morris index (μ*). (C) Sobol sensitivity analysis reveals importance of nonlinear interactions between network signaling reactions in ISO-specific context. The Total (checked red bars) and Main (hashed blue bars) Sobol indexes and pair synergic effects are illustrated for important reactions (35/174) in ISO-specific context obtained from Morris sensitivity analysis. The error bars represent 95% confidence interval of a mean. More difference between the Total and Main indexes demonstrates a more nonlinear effect. The pair synergic effect displays the ratio of the second-order effect of each reaction pair to the sum of higher-order effects of its reactions.

https://doi.org/10.1371/journal.pcbi.1008490.g005

In the second step, we fixed the varying parameter (EC50) of all non-important reactions identified by Morris sensitivity analysis and performed Sobol sensitivity analysis on the reduced version of the model. By employing the Sobol-Jansen formulation, we obtained Monte Carlo estimation of total, main (first-order), and second-order Sobol indices for each reaction. As shown in Fig 5C, ISO = >βAR and Gsα = >Gβγ reactions have the highest total (0.44) and main (0.17) effects on model validation percent in ISO-specific context, respectively. Regarding the synergy effect, pairs of (ISO = >βAR, βAR = > Gsα), (Gβγ = >PI3K, Akt = >GSK3β), (gp130LIFR = >JAK, JAK = >Ras) and (PLCB = >IP3, IP3 = > Calcium) have the largest relative synergy in the ISO-specific contexts with more than 30% synergy effect.

These results are consistent with Morris method outcomes and confirm the necessity of utilizing global sensitivity analysis to predict context-specific key reactions in large-scale signaling networks. The key reactions of model agreement with data in each context are shown in S7 Fig. Finally, after identifying the key single and pair reactions in ISO-specific context, we developed an ISO-specific model of cardiac hypertrophy by calibrating the EC50 of 12 main reactions making use of both qualitative data and semi-qualitative data (S3 Table) in ISO-specific context. The S8 Fig illustrates enhancement of model prediction accuracy (from 47% to 60%) in ISO-specific context.

Predicting new crosstalks in ISO-specific context

In general, developing a context-specific model from the prior-knowledge model needs revision of both parameters and structure of the model. In the previous section, we revised the prior model key parameters to obtain an ISO-specific model for cardiac hypertrophy. But there exist some experimental observations in β-adrenergic-induced hypertrophy that the calibrated model still cannot predict (S8 Fig). Therefore, we explored all possible new crosstalks in the network by systematic, one-by-one adding of new interactions between network nodes through the “OR” gate. We also explored adding new interactions from each network node to the existing reactions through the ‘AND” gate. Similar to the results of classified qualitative validation, we identified CaMKII and CaN as two source nodes with the highest potential in increasing the model agreement with data. In the ISO-specific context, the top-ranking reactions which improved model agreement with experimental data were new interactions from CaMKII and CaN to the upstream nodes of the Ras/Raf/MEK/ERK1/2 pathway.

Fig 6 illustrates the top nominations for crosstalks in ISO-specific context by different logic gates associated with CaMKII. As can be seen, individual crosstalks with the “AND” gate lead to higher increase (maximum +15%) in model prediction accuracy compared with crosstalks with the “OR” gate (maximum +10%). The model predicts that combinations of CaMKII with the upstream reactions of the Ras have the most positive effect on the model validation (8–15% increase). This prediction suggests a contribution of CaMKII in the activation of β-adrenergic non-classical pathways. Among the model predictions, reactions associated with JAK and Gβγ have the highest ranks. Among the combination of CaMKII with other network nodes through the “OR” gate, the positive feedback from CaMKII to Calcium has the highest rank. The other predictions indicate the role of CaMKII in the transactivation of other hypertrophic pathways in β-adrenergic-induced hypertrophy.

thumbnail
Download:
Fig 6. Exploring possible interactions nominates CaMKII as the key signaling node for missing crosstalks in ISO-specific network.

(A) Model predictions for crosstalks between CaMKII and network reactions through the “AND” gate. Each arrow represents a reaction that by adding an “AND” gate from CaMK to that reaction, we observed a significant improvement in model accuracy. (B) Model predictions for crosstalks between CaMKII and network nodes through the “OR” gate. The thicker red arrows illustrate reactions with higher positive effect on the validation percent in ISO-specific model.

https://doi.org/10.1371/journal.pcbi.1008490.g006

Exploring all possible “AND” gates in the hypertrophy network (18444 cases) showed that for the ISO-specific network, 12.1% of new interactions decreased model validation percent. 5.4% of cases increased model validation percent, of which only 1.4% resulted in improving prediction of more than one data point of 75 validation data points. As expected, 82.5% of cases did not affect model predictions. For the “OR” gate in ISO-specific context (11236 cases), 39% of cases resulted in a decrease in model validation percent. 2.7% of cases increased model validation percent which only 1.2% resulted in improving prediction of more than one data point. 58.3% of “OR” gate cases did not affect model validation percent.

In addition to the CaMKII, the model predicts a significant increase (8–9%) in model validation by adding new interactions through “AND” gate from CaN to the upstream reactions which activate Ras in the hypertrophy network (S9A Fig). Also, after adding the positive feedback from CaMKII to Calcium which was predicted by the model, we repeated the adding reaction process with the modified model and observed that interactions between CaN and other reactions through the “AND” gate have the highest rank in enhancing model validation percent (5–7%), see S9B Fig.

Validation of new crosstalks in β-adrenergic-induced hypertrophy

In the previous section, the model predicted three main crosstalks in ISO-specific context. First, a positive feedback from CaMKII to calcium through an “OR” gate has been predicted. This feedback signal has been verified by numerous previous studies on cardiac electrophysiology and calcium signaling [53–55]. For instance, Lee et al. [56] demonstrated that CaMKII is involved in the facilitation of L-type calcium channels in cardiac myocytes. Moreover, a review article by Grimm and Brown [57] describes CaMKII targets in β-adrenergic signaling, explains regulation of cardiac calcium signaling by CaMKII through a positive feedback with mediators like Ryanodine receptor, Phospholamban, and L-type Calcium channels [58].

The other two predictions of our model indicate the role of CaMKII and CaN in the activation of non-classical β-adrenergic pathways after ISO stimulation. While numerous studies have examined either the role of CaMKII and CaN in cardiac hypertrophy [13,59] or Ras/Raf/MEK/ERK1/2 pathway in cardiac hypertrophy [60–62], very few have examined cross-talk between calcium and MAPK pathways for cardiac hypertrophy. Zou et al. [63] investigated ISO-induced ERK1/2 phosphorylation in cultured cardiomyocytes, which they found to require CaN but not CaMKII. Sanna et al. [64] founded that calcineurin-NFAT and MEK1-ERK1/2 signaling pathways are interdependent and by forming a complex in cardiac myocytes coregulate the hypertrophic response.

To experimentally validate the model-predicted crosstalks between MAPK signaling and CaMKII/CaN, we assessed the effects of CaMKII, CaN, JAK2, MEK12 and Gβγ inhibition on ISO-induced ERK1/2 phosphorylation in primary neonatal rat cardiomyocytes (Fig 7A). Among nodes upstream of Ras, we chose to target Gβγ due to its predicted central role in ISO-specific context (Figs 5C, 6A and S7A) and JAK2 due to its mediatory position for high rank reactions (Fig 6A). ERK1/2 phosphorylation was assessed by Western blotting with or without isoproterenol for 10 minutes after 30 min pre-treatment with specified inhibitors. As expected, ISO increased ERK1/2 phosphorylation, and pretreatment with MEK12 inhibitor suppressed ISO-induced ERK12 phosphorylation (Fig 7B). Without ISO, only Gβγ inhibition with Gallein, and its combination with CaMKII inhibition (KN93), significantly diminished the basal level of ERK phosphorylation (80% decrease) compared to the control (DMSO). With ISO treatment, CaMKII inhibition by KN93 increased ERK1/2 phosphorylation significantly (140% increase in comparison with ISO), inhibition with Gβγ had no statistically significant effect, and combination of CaMKII and Gβγ inhibition significantly diminished ERK1/2 phosphorylation (70% decrease in comparison with ISO). These data indicate that combined CaMKII and Gβγ activation are required for full ISO-induced ERK1/2 phosphorylation. These experimental observations validate the model prediction of a missing crosstalk which connects CaMKII to Gβγ = >Ras reaction via an “AND” gate (Figs 6A and S9A).

thumbnail
Download:
Fig 7. Inhibition of CaMKII-Gβγ or CaN-Gβγ suppresses ISO-induced ERK1/2 activity.

(A) Effects of pretreatment for 30 min with five inhibitors including JAK-i (AG490), CaMKII-i (KN93), CaN-i (CsA), Gβγ-i (Gallein), and MEK-i (U0126) were examined on basal and ISO-induced ERK1/2 activity. (B) Upper panel shows western blotting of ERK1/2 phosphorylation after pretreatment with inhibitors in comparison with control (DMSO). Lower panel shows ERK1/2 phosphorylation after ISO stimulation with pretreatment with same inhibitors. *represents statistically significant change (p<0.05) in comparison with control (DMSO) for the upper panel, or with ISO in the lower panel. # represents significant change (p<0.05) compared with control (DMSO) in the lower panel. Data were collected from 3 independent experiments (Mean±SEM). (C) Immunofluorescence images show ISO-induced ERK1/2 phosphorylation and its suppression by MEK inhibitor in myocytes (alpha-actinin positive) and fibroblasts (alpha-actinin negative). The scale bar is 50 μm (D) Immunofluorescence imaging of ERK1/2 phosphorylation shows increase of ERK1/2 phosphorylation after CaMKII inhibition with KN93 and its suppression after pretreatment with Gβγ+CaMKII and Gβγ+CaN combined inhibitions. The scale bar is 200 μm.

https://doi.org/10.1371/journal.pcbi.1008490.g007

Calcineurin (CaN) inhibition also modified ISO-induced ERK1/2 phosphorylation, but to a lesser degree. CaN inhibition alone had no significant effect on ISO-induced ERK1/2 phosphorylation, but combined CaN and Gβγ inhibition reduced ERK1/2 activity significantly by 60% in comparison with ISO. Thus, like CaMKII, CaN contributes with Gβγ to enhance ERK1/2 activity.

The western blotting is not usually measuring only cardiomyocytes’ response because protein lysates comprise both cardiomyocytes and non-cardiomyocytes (mostly fibroblasts). Therefore, we repeated the above perturbations but measured ERK1/2 localization and phosphorylation by immunofluorescence, which can discriminate between responses in cardiomyocytes and non-cardiomyocytes (Fig 7C and 7D). In cardiomyocytes, ISO stimulation for 10 min increased ERK1/2 phosphorylation primarily in the perinuclear area. In contrast, fibroblasts exhibited much stronger ERK1/2 phosphorylation, consistent with nuclear translocation (Fig 7C). Pretreatment with MEK12 inhibitor suppressed ISO-induced ERK1/2 activity in both cardiomyocytes and fibroblasts. We further examined the influence of pretreatment with different inhibitors on ISO-induced ERK1/2 activity in both cardiomyocytes and fibroblasts by immunofluorescence (Fig 7D). Response of ERK1/2 phosphorylation in cardiomyocytes was qualitatively similar to the results from western blotting. Both cardiomyocytes and fibroblasts exhibited the strongest ERK1/2 phosphorylation with combined ISO stimulation and CaMKII inhibition by KN93. Like in western blotting, we found significant decrease in ISO-induced ERK1/2 activity in the combination of Gβγ inhibitor with CaMKII or CaN inhibitors. The results further support the finding of crosstalk between calcium and Gβγ signaling in both cardiomyocytes and fibroblasts. To quantify the difference between cell types in ISO-induced ERK1/2 response, we thresholded images based on alpha-actinin expression. As illustrated in S10 Fig, ISO-induced ERK1/2 activation was much stronger in fibroblasts than cardiomyocytes, but both cells had responses qualitatively consistent with the western blotting.

Based on the experimentally validated inferred interactions, we updated the hypertrophy model for the ISO-specific context (S4 Table) as illustrated in Fig 8. The red solid lines describe the CLASSED-inferred reactions that were experimentally validated in the ISO-specific context. While revision on model parameters enhanced model prediction accuracy from 48% (original model) to 72% in ISO-specific context, adding new interactions increases this prediction accuracy to 86.7% (65 of 75 predictions). Moreover, red dashed-lines in Fig 8 suggest additional new interactions that are consistent with the unanticipated effect of CaMKII inhibition on enhancing ISO-induced ERK1/2 phosphorylation. Consistent with the new suggested interactions in dashed lines, some studies with other cell types have shown that KN93 increases the MEK12/ERK1/2 activity through EGFR activation [65], but this activation is not through Ras/Raf1 [66] and directly influences ERK1/2 activity. Also, there are several studies which confirm the transactivation of other pathways like EGFR during β-adrenergic receptor activation and their Ras-independent influence on ERK1/2 activity [67].

thumbnail
Download:
Fig 8. Revised model of β-adrenergic cardiac hypertrophy.

The solid red arrows exhibit predicted and validated interactions in ISO-specific context. The dashed red arrows illustrate proposed interactions predicted by model based on the new experimental data on the effects of CaMKII inhibition on ISO-induced ERK1/2 activity.

https://doi.org/10.1371/journal.pcbi.1008490.g008

Discussion

The present study provides a systematic approach, CLASSED, for developing context-specific models for large-scale signaling networks. CLASSED has been designed to systematically revise a prior-knowledge network in accordance with context-specific perturbation data in qualitative or semi-quantitative formats. Given that the logic-based ODE formulation implemented in this approach could capture the dynamics of signaling networks, our method has the potential to accommodate quantitative data. The automated validation module generates LDE models and outputs model prediction accuracy and validation data for each and overall experimental conditions. In this regard, we employed expanded context-specific hypertrophic data curated from the literature and a large-scale hypertrophy network [4,42] to identify the key signaling reactions for model validity in four hypertrophic contexts. We also developed an ISO-specific hypertrophy model through the four-stage model revision process. The ISO-specific model predicted new crosstalks between the calcium/calmodulin-dependent pathways (CaMKII and CaN) and the upstream reactions of Ras, and a positive feedback from CaMKII to calcium. These new interactions were examined and validated by previous experimental data as well as new experiments with primary cardiomyocytes. Our experiments revealed a synergic effect of CaMKII-Gβγ and CaN-Gβγ pairs in regulating ISO-induced ERK1/2 phosphorylation which is an important player of cardiac hypertrophy [68,69].

Context-specific network model

Context-dependency is a central but poorly understood aspect of cell signaling [70]. Recently, this concept has been emphasized in developing new drug treatments, particularly personalized treatments for cancer patients. In several studies by the Saez-Rodriguez group, the authors showed that patient- or context-specific logic-based ODE models enhanced prediction of cancer cells’ response in comparison with non-context-specific models [40,71,72]. The Alberts group have highlighted the efficacy of discrete dynamic models to understand biological signaling and regulatory networks in different contexts [37]. In this regard, they developed different models to study T cell large granular lymphocyte (T-LGL) survival signaling [73,74], immune response to B. bronchiseptica during its infection of the lower respiratory system [75] and plant guard cell signaling [76] and demonstrated the model efficacy in different biological systems.

Apart from enabling dynamic predictions on large-scale signaling networks, logic-based ODE models have the capability to model context-dependent networks despite limited availability of patient- or context-specific data. Although in some fields data is more plentiful, such as with ex vivo high-throughput screening of cancer biopsies [40] or by widespread data sharing as in the DREAM challenges [77], most fields face challenge in insufficient data where prior-knowledge based models may be helpful. Furthermore, these approaches provide the opportunity to integrate the data from different labs and experimental studies to validate a model.

In our study, the prior-knowledge model reactions are supported with at least two independent studies [4]. For context-specific data, the majority of the experimental data are supported by at least two independent studies. Since diversity of data is very important for model prediction, in parts of the network that limited data are available, using one study that includes data from at least three replicates is acceptable as long as is not in complete disagreement with the other data related to that part of network.

CaMKII-Gβγ and CaN-Gβγ: new crosstalks in the ISO-specific context

Through systematic revision of the hypertrophy network model, we predicted a new synergistic cross-talk between CaMKII and Gβγ in ISO-induced hypertrophy. This prediction motivated new experiments, which found that while ISO-induced ERK1/2 phosphorylation was not suppressed by pharmacologic inhibition of CaMKII or Gβγ alone, combined inhibition significantly suppresses the ISO-induced ERK1/2 phosphorylation in both cardiac myocytes and fibroblasts. We also experimentally validated similar crosstalk between CaN and Gβγ, which had a more subtle quantitative influence on ERK1/2 phosphorylation. Moreover, while Zou et al. [63] reported no effect of CaMKII on ERK1/2 activity in ISO-specific context, our observations from western blot and immunofluorescence experiments showed an increase in the ERK1/2 activity after CaMKII inhibition with KN93 in both cardiomyocytes and fibroblasts.

The context-specificity of these interactions is emphasized by comparison of our results and those of Cipolletta et al. [62] on CaMKII/ERK1/2 interaction in PE-treated H9C2 cardiomyoblasts. While they found that CaMKII inhibition reduces PE-induced ERK1/2 activity and its nuclear accumulation, our data with ISO-induced ERK1/2 phosphorylation show enhancement by CaMKII inhibition alone or suppression with combined inhibition of CaMKII and Gβγ. In terms of cell type, while a number of studies found similarities in response for neonatal rat cardiomyocytes and H9C2 cardiomyoblasts [78], others have reported distinct responses for these two cell types [79]. In terms of stimuli, the model predicts more contribution of calcium to Ras/Raf1/MEK/ERK1/2 pathway through PKC in PE-specific context than ISO-specific context (S7 Fig). The effect of CaMKII inhibition on calcium could result in the suppression of ERK1/2 activity in PE-specific context. In another study by Illario et al. [80], the authors reported decrease in ERK1/2 phosphorylation after CaMKII inhibition for two non-myocyte cell lines after integrin stimulation by fibronectin. These studies indicate the context-specificity of ERK1/2 response.

In contrast to the few studies on crosstalks between CaMKII and ERK1/2, many studies have explored crosstalks between CaN and ERK1/2. For instance, expression of activated calcineurin transgene in mice showed ERK1/2 activation in the heart [81]. The results of Zou et al. study [63] also showed activation of ERK1/2 signaling through a calcineurin-dependent mechanism in ISO-specific context. On the other hand, several studies reported ERK1/2 signaling can regulate calcineurin–NFAT activity through stimulation or inhibition based on the context and present NFAT factors [60]. In our study, the ISO-specific model predicts the potential effect of CaN on ERK1/2 activity, and the validation experiment indicates the additional requirement of crosstalk with Gβγ.

Furthermore, our experiments showed a significant difference in ISO-induced ERK1/2 phosphorylation between cardiomyocytes and fibroblasts in terms of magnitude and localization (S10 Fig). While in this case both cardiomyocytes and fibroblasts had similar qualitative responses to inhibitors, improving cell purity or utilizing alternate methods like immunofluorescence may be critical in future studies of cell-specific responses.

Limitations and future directions

While in this study we initially used default parameters for major of signaling reactions, tuning of more reaction parameters can improve the model prediction accuracy as more context-specific data becomes available. The availability of high-quality databases for context-specific data like the current databases in cancer systems biology [82–84] could result in more accurate and robust model predictions. As richness and diversity of data affect model accuracy and its ability to confer new crosstalks, we tried to obtain diverse data from different pathways and various contexts as much as we could. Moreover, adding only a single reaction at a time in the final stage of CLASSED approach limits its capacity to infer new pair interactions with synergic effect in the signaling network.

Although applying the CLASSED approach to the hypertrophy signaling network in the β-adrenergic context resulted in predictions of new crosstalks that were validated by experiments, it is possible that in some cases, experiments do not support model predictions. This inconsistency can be a result of the following limitations: 1- The prediction power of the current approach is dependent on the accuracy of the prior knowledge network and, every false reaction could affect model predictions. 2- Diversity and accuracy of context-specific experimental data considerably affect model predictions. In parts of the network that we do not have enough context-specific data, model predictions should be considered with more caution. 3- Since this approach nominates considerable numbers of reactions as missing interactions, selecting the final cases for validation experiments is not always straightforward and requires comprehensive knowledge of the studied subject.

In this study, we employed our approach to develop an ISO-specific model of cardiac hypertrophy. For the future, curation of more context-specific and semi-quantitative data could enlighten the context-dependent characteristics of cardiac myocytes. Furthermore, validating the predictions of context-specific models with in vivo data could illustrate the contribution of other factors in the specificity of cell response in each context.

Supporting information

S1 Table. Prior-knowledge network model of cardiac hypertrophy.

https://doi.org/10.1371/journal.pcbi.1008490.s001

(XLSX)

S2 Table. Expanded context-specific qualitative data curated from literature.

https://doi.org/10.1371/journal.pcbi.1008490.s002

(XLSX)

S3 Table. Semi-quantitative ISO-specific data.

https://doi.org/10.1371/journal.pcbi.1008490.s003

(XLSX)

S4 Table. Revised network model of β-adrenergic cardiac hypertrophy.

https://doi.org/10.1371/journal.pcbi.1008490.s004

(XLSX)

S1 Fig. Prior-knowledge cardiac hypertrophy signaling network.

The network comprises 106 nodes, 191 reactions and 17 receptor inputs (green color) [4]. Seven nodes have been defined as model outputs (blue color).

https://doi.org/10.1371/journal.pcbi.1008490.s005

(TIF)

S2 Fig. Modifications of normalized-Hill LDE approach.

(A, B) Comparison of the original [39] and modified versions of the inhibition and (C, D) “AND” gate formula. (E) The sensitivity of model validation percent to in silico threshold for determination of a change in model outputs.

https://doi.org/10.1371/journal.pcbi.1008490.s006

(TIF)

S3 Fig. Classified qualitative validation of the model.

All qualitative data in four classes (A) input-output (B) input-intermediate (C) intermediate overexpression, and (D) intermediate inhibition are compared with model predictions except for ISO-specific context. The red, blue, and gray boxes illustrate increase, decrease, and no change, respectively. In the model, variations of the measured node activity greater than +1% or smaller than -1% have been considered as an increase or decrease, respectively. Statistically significant changes in comparison with control have been considered for variations in experimental data.

https://doi.org/10.1371/journal.pcbi.1008490.s007

(PDF)

S4 Fig. Key reactions identified by a single-reaction deletion technique.

(A) PE-specific and (B) Stretch-specific contexts have been visualized on the hypertrophy network map. The thicker and bolder red or blue arrows illustrate more decrease or increase in validation percent after removing each reaction, respectively. (C) Categorizing hypertrophy network reactions in three categories of “non-sensitive”, “direct”, and “shared” reactions in ISO, AngII, Stretch and PE-specific contexts.

https://doi.org/10.1371/journal.pcbi.1008490.s008

(TIF)

S5 Fig. Changing variant parameter from WR to EC50 and its effects on Morris sensitivity analysis.

The left diagram illustrates the Morris index (μ*) and its standard deviation (σ) for each reaction in ISO-specific network by variating the WR parameter. Consequently, important reactions of ISO-specific network are illustrated in the right diagram. Greater Morris index (μ*) illustrates more influence on model validation percent. Larger σ to μ* ratio for each reaction (above μ* = σ diagonal line) demonstrates a more nonlinear effect on model prediction accuracy (interaction with other reactions). The orange dots and blue stars illustrate the reactions with monotonic and non-monotonic effects in ISO-specific network, respectively. The vertical red line shows the significance level for identifying non-important reactions. In the network view, the thicker red arrows illustrate reactions with larger Morris index (μ*).

https://doi.org/10.1371/journal.pcbi.1008490.s009

(TIF)

S6 Fig. Morris sensitivity analysis identifies non-important reactions in each context.

(A) The Morris index (μ*) and standard deviation (σ) for each reaction (left) and the main reactions for model validity (right) are illustrated for AngII-specific, (B) PE-specific and (C) Stretch-specific contexts. Greater Morris index (μ*) illustrates more influence on model validation percent. Larger σ to μ* ratio for each reaction (above μ* = σ diagonal line) demonstrates a more nonlinear effect on model prediction accuracy (interaction with other reactions). The orange dots and blue stars illustrate the reactions with monotonic and non-monotonic effects in ISO-specific network, respectively. The vertical red line shows the significance level for identifying non- important reactions.

https://doi.org/10.1371/journal.pcbi.1008490.s010

(TIF)

S7 Fig. Key reactions identified by Sobol sensitivity analysis in four hypertrophic contexts.

(A) The main reactions regulating model validity in ISO-, (B) PE-, (C) AngII—and (D) Stretch-specific contexts are illustrated on the hypertrophy network map. Thicker red arrows illustrate the larger Total Sobol sensitivity index.

https://doi.org/10.1371/journal.pcbi.1008490.s011

(TIF)

S8 Fig. Validation of ISO-specific network model with qualitative and semi-qualitative data.

All data in ISO-specific context (75 qualitative and 100 semi-qualitative) are compared with model predictions before and after tuning of 12 reaction parameters (EC50). Semi-qualitative data has been categorized in 3 levels of increase (LH: Low High (1.01–2 folds), MH: Medium High (2–5 folds), HH: High High (higher than 5 folds)) and decrease (LL: Low Low (0.5–0.99 folds), ML: Medium Low (0.2–0.5 folds), HL: High Low (less than 0.2 folds)). See S3 Table for detailed information.

https://doi.org/10.1371/journal.pcbi.1008490.s012

(PDF)

S9 Fig. Nominations of Calcineurin crosstalks in ISO-specific network.

(A) Crosstalks identified by model between CaN and hypertrophy network reactions with the “AND” gate before and (B) after adding positive feedback predicted by model from CaMKII to Calcium. The thicker red arrows illustrate reactions with higher positive effect on ISO-specific model validation percent.

https://doi.org/10.1371/journal.pcbi.1008490.s013

(TIF)

S10 Fig. Discrepancy in ISO-induced ERK1/2 response for cardiac myocytes and fibroblasts.

(A) The immunofluorescence results display different ERK1/2 phosphorylation in ISO-specific context in terms of translocation and (B) magnitude between cardiac myocytes and fibroblasts. The single and double hashtags (#) display statistically significant changes (p<0.05) after ISO stimulation for cardiac fibroblasts and myocytes, respectively, in comparison with control (DMSO). The single and double stars (*) exhibit statistically significant changes (p<0.05) for cardiac fibroblasts and myocytes, respectively, in comparison with ISO. Data were collected from 3 independent experiments (Mean±SEM)

https://doi.org/10.1371/journal.pcbi.1008490.s014

(TIF)

Acknowledgments

The authors thank Dr. Chun and Dr. Woo for their helpful comments.

References

  1. 1. Schaub MC, Hefti MA, Harder BA, Eppenberger HM. Various hypertrophic stimuli induce distinct phenotypes in cardiomyocytes. J Mol Med (Berl). 1997;75(11–12):901–20. pmid:9428623
  2. 2. Nakamura M, Sadoshima J. Mechanisms of physiological and pathological cardiac hypertrophy. Nat Rev Cardiol. 2018 Jul;15(7):387–407. pmid:29674714
  3. 3. Ryall KA, Bezzerides VJ, Rosenzweig A, Saucerman JJ. Phenotypic screen quantifying differential regulation of cardiac myocyte hypertrophy identifies CITED4 regulation of myocyte elongation. J Mol Cell Cardiol. 2014 Jul;72:74–84. pmid:24613264
  4. 4. Ryall KA, Holland DO, Delaney KA, Kraeutler MJ, Parker AJ, Saucerman JJ. Network reconstruction and systems analysis of cardiac myocyte hypertrophy signaling. J Biol Chem. 2012 Dec;287(50):42259–68. pmid:23091058
  5. 5. Kang JH, Lee H-S, Kang Y-W, Cho K-H. Systems biological approaches to the cardiac signaling network. Brief Bioinform. 2016 May;17(3):419–28. pmid:26108228
  6. 6. Chung E, Heimiller J, Leinwand LA. Distinct cardiac transcriptional profiles defining pregnancy and exercise. PLoS One. 2012;7(7):e42297. pmid:22860109
  7. 7. Maillet M, van Berlo JH, Molkentin JD. Molecular basis of physiological heart growth: fundamental concepts and new players. Nat Rev Mol Cell Biol. 2013 Jan;14(1):38–48. pmid:23258295
  8. 8. Tamai J, Hori M, Kagiya T, Iwakura K, Iwai K, Kitabatake A, et al. Role of alpha 1-adrenoceptor activity in progression of cardiac hypertrophy in guinea pig hearts with pressure overload. Cardiovasc Res. 1989 Apr;23(4):315–22. pmid:2574075
  9. 9. Zhao M, Fajardo G, Urashima T, Spin JM, Poorfarahani S, Rajagopalan V, et al. Cardiac pressure overload hypertrophy is differentially regulated by β-adrenergic receptor subtypes. Am J Physiol Heart Circ Physiol. 2011 Oct;301(4):H1461–70. pmid:21705675
  10. 10. Dorn GW 2nd. Adrenergic pathways and left ventricular remodeling. J Card Fail. 2002 Dec;8(6 Suppl):S370–3. pmid:12555147
  11. 11. Zornoff LAM, Paiva SAR, Duarte DR, Spadaro J. Ventricular remodeling after myocardial infarction: concepts and clinical implications. Arq Bras Cardiol. 2009 Feb;92(2):150–64. pmid:19360249
  12. 12. Heineke J, Molkentin JD. Regulation of cardiac hypertrophy by intracellular signalling pathways. Nat Rev Mol Cell Biol. 2006 Aug;7(8):589–600. pmid:16936699
  13. 13. Molkentin JD, Lu JR, Antos CL, Markham B, Richardson J, Robbins J, et al. A calcineurin-dependent transcriptional pathway for cardiac hypertrophy. Cell. 1998 Apr;93(2):215–28. pmid:9568714
  14. 14. Baker KM, Chernin MI, Schreiber T, Sanghi S, Haiderzaidi S, Booz GW, et al. Evidence of a novel intracrine mechanism in angiotensin II-induced cardiac hypertrophy. Regul Pept. 2004 Aug;120(1–3):5–13. pmid:15177915
  15. 15. Morisco C, Zebrowski DC, Vatner DE, Vatner SF, Sadoshima J. Beta-adrenergic cardiac hypertrophy is mediated primarily by the beta(1)-subtype in the rat heart. J Mol Cell Cardiol. 2001 Mar;33(3):561–73. pmid:11181023
  16. 16. Shin S-Y, Kim T, Lee H-S, Kang JH, Lee JY, Cho K-H, et al. The switching role of β-adrenergic receptor signalling in cell survival or death decision of cardiomyocytes. Nat Commun. 2014 Dec;5:5777. pmid:25517116
  17. 17. Khalilimeybodi A, Daneshmehr A, Sharif-Kashani B. Investigating β-adrenergic-induced cardiac hypertrophy through computational approach: classical and non-classical pathways. J Physiol Sci [Internet]. 2018 Jul 3;68(4):503–20. Available from: http://link.springer.com/10.1007/s12576-017-0557-5 pmid:28674776
  18. 18. Yang JH, Saucerman JJ. Computational models reduce complexity and accelerate insight into cardiac signaling networks. Circ Res. 2011 Jan;108(1):85–97. pmid:21212391
  19. 19. Shin S-Y, Yang HW, Kim J-R, Heo W Do, Cho K-H. A hidden incoherent switch regulates RCAN1 in the calcineurin-NFAT signaling network. J Cell Sci. 2011 Jan;124(Pt 1):82–90. pmid:21172821
  20. 20. Khalilimeybodi A, Daneshmehr A, Kashani BS. Ca2+-dependent calcineurin/NFAT signaling in β-adrenergic-induced cardiac hypertrophy. Gen Physiol Biophys [Internet]. 2018;37(01):41–56. Available from: http://www.elis.sk/index.php?page=shop.product_details&flypage=flypage.tpl&product_id=5545&category_id=143&option=com_virtuemart pmid:29424351
  21. 21. Cooling MT, Hunter P, Crampin EJ. Sensitivity of NFAT cycling to cytosolic calcium concentration: implications for hypertrophic signals in cardiac myocytes. Biophys J. 2009 Mar;96(6):2095–104. pmid:19289036
  22. 22. Abou-Jaoudé W, Traynard P, Monteiro PT, Saez-Rodriguez J, Helikar T, Thieffry D, et al. Logical Modeling and Dynamical Analysis of Cellular Networks. Front Genet. 2016;7:94. pmid:27303434
  23. 23. Bass GT, Ryall KA, Katikapalli A, Taylor BE, Dang ST, Acton ST, et al. Automated image analysis identifies signaling pathways regulating distinct signatures of cardiac myocyte hypertrophy. J Mol Cell Cardiol. 2012 May;52(5):923–30. pmid:22142594
  24. 24. Chasman D, Fotuhi Siahpirani A, Roy S. Network-based approaches for analysis of complex biological systems. Curr Opin Biotechnol. 2016 Jun;39:157–66. pmid:27115495
  25. 25. Hill SM, Heiser LM, Cokelaer T, Unger M, Nesser NK, Carlin DE, et al. Inferring causal molecular networks: empirical assessment through a community-based effort. Nat Methods. 2016 Apr;13(4):310–8. pmid:26901648
  26. 26. Mattiazzi Usaj M, Styles EB, Verster AJ, Friesen H, Boone C, Andrews BJ. High-Content Screening for Quantitative Cell Biology. Trends Cell Biol. 2016 Aug;26(8):598–611. pmid:27118708
  27. 27. Rozum JC, Albert R. Identifying (un)controllable dynamical behavior in complex networks. PLoS Comput Biol. 2018 Dec;14(12):e1006630. pmid:30532150
  28. 28. Janes KA, Yaffe MB. Data-driven modelling of signal-transduction networks. Nat Rev Mol Cell Biol. 2006 Nov;7(11):820–8. pmid:17057752
  29. 29. Bhalla US. Understanding complex signaling networks through models and metaphors. Prog Biophys Mol Biol. 2003 Jan;81(1):45–65. pmid:12475569
  30. 30. Tian Y, Zhang B, Hoffman EP, Clarke R, Zhang Z, Shih I-M, et al. Knowledge-fused differential dependency network models for detecting significant rewiring in biological networks. BMC Syst Biol. 2014 Jul;8:87. pmid:25055984
  31. 31. Shu L, Chan KHK, Zhang G, Huan T, Kurt Z, Zhao Y, et al. Shared genetic regulatory networks for cardiovascular disease and type 2 diabetes in multiple populations of diverse ethnicities in the United States. PLoS Genet. 2017 Sep;13(9):e1007040. pmid:28957322
  32. 32. Levine JH, Simonds EF, Bendall SC, Davis KL, Amir ED, Tadmor MD, et al. Data-Driven Phenotypic Dissection of AML Reveals Progenitor-like Cells that Correlate with Prognosis. Cell. 2015 Jul;162(1):184–97. pmid:26095251
  33. 33. Molinelli EJ, Korkut A, Wang W, Miller ML, Gauthier NP, Jing X, et al. Perturbation biology: inferring signaling networks in cellular systems. PLoS Comput Biol. 2013;9(12):e1003290. pmid:24367245
  34. 34. Tan PM, Buchholz KS, Omens JH, McCulloch AD, Saucerman JJ. Predictive model identifies key network regulators of cardiomyocyte mechano-signaling. PLoS Comput Biol. 2017 Nov;13(11):e1005854. pmid:29131824
  35. 35. Samaga R, Klamt S. Modeling approaches for qualitative and semi-quantitative analysis of cellular signaling networks. Cell Commun Signal. 2013 Jun;11(1):43. pmid:23803171
  36. 36. Morris MK, Saez-Rodriguez J, Sorger PK, Lauffenburger DA. Logic-based models for the analysis of cell signaling networks. Biochemistry. 2010 Apr;49(15):3216–24. pmid:20225868
  37. 37. Albert R, Thakar J. Boolean modeling: a logic-based dynamic approach for understanding signaling and regulatory networks and for making useful predictions. Wiley Interdiscip Rev Syst Biol Med. 2014;6(5):353–69. pmid:25269159
  38. 38. Eduati F, Doldàn-Martelli V, Klinger B, Cokelaer T, Sieber A, Kogera F, et al. Drug Resistance Mechanisms in Colorectal Cancer Dissected with Cell Type-Specific Dynamic Logic Models. Cancer Res. 2017 Jun;77(12):3364–75. pmid:28381545
  39. 39. Kraeutler MJ, Soltis AR, Saucerman JJ. Modeling cardiac β-adrenergic signaling with normalized-Hill differential equations: comparison with a biochemical model. BMC Syst Biol. 2010 Nov;4:157. pmid:21087478
  40. 40. Eduati F, Jaaks P, Wappler J, Cramer T, Merten CA, Garnett MJ, et al. Patient-specific logic models of signaling pathways from screenings on cancer biopsies to prioritize personalized combination therapies. Mol Syst Biol. 2020 Feb;16(2):e8664. pmid:32073727
  41. 41. Henriques D, Villaverde AF, Rocha M, Saez-Rodriguez J, Banga JR. Data-driven reverse engineering of signaling pathways using ensembles of dynamic models. PLoS Comput Biol. 2017 Feb;13(2):e1005379. pmid:28166222
  42. 42. Frank DU, Sutcliffe MD, Saucerman JJ. Network-based predictions of in vivo cardiac hypertrophy. J Mol Cell Cardiol. 2018 Aug;121:180–9. pmid:30030017
  43. 43. Ugray Z, Lasdon L, Plummer J, Glover F, Kelly J, Martí R. Scatter Search and Local NLP Solvers: A Multistart Framework for Global Optimization. INFORMS J Comput [Internet]. 2007;19(3):328–40. Available from: https://doi.org/10.1287/ijoc.1060.0175
  44. 44. Khare YP, Muñoz-Carpena R, Rooney RW, Martinez CJ. A multi-criteria trajectory-based parameter sampling strategy for the screening method of elementary effects. Environ Model Softw. 2015;64:230–9.
  45. 45. Sobol IM. Global sensitivity indices for nonlinear mathematical models and their Monte Carlo estimates. Math Comput Simul. 2001;55(1–3):271–80.
  46. 46. Saltelli A, Annoni P, Azzini I, Campolongo F, Ratto M, Tarantola S. Variance based sensitivity analysis of model output. Design and estimator for the total sensitivity index. Comput Phys Commun [Internet]. 2010;181(2):259–70. Available from: http://www.sciencedirect.com/science/article/pii/S0010465509003087
  47. 47. Janes KA. An analysis of critical factors for quantitative immunoblotting. Sci Signal. 2015 Apr;8(371):rs2. pmid:25852189
  48. 48. Freund C, Schmidt-Ullrich R, Baurand A, Dunger S, Schneider W, Loser P, et al. Requirement of nuclear factor-kappaB in angiotensin II- and isoproterenol-induced cardiac hypertrophy in vivo. Circulation. 2005 May;111(18):2319–25. pmid:15870116
  49. 49. Weeks KL, McMullen JR. The athlete’s heart vs. the failing heart: can signaling explain the two distinct outcomes? Physiology (Bethesda). 2011 Apr;26(2):97–105. pmid:21487028
  50. 50. Lohse MJ, Engelhardt S, Eschenhagen T. What is the role of beta-adrenergic signaling in heart failure? Circ Res. 2003 Nov;93(10):896–906. pmid:14615493
  51. 51. Zhang G-X, Ohmori K, Nagai Y, Fujisawa Y, Nishiyama A, Abe Y, et al. Role of AT1 receptor in isoproterenol-induced cardiac hypertrophy and oxidative stress in mice. J Mol Cell Cardiol. 2007 Apr;42(4):804–11. pmid:17350036
  52. 52. Vidal M, Wieland T, Lohse MJ, Lorenz K. β-Adrenergic receptor stimulation causes cardiac hypertrophy via a Gβγ/Erk-dependent pathway. Cardiovasc Res. 2012 Nov;96(2):255–64. pmid:22843704
  53. 53. Yuan W, Bers DM. Ca-dependent facilitation of cardiac Ca current is due to Ca-calmodulin-dependent protein kinase. Am J Physiol. 1994 Sep;267(3 Pt 2):H982–93. pmid:8092302
  54. 54. Hudmon A, Schulman H, Kim J, Maltez JM, Tsien RW, Pitt GS. CaMKII tethers to L-type Ca2+ channels, establishing a local and dedicated integrator of Ca2+ signals for facilitation. J Cell Biol. 2005 Nov;171(3):537–47. pmid:16275756
  55. 55. Grueter CE, Colbran RJ, Anderson ME. CaMKII, an emerging molecular driver for calcium homeostasis, arrhythmias, and cardiac dysfunction. J Mol Med (Berl). 2007 Jan;85(1):5–14.
  56. 56. Lee T-S, Karl R, Moosmang S, Lenhardt P, Klugbauer N, Hofmann F, et al. Calmodulin kinase II is involved in voltage-dependent facilitation of the L-type Cav1.2 calcium channel: Identification of the phosphorylation sites. J Biol Chem. 2006 Sep;281(35):25560–7. pmid:16820363
  57. 57. Grimm M, Brown JH. Beta-adrenergic receptor signaling in the heart: role of CaMKII. J Mol Cell Cardiol. 2010 Feb;48(2):322–30. pmid:19883653
  58. 58. Soltis AR, Saucerman JJ. Synergy between CaMKII substrates and β-adrenergic signaling in regulation of cardiac myocyte Ca(2+) handling. Biophys J. 2010 Oct;99(7):2038–47. pmid:20923637
  59. 59. Mishra S, Ling H, Grimm M, Zhang T, Bers DM, Brown JH. Cardiac hypertrophy and heart failure development through Gq and CaM kinase II signaling. J Cardiovasc Pharmacol. 2010 Dec;56(6):598–603. pmid:20531218
  60. 60. Molkentin JD. Calcineurin-NFAT signaling regulates the cardiac hypertrophic response in coordination with the MAPKs. Cardiovasc Res. 2004 Aug;63(3):467–75. pmid:15276472
  61. 61. Lu Y-M, Shioda N, Han F, Kamata A, Shirasaki Y, Qin Z-H, et al. DY-9760e inhibits endothelin-1-induced cardiomyocyte hypertrophy through inhibition of CaMKII and ERK activities. Cardiovasc Ther. 2009;27(1):17–27. pmid:19207476
  62. 62. Cipolletta E, Rusciano MR, Maione AS, Santulli G, Sorriento D, Del Giudice C, et al. Targeting the CaMKII/ERK Interaction in the Heart Prevents Cardiac Hypertrophy. PLoS One. 2015;10(6):e0130477. pmid:26110816
  63. 63. Zou Y, Yao A, Zhu W, Kudoh S, Hiroi Y, Shimoyama M, et al. Isoproterenol activates extracellular signal-regulated protein kinases in cardiomyocytes through calcineurin. Circulation. 2001 Jul;104(1):102–8. pmid:11435346
  64. 64. Sanna B, Bueno OF, Dai Y-S, Wilkins BJ, Molkentin JD. Direct and indirect interactions between calcineurin-NFAT and MEK1-extracellular signal-regulated kinase 1/2 signaling pathways regulate cardiac gene expression and cellular growth. Mol Cell Biol. 2005 Feb;25(3):865–78. pmid:15657416
  65. 65. Tebar F, Villalonga P, Sorkina T, Agell N, Sorkin A, Enrich C. Calmodulin regulates intracellular trafficking of epidermal growth factor receptor and the MAPK signaling pathway. Mol Biol Cell. 2002 Jun;13(6):2057–68. pmid:12058069
  66. 66. Salzano M, Rusciano MR, Russo E, Bifulco M, Postiglione L, Vitale M. Calcium/calmodulin-dependent protein kinase II (CaMKII) phosphorylates Raf-1 at serine 338 and mediates Ras-stimulated Raf-1 activation. Cell Cycle. 2012 Jun;11(11):2100–6. pmid:22592532
  67. 67. Grisanti LA, Guo S, Tilley DG. Cardiac GPCR-Mediated EGFR Transactivation: Impact and Therapeutic Implications. J Cardiovasc Pharmacol. 2017 Jul;70(1):3–9. pmid:28059858
  68. 68. Kehat I, Molkentin JD. Extracellular signal-regulated kinase 1/2 (ERK1/2) signaling in cardiac hypertrophy. Ann N Y Acad Sci. 2010 Feb;1188:96–102. pmid:20201891
  69. 69. Gallo S, Vitacolonna A, Bonzano A, Comoglio P, Crepaldi T. ERK: A Key Player in the Pathophysiology of Cardiac Hypertrophy. Int J Mol Sci. 2019 May;20(9). pmid:31052420
  70. 70. Robin X, Creixell P, Radetskaya O, Santini CC, Longden J, Linding R. Personalized network-based treatments in oncology. Clin Pharmacol Ther. 2013 Dec;94(6):646–50. pmid:23995267
  71. 71. Saez-Rodriguez J, Blüthgen N. Personalized signaling models for personalized treatments. Mol Syst Biol. 2020 Jan;16(1):e9042. pmid:32129942
  72. 72. Razzaq M, Paulevé L, Siegel A, Saez-Rodriguez J, Bourdon J, Guziolowski C. Computational discovery of dynamic cell line specific Boolean networks from multiplex time-course data. PLoS Comput Biol. 2018 Oct;14(10):e1006538. pmid:30372442
  73. 73. Zhang R, Shah MV, Yang J, Nyland SB, Liu X, Yun JK, et al. Network model of survival signaling in large granular lymphocyte leukemia. Proc Natl Acad Sci U S A. 2008 Oct;105(42):16308–13. pmid:18852469
  74. 74. Saadatpour A, Wang R-S, Liao A, Liu X, Loughran TP, Albert I, et al. Dynamical and structural analysis of a T cell survival network identifies novel candidate therapeutic targets for large granular lymphocyte leukemia. PLoS Comput Biol. 2011 Nov;7(11):e1002267. pmid:22102804
  75. 75. Thakar J, Saadatpour-Moghaddam A, Harvill ET, Albert R. Constraint-based network model of pathogen-immune system interactions. J R Soc Interface. 2009 Jul;6(36):599–612. pmid:18952547
  76. 76. Maheshwari P, Du H, Sheen J, Assmann SM, Albert R. Model-driven discovery of calcium-related protein-phosphatase inhibition in plant guard cell signaling. PLoS Comput Biol. 2019 Oct;15(10):e1007429. pmid:31658257
  77. 77. Stolovitzky G, Monroe D, Califano A. Dialogue on reverse-engineering assessment and methods: the DREAM of high-throughput pathway inference. Ann N Y Acad Sci. 2007 Dec;1115:1–22. pmid:17925349
  78. 78. Watkins SJ, Borthwick GM, Arthur HM. The H9C2 cell line and primary neonatal cardiomyocyte cells show similar hypertrophic responses in vitro. In Vitro Cell Dev Biol Anim. 2011 Feb;47(2):125–31. pmid:21082279
  79. 79. Branco AF, Pereira SL, Moreira AC, Holy J, Sardão VA, Oliveira PJ. Isoproterenol cytotoxicity is dependent on the differentiation state of the cardiomyoblast H9c2 cell line. Cardiovasc Toxicol. 2011 Sep;11(3):191–203. pmid:21455642
  80. 80. Illario M, Cavallo AL, Bayer KU, Di Matola T, Fenzi G, Rossi G, et al. Calcium/calmodulin-dependent protein kinase II binds to Raf-1 and modulates integrin-stimulated ERK activation. J Biol Chem. 2003 Nov;278(46):45101–8. pmid:12954639
  81. 81. De Windt LJ, Lim HW, Haq S, Force T, Molkentin JD. Calcineurin promotes protein kinase C and c-Jun NH2-terminal kinase activation in the heart. Cross-talk between cardiac hypertrophic signaling pathways. J Biol Chem. 2000 May;275(18):13571–9. pmid:10788473
  82. 82. Hudson TJ, Anderson W, Artez A, Barker AD, Bell C, Bernabé RR, et al. International network of cancer genome projects. Nature. 2010 Apr;464(7291):993–8. pmid:20393554
  83. 83. Jeske L, Placzek S, Schomburg I, Chang A, Schomburg D. BRENDA in 2019: a European ELIXIR core data resource. Nucleic Acids Res. 2019 Jan;47(D1):D542–9. pmid:30395242
  84. 84. Krempel R, Kulkarni P, Yim A, Lang U, Habermann B, Frommolt P. Integrative analysis and machine learning on cancer genomics data using the Cancer Systems Biology Database (CancerSysDB). BMC Bioinformatics. 2018 Apr;19(1):156. pmid:29699486