Attractor landscape analysis of the cardiac signaling network reveals mechanism-based therapeutic strategies for heart failure

Attractor landscape analysis of the cardiac signaling network reveals mechanism-based therapeutic... Abstract Apoptosis and hypertrophy of cardiomyocytes are the primary causes of heart failure (HF), a global leading cause of death, and are regulated through the complicated intracellular signaling network, limiting the development of effective treatments due to its complexity. To identify effective therapeutic strategies for HF at a system level, we develop a large-scale comprehensive mathematical model of the cardiac signaling network by integrating all available experimental evidence. Attractor landscape analysis of the network model identifies distinct sets of control nodes that effectively suppress apoptosis and hypertrophy of cardiomyocytes under ischemic or pressure overload-induced HF, the two major types of HF. Intriguingly, our system-level analysis suggests that intervention of these control nodes may increase the efficacy of clinical drugs for HF and, of most importance, different combinations of control nodes are suggested as potentially effective candidate drug targets depending on the types of HF. Our study provides a systematic way of developing mechanism-based therapeutic strategies for HF. heart failure, mathematical modeling, cardiac signaling network, systems analysis, disease mechanism, mechanism-based therapeutic strategy, systems biology Introduction Heart failure (HF) is a clinical syndrome characterized by typical symptoms such as shortness of breathing and peripheral swelling resulting from a reduced cardiac output (Remme and Swedberg, 2001; Hilfiker-Klemer et al., 2006). The disease is caused by a structural and/or functional cardiac abnormality that arises from myocardial ischemia (e.g. angina, myocardial infarction) or sustained pressure overload (e.g. hypertension) (Remme and Swedberg, 2001; Hilfiker-Klemer et al., 2006). It is one of the global leading causes of death, and its incidence and prevalence are increasing over time (Stewart et al., 2003). To date, several drugs are available for HF treatment that are effective at an early stage in clinics (ACE Inhibitor Myocardial Infarction Collaborative Group, 1998), but yet they are reported to have limited efficacy in blocking the disease progression in the long-term treatment, which highlights the urgent need to develop effective therapeutic strategies. Apoptosis and hypertrophy of cardiomyocytes are the main causes of HF and they are the typically observed molecular phenotypes in cardiomyocytes during the disease progression (Remme and Swedberg, 2001). Cardiomyocyte apoptosis occurs in diverse pathological conditions, including ischemia or pressure overload, and causes the permanent loss of the functioning units in the myocardium, thereby contributing to the development of HF (Narula et al., 1996; Gonzalez et al., 2003). Meanwhile, cardiomyocyte hypertrophy occurs to compensate for the loss of myocytes following injury and to increase the cardiac pumping function temporally (Kang et al., 2016). However, over time, myocardial hypertrophy induces pathological changes in the heart that eventually lead to a functional deterioration, which also promotes HF (Nadal-Ginard et al., 2003). Apoptosis and hypertrophy of cardiomyocytes are known to be regulated by the complicated interactions in the intracellular cardiac signaling network (Ryall et al., 2012; Kang et al., 2016). Previous experimental researches were conducted to understand how the cardiac signal transduction pathways are regulated and how they are related to pathological phenotypes of cardiomyocytes (Saraste et al., 1997; Bernardo et al., 2010). However, there is a fundamental limitation in applying conventional experimental approaches to investigate the underlying mechanisms of HF because the cardiac signaling network involves complicated feedback regulation and crosstalk (Kitano, 2002), which makes it difficult to intuitively understand the cellular behavior and develop effective treatments for HF patients. To overcome such limitations and investigate effective therapeutic strategies for HF treatment, we have employed systems biological approaches based on mathematical modeling and simulation analysis that have been proven to facilitate the investigation of the complex underlying mechanisms of cellular signaling networks (Kitano, 2002; Kim et al., 2007; Murray et al., 2010). To this end, we have constructed a very large-scale comprehensive cardiac signaling network and developed its mathematical model by integrating all relevant experimental information that is closely associated with cardiomyocytes. Attractor landscape analysis, which aims to understand a cellular behavior by investigating a specific converging network state determined by the interaction of signaling molecules, has been widely used to explore the dynamics of a complex signaling network (Choi et al., 2012; Cho et al., 2016). From attractor landscape analysis of the cardiac signaling network model, we have identified distinct sets of control nodes that can effectively suppress apoptosis and hypertrophy of cardiomyocytes under ischemic or pressure overload-induced HF condition, the two major types of HF. Intriguingly, we found that intervention of these control nodes may increase the efficacy of currently available clinical drugs for HF patients. Most importantly, different combinations of control nodes are suggested as potentially effective candidate drug targets depending on the types of HF. Our study highlights the significance of network-level analyses for unraveling molecular regulatory mechanisms, and provides a systematic way of developing mechanism-based therapeutic strategies for HF. Results A mathematical model of the cardiac signaling network To understand the complicated dynamics of the underlying intracellular signaling process of cardiomyocytes, we first reconstructed a very large-scale comprehensive cardiac signaling network by integrating all relevant evidence on molecules and their interactions through an extensive search over literature and database such as Kyoto Encyclopedia of Genes and Genomes (KEGG) and Pathway Interaction Database (PID) (Kanehisa and Goto, 2000; Schaefer et al., 2009). The reconstructed cardiac signaling network, the most comprehensive cardiac signaling network to the best of our knowledge at present, consists of 141 nodes and 274 directed links, including 17 external input nodes and 15 output nodes (Figure 1 and Supplementary Dataset S1). The input nodes in the cardiac signaling network represent diverse external stimuli that activate downstream signaling pathways to transmit signals from the outer surface of cardiomyocytes to cell interior and therefore generate appropriate cellular responses (Figure 1 and Supplementary Dataset S1). The network involves various signaling pathways, including α- or β-adrenergic (AR) (Kanehisa and Goto, 2000), calcium (Ca) (Kehat and Molkentin, 2010), cyclic guanosine monophosphate (cGMP)-protein kinase G (PKG) (Francis, 2010), mitogen-activated protein kinase (MAPK) (Wang, 2007), phosphatidylinositol 3-kinase (PI3K)-Akt (Matsui and Rosenzweig, 2005), transforming growth factor beta (TGFβ) (Bujak and Frangogiannis, 2007), Src-focal adhesion kinase (FAK) (Kuramochi et al., 2006), Janus kinase (JAK)-signal transducer and activator of transcription (STAT) (Pan et al., 1999), tumor necrosis factor alpha (TNFα) (Condorelli et al., 2002), and nuclear factor kappa-light-chain-enhancer of activated B cells (NFκB) pathways (Ling et al., 2013), which play important roles in coordinating the cellular functions of cardiomyocytes. These signaling pathways relay and integrate biochemical signals, and transmit them to the nucleus to activate enzymes and transcription factors that regulate apoptosis or hypertrophy of cardiomyocytes. In the cardiac signaling network, caspase 3 (Casp3) was used as an indicator of apoptosis, while the remaining 14 output nodes were taken as indicators of hypertrophy on the basis of their crucial roles for the regulation of major cellular functions in cardiomyocytes (Figure 1 and Supplementary Table S1). Of note, several nodes that are involved multiple signaling pathways function to mediate crosstalks between different cascades and they are known to play a key role in coordinating cellular responses of cardiomyocytes under diverse pathophysiological conditions (Heineke and Molkentin, 2006; Kang et al., 2016) (Figure 1). We further found that the cardiac signaling network has a higher clustering coefficient and neighborhood connectivity than random networks (Supplementary Table S2) and that the in- and out-degree distributions of nodes in the network follow an approximate power-law distribution (Supplementary Figure S1), which are in accord with the general characteristics of biological networks (Watts and Strogatz, 1998; Barabasi and Oltvai, 2004; Barabasi et al., 2011; Ryall et al., 2012; Lu et al., 2015). Figure 1 View largeDownload slide A schematic diagram of the cardiac signaling network. The reconstructed cardiac signaling network consists of 141 nodes and 274 links; 212 links are activation links (red pointed arrows) and 62 links are inhibition links (blue blunted arrows). Of the 141 nodes, there are 17 input nodes and 15 output nodes. Nodes that are involved in two or more signaling pathways according to the KEGG classification are highlighted with multiple colors that represent the corresponding pathways. Detailed information for each link and assigned logic table for each node is shown in Supplementary Datasets S1 and S2. Figure 1 View largeDownload slide A schematic diagram of the cardiac signaling network. The reconstructed cardiac signaling network consists of 141 nodes and 274 links; 212 links are activation links (red pointed arrows) and 62 links are inhibition links (blue blunted arrows). Of the 141 nodes, there are 17 input nodes and 15 output nodes. Nodes that are involved in two or more signaling pathways according to the KEGG classification are highlighted with multiple colors that represent the corresponding pathways. Detailed information for each link and assigned logic table for each node is shown in Supplementary Datasets S1 and S2. To validate whether our network model components are all expressed in cardiomyocytes, we have further evaluated the gene or protein expression levels of the network components using mRNA-sequencing or biochemical experimental data obtained from previous studies (Nakamura et al., 2000; Cuenca et al., 2006; Yang et al., 2006; Loan Le et al., 2012; Song et al., 2012; Ucar et al., 2012; Zhao et al., 2016). Note that for the case of mRNA-sequencing data analysis, genes whose reads per kilobase per million mapped reads (RPKM) values are greater than 1 were determined to be expressed as in the previous studies (Phanstiel et al., 2011; Handler et al., 2013). As a result, we found that all the nodes (except input ligands that are derived from extracellular microenvironment of cardiomyocytes, and the chemical elements that are not transcriptionally or translationally expressed/activated such as Ca, cyclic adenosine monophosphate (cAMP), cGMP, and inositol trisphosphate (IP3)) are specifically expressed in cardiomyocytes (Supplementary Table S3), further supporting that the reconstructed network can comprehensively represent cardiomyocyte-specific signaling mechanisms. Due to the limited experimental information on detailed reaction kinetics of the interactions between biomolecules, continuous dynamic models (e.g. an ordinary differential equation-based model) have difficulties in estimating their parameter values (Park et al., 2006; Helikar et al., 2008). To avoid such a problem, parameter-free discrete Boolean network modeling approach has been widely used in modeling various biological systems, which is proven to be highly useful in investigating key dynamic properties of complex intracellular signaling networks (Helikar et al., 2008). This modeling method has been shown to successfully capture the well-known biological features and provides reliable predictions of the potential effects of drug treatments (Choi et al., 2012; Lee et al., 2015). Therefore, we have developed a large-scale Boolean model of the cardiac signaling network based on the qualitative properties of activation and/or inhibition of the signaling proteins in cardiomyocytes (Supplementary Dataset S2). In the cardiac signaling network model, the state value of a distinct node represents its activity, and is reflected as ‘1’ or ‘0’ for each active or inactive state, respectively. Logic tables for nodes were created based on experimental information from the relevant literature and database search (Supplementary Datasets S1 and S2). To investigate the optimal input setting that can reasonably represent the physiological condition of cardiomyocytes, we have performed extensive simulations for the analysis of input–output relationships of network components (see Materials and methods). To this end, 100 million sets of randomly sampled input intensities were applied to the network and the resulting cellular responses were monitored (see Materials and methods). As a result, the intensity of input nodes that can most closely reproduce the well-known qualitative activities of network components was determined as an optimal input setting for the physiological condition (Supplementary Table S4). The simulation results of the input–output relationship analysis are shown in Figure 2, Supplementary Figures S2 and S3, and Table S5, and these results indicate that the cardiac signaling network model can well replicate the cellular responses of cardiomyocytes that are experimentally observed (Aikawa et al., 1997; Condorelli et al., 2002; Chen et al., 2003; Nakaoka et al., 2003, 2007; Fu et al., 2004; Vega et al., 2004; Zarubin and Han, 2005; Schaub and Hefti, 2007; Colella et al., 2008; Fischer and Hilfiker-Kleiner, 2008; Kawasaki et al., 2008; Guo and Wang, 2009; Ronca et al., 2009; Backs et al., 2011; Liu et al., 2012; Pennanen et al., 2014; Roberge et al., 2014; Wang et al., 2015; Ponnusamy et al., 2017). Figure 2 View largeDownload slide Qualitative input–output relationships in the cardiac signaling network model using last 100 average steps for the activity calculation. (A) Positive relationship between TNFα and Casp8 activation (Roberge et al., 2014). (B) Positive relationship between TNFα and p38 activation (Chen et al., 2003). (C) Positive relationship between TNFα and ANP expression (Condorelli et al., 2002). (D) Positive relationship between TNFα and cJun activation (Condorelli et al., 2002). (E) Positive relationship between NE and Ca activation (Pennanen et al., 2014). (F) Positive relationship between NE and NFAT activation (Colella et al., 2008). (G) Positive relationship between NE and MEF2 activation (Backs et al., 2011). (H) Positive relationship between NE and Casp3 activation (Fu et al., 2004). (I) Positive relationship between NRG1 and SHP2 activation (Nakaoka et al., 2007). (J) Positive relationship between PGE2 and PKA activation (Liu et al., 2012). (K) Positive relationship between LIF and GAB1 activation (Nakaoka et al., 2003). (L) Positive relationship between LIF and STAT activation (Schaub and Hefti, 2007). (M) Positive relationship between TGFβ and Ras activation (Guo and Wang, 2009). (N) Positive relationship between TGFβ and p38 activation (Zarubin and Han, 2005). (O) Positive relationship between TGFβ and BNP expression (Vega et al., 2004). (P) Positive relationship between TGFβ and ANP expression (Vega et al., 2004). (Q) Positive relationship between FGF and Ras activation (Ronca et al., 2009). (R) Positive relationship between AngII and Src activation (Aikawa et al., 1997). (S) Positive relationship between CT1 and JAK activation (Fischer and Hilfiker-Kleiner, 2008). (T) Positive relationship between CT1 and Ras activation (Fischer and Hilfiker-Kleiner, 2008). (U) Positive relationship between HGF and SHP2 activation (Wang et al., 2015). (V) Positive relationship between VEGF and Ras activation (Kawasaki et al., 2008). (W) Negative relationship between FGF and FOXO activation (Ponnusamy et al., 2017). (X) Positive relationship between EGF and GAB1 activation (Nakaoka et al., 2007). The average ‘1’s over last 100 time steps were calculated and used for stable state activity on each node. The dose–response curve presented here is intended to illustrate how the cardiac signaling network model can qualitatively reproduce the referenced input–output relationship for a wide range of input signals. The simulation was repeated (n = 10) for the input at 2% noise (see Materials and methods). Figure 2 View largeDownload slide Qualitative input–output relationships in the cardiac signaling network model using last 100 average steps for the activity calculation. (A) Positive relationship between TNFα and Casp8 activation (Roberge et al., 2014). (B) Positive relationship between TNFα and p38 activation (Chen et al., 2003). (C) Positive relationship between TNFα and ANP expression (Condorelli et al., 2002). (D) Positive relationship between TNFα and cJun activation (Condorelli et al., 2002). (E) Positive relationship between NE and Ca activation (Pennanen et al., 2014). (F) Positive relationship between NE and NFAT activation (Colella et al., 2008). (G) Positive relationship between NE and MEF2 activation (Backs et al., 2011). (H) Positive relationship between NE and Casp3 activation (Fu et al., 2004). (I) Positive relationship between NRG1 and SHP2 activation (Nakaoka et al., 2007). (J) Positive relationship between PGE2 and PKA activation (Liu et al., 2012). (K) Positive relationship between LIF and GAB1 activation (Nakaoka et al., 2003). (L) Positive relationship between LIF and STAT activation (Schaub and Hefti, 2007). (M) Positive relationship between TGFβ and Ras activation (Guo and Wang, 2009). (N) Positive relationship between TGFβ and p38 activation (Zarubin and Han, 2005). (O) Positive relationship between TGFβ and BNP expression (Vega et al., 2004). (P) Positive relationship between TGFβ and ANP expression (Vega et al., 2004). (Q) Positive relationship between FGF and Ras activation (Ronca et al., 2009). (R) Positive relationship between AngII and Src activation (Aikawa et al., 1997). (S) Positive relationship between CT1 and JAK activation (Fischer and Hilfiker-Kleiner, 2008). (T) Positive relationship between CT1 and Ras activation (Fischer and Hilfiker-Kleiner, 2008). (U) Positive relationship between HGF and SHP2 activation (Wang et al., 2015). (V) Positive relationship between VEGF and Ras activation (Kawasaki et al., 2008). (W) Negative relationship between FGF and FOXO activation (Ponnusamy et al., 2017). (X) Positive relationship between EGF and GAB1 activation (Nakaoka et al., 2007). The average ‘1’s over last 100 time steps were calculated and used for stable state activity on each node. The dose–response curve presented here is intended to illustrate how the cardiac signaling network model can qualitatively reproduce the referenced input–output relationship for a wide range of input signals. The simulation was repeated (n = 10) for the input at 2% noise (see Materials and methods). Reflecting molecular characteristics under HF conditions to the cardiac signaling network model The aim of this study is to investigate most effective therapeutic strategies for HF depending on conditions. Although HF is caused by multiple pathologies (Hilfiker-Klemer et al., 2006), myocardial ischemia or sustained pressure overload is known to be its leading cause according to previous studies and clinical evidence (Blankesteijn et al., 2008; Thygesen et al., 2012). For this reason, we focused on investigating the signaling mechanisms of the two major types of HF (i.e. ischemic HF and pressure overload-induced HF). The heart function is regulated by various growth factors (i.e. insulin-like growth factor 1 (IGF1), epidermal growth factor (EGF), vascular endothelial growth factor (VEGF), TGFβ, and fibroblast growth factor (FGF)), catecholamines (i.e. norepinephrine (NE) and phenylephrine (PE)) or cytokines (i.e. TNFα, endothelin-1 (ET1), and interleukin (IL)) of which the secretion is modulated by autonomic nerve or endocrine system (Iwamoto et al., 2003; Florea and Cohn, 2014). During the development of HF due to the cardiac dysfunction, the sympathetic nervous system or the adrenal gland is activated as an early compensatory mechanism that supports inotropic action and maintains cardiac output (Jackson et al., 2000). This activated sympathetic nervous system or the adrenal gland changes the secretion profiles of various growth factors, catecholamines, or cytokines, which show different aspects from those of physiological condition (Kang et al., 2016). The changes in the secretion of those extracellular stimuli augment cardiac contractility and relaxation that increase the heart pumping function in the early phase of HF (Nadal-Ginard et al., 2003). However, in the long run, these stimuli activate various signaling pathways that play pivotal roles in the regulation of cardiomyocyte apoptosis and hypertrophy, which eventually contribute to the progression of HF (MacLellan and Schneider, 2000). Considering that the change in stimuli secretion profiles is one of the major causes of cardiomyocyte apoptosis and hypertrophy (MacLellan and Schneider, 2000), we have reflected ischemic or pressure overload-induced HF condition to the cardiac signaling network model by applying the distinct input intensities of individual HF condition. For this purpose, we integrated all relevant information on upregulated or downregulated molecular activities in each HF condition through an extensive search over literature and database (Supplementary Table S6). These collected observations include not only the changed activity profiles of the external ligands but also those of intracellular signaling molecules and transcription factors under HF conditions. Based on the fact that the changed stimuli secretion profiles in HF conditions cause the changes in the activity of downstream molecules (MacLellan and Schneider, 2000), we have performed extensive simulation analyses while randomly increasing (or decreasing) the intensity of input nodes of which the activity is enhanced (or suppressed) under distinct HF conditions in a wide range for 100 million times (Supplementary Table S6), and monitored the activity profiles of downstream signaling components in the cardiac signaling network. As a result, we found the optimal input setting that can reasonably represent distinct HF conditions of cardiomyocytes (Figure 3). To determine the robustness of the optimized input setting for HF conditions, additional computer simulations were performed while randomly perturbing the intensity of upregulated or downregulated inputs in the range of 5% or 20% for 100 million times. As a result, activity profiles of the molecules were well reproduced regardless of the perturbation range (Supplementary Figure S4). This result indicates that the mathematical model is adequate for studying the molecular mechanisms of the cardiac signaling network under HF conditions. Figure 3 View largeDownload slide The identified optimal input setting representing ischemic or pressure overload-induced HF condition of cardiomyocytes. The optimal input setting of cardiomyocytes was determined for ischemic or pressure overload-induced HF condition that reasonably reproduces qualitative characteristics of the activity profile for the network components under distinct HF conditions. Detailed information about the optimal input setting is shown in the main text. Blue or red nodes indicate differentially activated proteins (cutoffs: fold change >1.5, P < 0.05) under ischemic or pressure overload-induced HF, respectively. Green nodes indicate upregulated components in signaling pathways that contain the differentially activated proteins. Figure 3 View largeDownload slide The identified optimal input setting representing ischemic or pressure overload-induced HF condition of cardiomyocytes. The optimal input setting of cardiomyocytes was determined for ischemic or pressure overload-induced HF condition that reasonably reproduces qualitative characteristics of the activity profile for the network components under distinct HF conditions. Detailed information about the optimal input setting is shown in the main text. Blue or red nodes indicate differentially activated proteins (cutoffs: fold change >1.5, P < 0.05) under ischemic or pressure overload-induced HF, respectively. Green nodes indicate upregulated components in signaling pathways that contain the differentially activated proteins. To further investigate whether the signaling dynamics are differentially regulated under ischemic or pressure overload-induced HF condition, we performed pathway analysis of the signaling proteins that are significantly differentially activated under each HF condition using g:Profiler (Reimand et al., 2016). As a result, the differentially activated proteins (cutoffs: fold change >1.5, P < 0.05) between each individual HF type were found to be significantly enriched in different signaling pathways (Figure 3 and Table 1). These simulation results are consistent with previous experiments: for instance, adenylyl cyclase (AC)-cAMP signaling promotes the progression of ischemic HF (Appukuttan et al., 2012; Bravo et al., 2016), while TGFβ signaling is found to be important in the development of pressure overload-induced HF (Lim et al., 2005; Koitabashi et al., 2011). This further supports that the underlying signaling mechanisms may differ depending on the HF types and, therefore, different signaling pathways should be targeted to achieve optimal therapeutic responses, which can be useful for the development of mechanism-based treatment strategy for HF. Table 1 Pathway analysis of differentially activated proteins between ischemic and pressure overload-induced HF. Category Term Ischemic HF Pressure overload-induced HF Count Percent P-value Count Percent P-value KEGG Adrenergic signaling in cardiomyocytes 10 6.9 <0.001 0 0.0 N.S. Calcium signaling pathway 9 5.0 <0.001 1 0.6 N.S. cAMP signaling pathway 8 4.0 <0.001 1 0.5 N.S. cGMP–PKG signaling pathway 8 4.9 <0.001 2 1.2 N.S. Jak–STAT signaling pathway 5 3.2 0.020 0 0.0 N.S. Ras signaling pathway 6 2.6 0.013 3 1.3 N.S Renin–angiotensin system 3 13 0.006 0 0.0 N.S. Focal adhesion 2 1.0 N.S. 6 3.0 <0.001 FoxO signaling pathway 0 0.0 N.S. 5 3.8 <0.001 MAPK signaling pathway 5 2.0 N.S. 6 2.4 <0.001 TGF-β signaling pathway 1 1.2 N.S. 6 7.1 <0.001 Category Term Ischemic HF Pressure overload-induced HF Count Percent P-value Count Percent P-value KEGG Adrenergic signaling in cardiomyocytes 10 6.9 <0.001 0 0.0 N.S. Calcium signaling pathway 9 5.0 <0.001 1 0.6 N.S. cAMP signaling pathway 8 4.0 <0.001 1 0.5 N.S. cGMP–PKG signaling pathway 8 4.9 <0.001 2 1.2 N.S. Jak–STAT signaling pathway 5 3.2 0.020 0 0.0 N.S. Ras signaling pathway 6 2.6 0.013 3 1.3 N.S Renin–angiotensin system 3 13 0.006 0 0.0 N.S. Focal adhesion 2 1.0 N.S. 6 3.0 <0.001 FoxO signaling pathway 0 0.0 N.S. 5 3.8 <0.001 MAPK signaling pathway 5 2.0 N.S. 6 2.4 <0.001 TGF-β signaling pathway 1 1.2 N.S. 6 7.1 <0.001 Table showing the pathway analysis results of differentially activated proteins (cutoffs: fold change >1.5, P < 0.05) between ischemic and pressure overload-induced HF using g:Profiler (Reimand et al., 2016). N.S., no significance. Table 1 Pathway analysis of differentially activated proteins between ischemic and pressure overload-induced HF. Category Term Ischemic HF Pressure overload-induced HF Count Percent P-value Count Percent P-value KEGG Adrenergic signaling in cardiomyocytes 10 6.9 <0.001 0 0.0 N.S. Calcium signaling pathway 9 5.0 <0.001 1 0.6 N.S. cAMP signaling pathway 8 4.0 <0.001 1 0.5 N.S. cGMP–PKG signaling pathway 8 4.9 <0.001 2 1.2 N.S. Jak–STAT signaling pathway 5 3.2 0.020 0 0.0 N.S. Ras signaling pathway 6 2.6 0.013 3 1.3 N.S Renin–angiotensin system 3 13 0.006 0 0.0 N.S. Focal adhesion 2 1.0 N.S. 6 3.0 <0.001 FoxO signaling pathway 0 0.0 N.S. 5 3.8 <0.001 MAPK signaling pathway 5 2.0 N.S. 6 2.4 <0.001 TGF-β signaling pathway 1 1.2 N.S. 6 7.1 <0.001 Category Term Ischemic HF Pressure overload-induced HF Count Percent P-value Count Percent P-value KEGG Adrenergic signaling in cardiomyocytes 10 6.9 <0.001 0 0.0 N.S. Calcium signaling pathway 9 5.0 <0.001 1 0.6 N.S. cAMP signaling pathway 8 4.0 <0.001 1 0.5 N.S. cGMP–PKG signaling pathway 8 4.9 <0.001 2 1.2 N.S. Jak–STAT signaling pathway 5 3.2 0.020 0 0.0 N.S. Ras signaling pathway 6 2.6 0.013 3 1.3 N.S Renin–angiotensin system 3 13 0.006 0 0.0 N.S. Focal adhesion 2 1.0 N.S. 6 3.0 <0.001 FoxO signaling pathway 0 0.0 N.S. 5 3.8 <0.001 MAPK signaling pathway 5 2.0 N.S. 6 2.4 <0.001 TGF-β signaling pathway 1 1.2 N.S. 6 7.1 <0.001 Table showing the pathway analysis results of differentially activated proteins (cutoffs: fold change >1.5, P < 0.05) between ischemic and pressure overload-induced HF using g:Profiler (Reimand et al., 2016). N.S., no significance. Attractor landscape analysis of the cardiac signaling network model To investigate the dynamic behavior of the intracellular signaling process of cardiomyocytes and to identify effective molecular targets depending on the types of HF, we applied the state space analysis to the cardiac signaling network model, which has been proved to be effective for investigating the dynamic characteristics of a cellular system (Choi et al., 2012; Kim et al., 2014). It is known that complex signaling networks exhibit ordered (stable) dynamics, and that distinct cellular phenotypes can be described by attractor states (Choi et al., 2012; Kim et al., 2014). Attractors are stable states in which a dynamic system (i.e. the cardiac signaling network) tends to converge over time (Choi et al., 2012). The dynamic characteristics of a complex signaling network could be mapped to an attractor landscape (Lee et al., 2015). Every point in the attractor landscape corresponds to a specific network state which is defined by the activity state profile of every molecule in the network (Huang et al., 2009). The state space of the Boolean network with N nodes consists of 2N different states and the dynamics of the state transitions are characterized by the flow of sequential states which eventually converge towards attractors (Lee et al., 2015). The attractor is a particular subset of states, either a fixed point (a point attractor) which remains stable as a single network state, or a limited cycle of period p (a cyclic attractor), consisting of p states that are repeatedly visited by the network dynamics (Huang et al., 2005). The set of initial states with trajectories converging to a given attractor is the basin of attraction (Huang et al., 2005). Based on the state space analysis, distinct attractors of the cardiac signaling network model were classified into three groups characterized by particular cellular phenotypes. The apoptosis or hypertrophy phenotype was defined as an attractor in which any of output nodes used as indicators for corresponding phenotypes is ‘ON’ at least once in its cyclic state transition. The remaining attractors where all the output nodes maintain ‘OFF’ state were determined as a normal phenotype. During the attractor landscape analysis, physiological, ischemic HF, or pressure overload-induced HF condition was reflected to the network model by assigning the state value of particular input nodes to ‘1’ with an intensity determined in its corresponding input setting (Figure 3) while sampling a given number of initial states. From attractor landscape analysis of the cardiac signaling network model, we revealed that the basin sizes of attractors for normal phenotype were highly decreased whereas those for apoptosis or hypertrophy were dramatically increased under HF conditions compared to the case of physiological condition (Supplementary Figure S5 and Table S7). We note that distributions of the relative basin sizes of attractors are maintained regardless of the sampling size of the initial states (Supplementary Figures S5 and S6). In summary, the results suggest that apoptosis or hypertrophy of cardiomyocytes occurs through the coordinated activation of various pathways in the cardiac signaling network under HF conditions. Identification of distinct control nodes that effectively suppress apoptosis and hypertrophy of cardiomyocytes under HF conditions The cellular behavior is controlled by complicated regulatory interactions of various molecules that constitute the intracellular signaling network (Barabasi and Oltvai, 2004). Previous studies showed that intervention of only a small portion of the molecules, which we call control nodes, is needed to govern such complicated network dynamics and to control toward a desired cellular phenotype (Kim et al., 2011, 2013). To identify a set of nodes that can effectively suppress apoptosis and hypertrophy of cardiomyocytes under ischemic or pressure overload-induced HF condition, we employed genetic algorithm (GA) which is widely used to solve complex optimization problems (see Materials and methods for details). As a result, distinct sets of control nodes were found to be important for effectively suppressing apoptosis and hypertrophy of cardiomyocytes: G proteins (q) alpha subunit and G proteins (γ) alpha subunit (Gαq/11), TNF receptor associated factor (TRAF), PI3K, AC, and the growth factor receptor-bound protein 2 (Grb2) and Son of Sevenless (SOS) complex (G/S) under ischemic HF condition and TGFβ receptor (TGFR), Integrins, PI3K, Gαq/11, TRAF, Ras, and G proteins (s) alpha subunit (Gsα) under pressure overload-induced HF condition. These model predictions are supported by previous experimental findings: for instance, inhibition of AC, an identified control node under ischemic HF condition, was shown to suppress apoptosis of cardiomyocytes in response to ischemia injury (Appukuttan et al., 2012; Bravo et al., 2016). Likewise, myocyte knockdown of TGFR, an identified control node under pressure overload-induced HF condition, significantly reduced maladaptive hypertrophy of cardiomyocytes in transverse aortic constriction (an experimental model used for inducing pressure overload-induced cardiac hypertrophy and apoptosis in animals)-operated mice (Koitabashi et al., 2011). We further found that simultaneous inhibition of the control nodes highly decreases the basin sizes of attractors for apoptosis and hypertrophy under HF conditions, implying that the activities of these nodes should be negatively controlled to exhibit therapeutic effects (Figure 4A). To investigate the correlation between the number of controlled nodes and the basin sizes of attractors for normal phenotype, apoptosis, and hypertrophy, we have investigated the cellular behavior while controlling activities of different numbers of randomly sampled nodes. As a result, we found that the maximal basin sizes of attractors for normal phenotype were increased when a more number of sampled nodes were controlled, but controlling a few molecular targets was sufficient to regulate the cellular responses of cardiomyocytes (Supplementary Figure S7). Taken together, these results suggest that only a small fraction of network nodes is required to be controlled to obtain desired network states, which is consistent with the previous findings in other biological networks (Liu et al., 2011; Kim et al., 2013). Figure 4 View largeDownload slide Effect of inhibition of the control nodes on the cellular responses of cardiomyocytes under HF conditions. (A) Distribution of the basin sizes of attractors before and after inhibition of the control nodes under ischemic or pressure overload-induced HF condition. The colored boxes in the right panel describe the classified attractors representing different cellular phenotypes. White or black boxes represent cell phenotype as ‘off’ or ‘on’, respectively. Detailed information about the attractor classification is shown in the main text. (B−D) Distribution of the basin sizes of attractors for normal phenotype (B), apoptosis (C), or hypertrophy (D) after inhibiting the activities of a fixed number of sampled control nodes or random nodes under ischemic or pressure overload-induced HF condition. Simulation analysis was performed while controlling activities of possible combinations of control nodes or randomly sampled nodes and was repeated for 10 times using different random seeds of one million initial states following the same procedure. Data represent mean + SEM of all simulation results. *P < 0.05; **P < 0.01; ***P < 0.001; ****P < 0.0001; Student’s t-test. Figure 4 View largeDownload slide Effect of inhibition of the control nodes on the cellular responses of cardiomyocytes under HF conditions. (A) Distribution of the basin sizes of attractors before and after inhibition of the control nodes under ischemic or pressure overload-induced HF condition. The colored boxes in the right panel describe the classified attractors representing different cellular phenotypes. White or black boxes represent cell phenotype as ‘off’ or ‘on’, respectively. Detailed information about the attractor classification is shown in the main text. (B−D) Distribution of the basin sizes of attractors for normal phenotype (B), apoptosis (C), or hypertrophy (D) after inhibiting the activities of a fixed number of sampled control nodes or random nodes under ischemic or pressure overload-induced HF condition. Simulation analysis was performed while controlling activities of possible combinations of control nodes or randomly sampled nodes and was repeated for 10 times using different random seeds of one million initial states following the same procedure. Data represent mean + SEM of all simulation results. *P < 0.05; **P < 0.01; ***P < 0.001; ****P < 0.0001; Student’s t-test. To further evaluate whether the control nodes could effectively suppress apoptosis and hypertrophy of cardiomyocytes, we have investigated the cellular response while inhibiting all possible combinations of a fixed number of sampled control nodes under HF conditions. As a result, we found that the inhibition of identified control nodes significantly decreased the basin sizes of attractors for apoptosis and hypertrophy, but increased those for normal phenotype compared to the cases when the same number of nodes was randomly chosen (Figure 4B–D). Moreover, the efficacy for regulating the cellular behavior was improved when an increasing number of control nodes were inhibited (Figure 4B–D). These results were still consistent under the given model perturbations (Supplementary Figure S8, see Materials and methods for details). In summary, the result suggests that the identified sets of control nodes are important to effectively suppress apoptosis and hypertrophy of cardiomyocytes under HF conditions. Targeting control nodes increases the efficacy of currently available clinical drugs for HF treatment Previous studies suggested that nodes which are used as drug targets globally influence the cellular behavior (Spiro et al., 2008). Considering that the identified control nodes are found to be important for effectively suppressing apoptosis and hypertrophy of cardiomyocytes under HF conditions, we hypothesized that the control nodes might serve as drug targets. To this end, we have collected all approved or potential drug targets for HF from DrugBank database and literature search (Supplementary Table S8). Next, we compared the ratio of drug targets among the nodes of the control nodes with that of non-control nodes for the cardiac signaling network. Intriguingly, we found that the proportion of drug targets from control nodes was significantly higher than that from non-control nodes (Figure 5A). This suggests that the identified control nodes may be closely related to drug targets. Figure 5 View largeDownload slide Intervention of control nodes increases the efficacy of clinical drugs for HF treatment. (A) The fraction of drug targets in control nodes is significantly higher than that in non-control nodes under both ischemic and pressure overload-induced HF conditions. *P < 0.05; **P < 0.01; Fisher’s exact test. (B−D) Distribution of the increased (B) or decreased (C and D) basin sizes of attractors for normal phenotype (B), apoptosis (C), or hypertrophy (D) after inhibiting a fixed number of sampled control nodes–drug targets when combined with β-blocker treatment under ischemic or pressure overload-induced HF condition. Simulation analysis was performed while controlling activities of possible combinations of control nodes–drug targets or randomly sampled drug targets and was repeated for 10 times using different random seeds of one million initial states following the same procedure. Data represent mean + SEM of all simulation results. *P < 0.05; **P < 0.01; ***P < 0.001; ****P < 0.0001; Student’s t-test. See also Supplementary Table S9. Figure 5 View largeDownload slide Intervention of control nodes increases the efficacy of clinical drugs for HF treatment. (A) The fraction of drug targets in control nodes is significantly higher than that in non-control nodes under both ischemic and pressure overload-induced HF conditions. *P < 0.05; **P < 0.01; Fisher’s exact test. (B−D) Distribution of the increased (B) or decreased (C and D) basin sizes of attractors for normal phenotype (B), apoptosis (C), or hypertrophy (D) after inhibiting a fixed number of sampled control nodes–drug targets when combined with β-blocker treatment under ischemic or pressure overload-induced HF condition. Simulation analysis was performed while controlling activities of possible combinations of control nodes–drug targets or randomly sampled drug targets and was repeated for 10 times using different random seeds of one million initial states following the same procedure. Data represent mean + SEM of all simulation results. *P < 0.05; **P < 0.01; ***P < 0.001; ****P < 0.0001; Student’s t-test. See also Supplementary Table S9. Presently, several drugs (β-, α-, angiotensin II receptor (AT1R)-, and ET1 receptor (ET1R)-blockers) are clinically used for the HF patients, but their efficacy is often limited for blocking the disease progression (Rehsia and Dhalla, 2010). Since control nodes are largely overlapped with the drug targets, we further investigated whether the intervention of the control nodes that are drug targets (hereafter control nodes–drug targets) could improve the therapeutic effects of currently available clinical drugs. For this purpose, we have investigated the cellular response while inhibiting all possible combinations of a fixed number of sampled control nodes–drug targets when combined with each clinically available drug under HF conditions. Note that during the mathematical simulation, the effect of the clinical drugs (i.e. β-, α-, AT1R-, and ET1R-blockers) was reflected by pinning the state of their corresponding target nodes, β1/β2AR, αAR, AT1R, and ET1R, to ‘OFF’, respectively. As a result, inhibition of sampled control nodes–drug targets decreased the basin sizes of attractors for apoptosis and hypertrophy, but increased those for normal phenotype when combined with each clinical drug (Figure 5B–D and Supplementary Table S9). These results were still consistent under the given model perturbations (Supplementary Figures S9–S12, see Materials and methods for details). Intriguingly, optimal combinations of control nodes–drug targets that are most effective in increasing the basin sizes of attractors for normal phenotype were different depending on the types of HF or clinical drugs, suggesting the importance of a system-level development of mechanism-based therapeutic strategies for HF (Table 2). Table 2 The optimal combinations of control nodes for increasing the basin sizes of attractors for normal phenotype under each distinct HF condition. The optimal combinations of control nodes–drug targets that are most effective in increasing the basin sizes of attractors for normal phenotype when combined with each clinical drug under distinct HF condition. Black boxes indicate corresponding optimal combinations of control nodes. Table 2 The optimal combinations of control nodes for increasing the basin sizes of attractors for normal phenotype under each distinct HF condition. The optimal combinations of control nodes–drug targets that are most effective in increasing the basin sizes of attractors for normal phenotype when combined with each clinical drug under distinct HF condition. Black boxes indicate corresponding optimal combinations of control nodes. To further evaluate whether G/S or Gsα, the identified control node that is not included in drug targets, could also improve the therapeutic effects of currently available clinical drugs, we have investigated the cellular behavior while perturbing activities of G/S and Gsα when combined with each clinical drug under ischemic and pressure overload-induced HF conditions, respectively. As a result, we found that the inhibition of G/S or Gsα together with combinations of sampled control nodes–drug targets decreased the basin sizes of attractors for apoptosis and hypertrophy, but increased those for normal phenotype when combined with each clinical drug compared to the cases when the same number of drug targets was randomly chosen (Supplementary Table S10). This result suggests that G/S or Gsα may have a potential to serve as a therapeutic target for HF treatment. Discussion HF is a pathophysiological state in which the heart is unable to pump out a sufficient amount of blood throughout the body (Hilfiker-Klemer et al., 2006). Apoptosis and hypertrophy of cardiomyocytes are the leading causes of HF and they are the typically observed pathological phenotypes in cardiomyocytes during the progression of HF (Remme and Swedberg, 2001). A number of previous studies have used in silico simulation analysis to understand the complicated regulatory mechanisms in cardiomyocytes and to develop the effective therapeutic strategies for HF (Yang and Saucerman, 2011). However, those studies still exhibited limitations because their network models were established focusing only on either apoptosis or hypertrophy separately, and only a few relevant signaling pathways were partially considered in the models (Shin et al., 2011, 2014; Yang and Saucerman, 2011; Ryall et al., 2012; Kang et al., 2017). To overcome such limitations, we have developed a novel mathematical model of the cardiac signaling network, the most comprehensive network among currently existing models to the best of our knowledge, by collecting all available relevant experimental information, and investigated effective therapeutic strategies for HF treatment. Our novel findings obtained from these efforts led to the following conclusions: (i) we have identified distinct sets of control nodes that can effectively suppress apoptosis and hypertrophy of cardiomyocytes under ischemic or pressure overload-induced HF condition. (ii) Intervention of the control nodes may increase the efficacy of clinical drugs for HF. (iii) Different combinations of control nodes are suggested to serve as potentially effective candidate drug targets depending on the types of HF. In general, HF is not a single disease entity but can be considered as a wide disease spectrum including various pathological phenotypes (Ponikowski et al., 2016). In this regard, the underlying molecular pathophysiology can be different between HF patients although they have the same clinical features (De Keulenaer and Brutsaert, 2007). In other words, individual HF patient follows a unique disease trajectory depending on his/her clinical presentations, comorbidities, as well as the molecular pathophysiology including complex interplay of signaling dynamics at the molecular level (De Keulenaer and Brutsaert, 2007). However, the current guideline for HF treatment is largely based on the clinical features of the disease such as symptoms and signs rather than the molecular level understanding (Yancy et al., 2013; Ponikowski et al., 2016), and, therefore, the efficacy of currently available drugs considerably varies between HF patients and the drug treatment may even cause unwanted adverse effects (Yancy et al., 2013). Here, intriguingly, distinct combinations of control nodes were suggested as potentially effective candidate drug targets depending on HF types. This finding implies that individualized drug treatment strategies based on the molecular mechanism of HF can be used to achieve optimal therapeutic responses, underscoring the urgent need to develop mechanism-based therapeutic strategies for HF. We also note that the drugs that are widely used for HF patients in the clinics are mostly receptor antagonists (Law et al., 2014). Considering that the identified control nodes comprise not only receptors but also intracellular kinases, the results suggest a rationale for the development of the downstream kinase inhibitors for HF treatment (Law et al., 2014). Collectively, our findings may contribute to the improvement of the current guideline for HF treatment and the development of effective therapeutic strategies. Epigenetic landscape approach has been widely employed to study the complex regulatory dynamics of various biological phenomena (Wang et al., 2010; Li and Wang, 2013, 2014a, b; Li, 2017). This method attempts to uncover the functional mechanism of transition between distinct cellular states by analyzing the landscape topography and comparing the height of barriers that separate different basins of attraction in biological network systems (Wang et al., 2010; Li and Wang, 2013, 2014a, b; Li, 2017). Previous studies have proven the usefulness of landscape topography analysis in investigating biological processes and cellular functions since it provides a quantitative way of understanding the dynamic behavior of complex intracellular signaling networks (Wang et al., 2010; Li and Wang, 2013, 2014a, b; Li, 2017). Based on the epigenetic landscape concept and topography analysis of continuous mathematical models of specific gene regulatory systems, including stem cell differentiation, cancer, and cancer immunity (Li and Wang, 2013, 2014b; Li, 2017), they revealed key factors that are critically responsible for determining the barrier height between different attractors and identified optimal intervention strategies that may produce desired outcomes by controlling the barrier height (Li and Wang, 2013, 2014b; Li, 2017). In this study, we have developed a large-scale discrete Boolean model of cardiac signaling network and focused on the investigation of distinct sets of control nodes that can effectively maximize the basin sizes of attractors for normal phenotype under HF conditions, rather than analyzing the landscape topography of the network model. A future study employing the epigenetic landscape approach based on topography analysis of our cardiac signaling network model will further advance the understanding of regulatory dynamics of HF and may provide a rationale for a novel therapeutic strategy for HF treatment. In this study, we have developed a mathematical model of a comprehensive cardiac signaling network by integrating available cardiomyocyte-specific experimental evidence through an extensive search over various literatures and databases. However, there were some difficulties in establishing a large-scale network model by only using mechanical information obtained from such literatures and databases due to the intrinsic biological noise and false positive results of biochemical experiments, and the lack of sufficient available information about the cooperativity among the incoming regulatory signals of biomolecules (Macarthur et al., 2009). To overcome such challenge, we used GA, an optimization algorithm (Goldberg, 1989), in assigning the most appropriate logic tables to the nodes that have multiple positive and/or negative upstream regulators as such that the mathematical model can well reproduce the cellular responses of cardiomyocytes (see Materials and methods). Sufficient accumulation of biological data that can provide accurate mechanistic information on functional regulation of genes and proteins based on advanced experimental techniques would facilitate the development of more reliable mathematical model. Moreover, we mainly focused on apoptosis and hypertrophy of cardiomyocytes as key molecular phenotypic features of HF and investigated optimal therapeutic strategies that can effectively suppress such maladaptive cellular responses. However, myocardial fibrosis developed by complex heterocellular interactions between cardiomyocytes and other various cell types, including cardiac fibroblasts, is also known to play an important role in the HF progression (Kong et al., 2014; Leask, 2015; Gourdie et al., 2016; Travers et al., 2016). Development and analysis of mathematical models that can comprehensively encompass both intra- and intercellular signaling mechanisms of pathophysiology of HF may provide more improved therapeutic approaches. In summary, from attractor landscape analysis of the cardiac signaling network model, we have identified distinct sets of control nodes that can effectively suppress apoptosis and hypertrophy of cardiomyocytes under ischemic or pressure overload-induced HF condition. We further found that intervention of these control nodes may increase the efficacy of clinical drugs for HF treatment and that different combinations of control nodes are suggested as potentially effective candidate drug targets depending on the types of HF. Our study provides a systematic way of developing effective mechanism-based therapeutic strategies for HF. Materials and methods Development of the cardiac signaling network model The cardiac signaling network was reconstructed by incorporating all relevant information on the complicated interactions between individual proteins in cardiomyocytes through an extensive search over literature and database, and it comprises 141 nodes and 274 directed links. The 274 directed links comprising the cardiac signaling network are shown in Supplementary Dataset S1. The column ‘Source’ indicates the molecules that regulate downstream ‘Target’ molecules; and the column ‘Target’ indicates the molecules being activated or inhibited by the upstream ‘Source’ molecules. In the column ‘Interaction’, ‘+’ and ‘−’ represent activation and inhibition, respectively. The activity states of intracellular signaling nodes except for the input nodes were determined by their assigned logic tables as shown in Supplementary Dataset S2. The logic table for each node was determined based on the biological information on activation and/or inhibition of distinct molecules obtained from the related literature and database search (Supplementary Dataset S2). For the cases of nodes that have two or more positive and/or negative upstream regulators, however, there were some difficulties in generating their logic tables based only on the literature mining due to the lack of sufficient available information about the cooperativity among the incoming regulatory signals. Therefore, we used the weighted sum Boolean function to assign the most appropriate logic tables to the nodes with at least two positive and/or negative input nodes as in the previous study (Fumia and Martins, 2013). In this weighted sum modeling method, each node i has only two states, Si = 1 and Si = 0, representing the active and the inactive state of a node, respectively. The state values of nodes in the next time step (t + 1) are determined by the state values of their input nodes in the present time step (t) according to the following rule: Si(t+1)=sgn(∑j=1kin(i)JjiSj(t)−θi) Here, Jji is the interaction weight from input j on node i, where activating and inhibitory interactions are positive and negative numbers, respectively. kin(i) indicates the number of input nodes that regulate node i and determine its activation state. The threshold function sgn(x) is a unitary step function (sgn(x) = 1 if x > 0, sgn(x) = 0 if x ≤ 0). θi is the activation threshold of node i. Therefore, if a weighted sum of input signals of node i exceeds its activation threshold θi, the node is activated or stays active if it was already active; otherwise, it becomes inactive or stays inactive. We employed GA (Goldberg, 1989), which is widely used to solve complex optimization problems, to identify the optimal values of interaction weight Jji and activation threshold θi that can reasonably recapitulate the dynamic activity profiles of network components as reported from previous experiments (Figure 2). All the Boolean simulation analyses of the cardiac signaling network model were implemented using Matlab® 2014b in a Window Cluster system composed of 288 CPUs in parallel. In silico simulation analysis of input–output relationships of network components While the state space of the cardiac signaling network is huge (2141 states), we found that each network state eventually converges to a relatively small periodic cycle within 300 time steps from numerous iterative simulations. Thus, we generated Boolean function updates for 1000 time steps to measure the stable state activities of nodes using various levels of input stimulations. To represent the input intensity of each simulation run, we placed the input node in the cycle that yields the desired percentage of ON (Lee et al., 2015). For instance, if an input is set to 50% ON, it would be put on a cycle of ‘1010’ or ‘0101’, and this condition is fixed during the mathematical simulation. A small percentage of background noise (5%, described below) was added to the inputs to produce biologically reliable characteristics. For instance, if an input is set to 40% ON for a given run under 5% noise, the input changes randomly from 35% to 45% ON. The initial condition of each node (except input nodes) was randomized in every run (50% ‘0’, 50% ‘1’). The average ‘1’s over the last 50, 100, or 200 time steps were calculated and used to obtain the steady-state activity of each node. The results showed that the input–output relationships of the candidate nodes are similar regardless of the number of last average steps (Figure 2; Supplementary Figures S2 and S3). Note that, during the extensive simulation for the input–output relationship analysis of the cardiac signaling network model, the input intensity of a particular input node of interest was set to 0%−100% ON with 1% increments while intensities of remaining input nodes were fixed to the optimal input settings, as previously described (Helikar et al., 2008) (Supplementary Table S4). Attractor landscape analysis of the cardiac signaling network model To understand the cellular characteristics of the cardiac signaling network model, we performed the attractor landscape analysis (Choi et al., 2012). The attractor landscape is composed of attractors, the stable states of the network (Cho et al., 2016). The set of network states that converge to the corresponding attractors is called the basin of attraction (Huang et al., 2005). For a given input condition, we obtained the attractor landscape of the cardiac signaling network model using one million randomly sampled initial states until the state trajectory of an initial state reaches a point attractor with a single network state or a cyclic attractor with repetitive states (Choi et al., 2012; Cho et al., 2016). Note that the simulation result was maintained regardless of the sampling size of initial states (Supplementary Figures S5 and S6). During the attractor landscape analysis, physiological, ischemic HF, or pressure overload-induced HF condition was reflected to the network model by assigning the state value of particular input nodes to ‘1’ with an intensity determined in its corresponding input setting (Figure 3) while sampling a given number of initial states. To examine the robustness of the simulation results, we have performed additional computer simulation experiments while randomly perturbing the network structure and parameter settings to confirm the robustness of the dynamic analysis. Of note, two forms of network model perturbations were given in a combined manner during the attractor landscape analysis: (i) we have randomly rewired indicated percentage of the total number of links comprising the cardiac signaling network and (ii) logic tables of the indicated percentage of all internal regulatory nodes in the network were randomly perturbed. Identification of control nodes To identify the control nodes for minimizing the basin sizes of attractors for apoptosis and hypertrophy under ischemic or pressure overload-induced HF condition, we adopted the previously developed algorithm which identifies a set of nodes needed to be regulated to control the network states to converge to a desired attractor (Cho et al., 2016). During the mathematical simulation, the state value of each candidate control node was pinned to either 0 or 1. We employed GA that is used for optimization and search problems by artificially mutating, evolving, and selecting the most promising chromosomes containing all possible solutions (Kim et al., 2013). The fitness function was adopted from Cho et al. (2016), and is as follows: Fitness={B3×(n−∑i=1nXi)2×2,ifB=1B3×(n−∑i=1nXi)2,otherwise where B is the basin sizes of attractors for a desired phenotype, n is the number of nodes in the cardiac signaling network, and Xi is either 1 or 0 when node i is included or not in the chromosome X, respectively (Cho et al., 2016). Specifically, the basin size B is defined as follows: the basin of attraction converging to the attractors for normal phenotype when the candidate control nodes were perturbed under ischemic or pressure overload-induced HF condition. Here, attractors where all the output nodes used as indicators of hypertrophy or apoptosis of cardiomyocytes maintain ‘OFF’ state were determined as a normal phenotype. To measure the basin size B, we performed the attractor analysis of the cardiac signaling network model using one million randomly sampled initial states until the state trajectory of an initial state reaches a point attractor with a single network state or a cyclic attractor with repetitive states. Then, the fraction of initial states that converge to the attractors for normal phenotype was determined as the basin size B. The node set with the highest fitness value (i.e. the node set that can maximize the basin size B under distinct HF conditions) in the chromosome is chosen as a set of control nodes (Cho et al., 2016). Supplementary material Supplementary material is available at Journal of Molecular Cell Biology online. Funding This work was supported by the National Research Foundation of Korea (NRF) grants funded by the Korea Government, the Ministry of Science and ICT (2017R1A2A1A17069642, 2015M3A9A7067220, and 2013M3A9A7046303). Conflict of interest none declared. References ACE Inhibitor Myocardial Infarction Collaborative Group . ( 1998 ). Indications for ACE inhibitors in the early treatment of acute myocardial infarction: systematic overview of individual data from 100000 patients in randomized trials . Circulation 97 , 2202 – 2212 . CrossRef Search ADS PubMed Aikawa , R. , Komuro , I. , Yamazaki , T. , et al. . ( 1997 ). Oxidative stress activates extracellular signal-regulated kinases through Src and Ras in cultured cardiac myocytes of neonatal rats . J. Clin. Invest. 100 , 1813 – 1821 . Google Scholar CrossRef Search ADS PubMed Appukuttan , A. , Kasseckert , S.A. , Micoogullari , M. , et al. . ( 2012 ). Type 10 adenylyl cyclase mediates mitochondrial Bax translocation and apoptosis of adult rat cardiomyocytes under simulated ischaemia/reperfusion . Cardiovasc. Res. 93 , 340 – 349 . Google Scholar CrossRef Search ADS PubMed Backs , J. , Worst , B.C. , Lehmann , L.H. , et al. . ( 2011 ). Selective repression of MEF2 activity by PKA-dependent proteolysis of HDAC4 . J. Cell Biol. 195 , 403 – 415 . Google Scholar CrossRef Search ADS PubMed Barabasi , A.L. , Gulbahce , N. , and Loscalzo , J. ( 2011 ). Network medicine: a network-based approach to human disease . Nat. Rev. Genet. 12 , 56 – 68 . Google Scholar CrossRef Search ADS PubMed Barabasi , A.L. , and Oltvai , Z.N. ( 2004 ). Network biology: understanding the cell’s functional organization . Nat. Rev. Genet. 5 , 101 – 113 . Google Scholar CrossRef Search ADS PubMed Bernardo , B.C. , Weeks , K.L. , Pretorius , L. , et al. . ( 2010 ). Molecular distinction between physiological and pathological cardiac hypertrophy: experimental findings and therapeutic strategies . Pharmacol. Ther. 128 , 191 – 227 . Google Scholar CrossRef Search ADS PubMed Blankesteijn , W.M. , van de Schans , V.A. , ter Horst , P. , et al. . ( 2008 ). The Wnt/frizzled/GSK-3β pathway: a novel therapeutic target for cardiac hypertrophy . Trends Pharmacol. Sci. 29 , 175 – 180 . Google Scholar CrossRef Search ADS PubMed Bravo , C.A. , Vatner , D.E. , Pachon , R. , et al. . ( 2016 ). A food and drug administration-approved antiviral agent that inhibits adenylyl cyclase type 5 protects the ischemic heart even when administered after reperfusion . J. Pharmacol. Exp. Ther. 357 , 331 – 336 . Google Scholar CrossRef Search ADS PubMed Bujak , M. , and Frangogiannis , N.G. ( 2007 ). The role of TGF-β signaling in myocardial infarction and cardiac remodeling . Cardiovasc. Res. 74 , 184 – 195 . Google Scholar CrossRef Search ADS PubMed Chen , Y. , Ke , Q. , Yang , Y. , et al. . ( 2003 ). Cardiomyocytes overexpressing TNF-α attract migration of embryonic stem cells via activation of p38 and c-Jun amino-terminal kinase . FASEB J. 17 , 2231 – 2239 . Google Scholar CrossRef Search ADS PubMed Cho , S.H. , Park , S.M. , Lee , H.S. , et al. . ( 2016 ). Attractor landscape analysis of colorectal tumorigenesis and its reversion . BMC Syst. Biol. 10 , 96 . Google Scholar CrossRef Search ADS PubMed Choi , M. , Shi , J. , Jung , S.H. , et al. . ( 2012 ). Attractor landscape analysis reveals feedback loops in the p53 network that control the cellular response to DNA damage . Sci. Signal. 5 , ra83 . Google Scholar CrossRef Search ADS PubMed Colella , M. , Grisan , F. , Robert , V. , et al. . ( 2008 ). Ca2+ oscillation frequency decoding in cardiac cell hypertrophy: role of calcineurin/NFAT as Ca2+ signal integrators . Proc. Natl Acad. Sci. USA 105 , 2859 – 2864 . Google Scholar CrossRef Search ADS Condorelli , G. , Morisco , C. , Latronico , M.V. , et al. . ( 2002 ). TNF-α signal transduction in rat neonatal cardiac myocytes: definition of pathways generating from the TNF-α receptor . FASEB J. 16 , 1732 – 1737 . Google Scholar CrossRef Search ADS PubMed Cuenca , J. , Martin-Sanz , P. , Alvarez-Barrientos , A.M. , et al. . ( 2006 ). Infiltration of inflammatory cells plays an important role in matrix metalloproteinase expression and activation in the heart during sepsis . Am. J. Pathol. 169 , 1567 – 1576 . Google Scholar CrossRef Search ADS PubMed De Keulenaer , G.W. , and Brutsaert , D.L. ( 2007 ). Systolic and diastolic heart failure: different phenotypes of the same disease? Eur. J. Heart Fail . 9 , 136 – 143 . Google Scholar CrossRef Search ADS PubMed Fischer , P. , and Hilfiker-Kleiner , D. ( 2008 ). Role of gp130-mediated signalling pathways in the heart and its impact on potential therapeutic aspects . Br. J. Pharmacol. 153(Suppl 1) , S414 – S427 . Google Scholar PubMed Florea , V.G. , and Cohn , J.N. ( 2014 ). The autonomic nervous system and heart failure . Circ. Res. 114 , 1815 – 1826 . Google Scholar CrossRef Search ADS PubMed Francis , S.H. ( 2010 ). The role of cGMP-dependent protein kinase in controlling cardiomyocyte cGMP . Circ. Res. 107 , 1164 – 1166 . Google Scholar CrossRef Search ADS PubMed Fu , Y.C. , Chi , C.S. , Yin , S.C. , et al. . ( 2004 ). Norepinephrine induces apoptosis in neonatal rat cardiomyocytes through a reactive oxygen species-TNFα-caspase signaling pathway . Cardiovasc. Res. 62 , 558 – 567 . Google Scholar CrossRef Search ADS PubMed Fumia , H.F. , and Martins , M.L. ( 2013 ). Boolean network model for cancer pathways: predicting carcinogenesis and targeted therapy outcomes . PLoS One 8 , e69008 . Google Scholar CrossRef Search ADS PubMed Goldberg , D.E. ( 1989 ). Genetic Algorithms in Search, Optimization, and Machine Learning (1st edn) . Reading, MA : Addison-Wesley Professional . Gonzalez , A. , Fortuno , M.A. , Querejeta , R. , et al. . ( 2003 ). Cardiomyocyte apoptosis in hypertensive cardiomyopathy . Cardiovasc. Res. 59 , 549 – 562 . Google Scholar CrossRef Search ADS PubMed Gourdie , R.G. , Dimmeler , S. , and Kohl , P. ( 2016 ). Novel therapeutic strategies targeting fibroblasts and fibrosis in heart disease . Nat. Rev. Drug Discov. 15 , 620 – 638 . Google Scholar CrossRef Search ADS PubMed Guo , X. , and Wang , X.F. ( 2009 ). Signaling cross-talk between TGF-β/BMP and other pathways . Cell Res. 19 , 71 – 88 . Google Scholar CrossRef Search ADS PubMed Handler , D. , Meixner , K. , Pizka , M. , et al. . ( 2013 ). The genetic makeup of the Drosophila piRNA pathway . Mol. Cell 50 , 762 – 777 . Google Scholar CrossRef Search ADS PubMed Heineke , J. , and Molkentin , J.D. ( 2006 ). Regulation of cardiac hypertrophy by intracellular signalling pathways . Nat. Rev. Mol. Cell Biol. 7 , 589 – 600 . Google Scholar CrossRef Search ADS PubMed Helikar , T. , Konvalina , J. , Heidel , J. , et al. . ( 2008 ). Emergent decision-making in biological signal transduction networks . Proc. Natl Acad. Sci. USA 105 , 1913 – 1918 . Google Scholar CrossRef Search ADS Hilfiker-Klemer , D. , Landmesser , U. , and Drexler , H. ( 2006 ). Molecular mechanisms in heart failure—focus on cardiac hypertrophy, inflammation, angiogenesis, and apoptosis . J. Am. Coll. Cardiol. 48 , A56 – A66 . Google Scholar CrossRef Search ADS Huang , S. , Eichler , G. , Bar-Yam , Y. , et al. . ( 2005 ). Cell fates as high-dimensional attractor states of a complex gene regulatory network . Phys. Rev. Lett. 94 , 128701 . Google Scholar CrossRef Search ADS PubMed Huang , S. , Ernberg , I. , and Kauffman , S. ( 2009 ). Cancer attractors: a systems view of tumors from a gene network dynamics and developmental perspective . Semin. Cell Dev. Biol. 20 , 869 – 876 . Google Scholar CrossRef Search ADS PubMed Iwamoto , R. , Yamazaki , S. , Asakura , M. , et al. . ( 2003 ). Heparin-binding EGF-like growth factor and ErbB signaling is essential for heart function . Proc. Natl Acad. Sci. USA 100 , 3221 – 3226 . Google Scholar CrossRef Search ADS Jackson , G. , Gibbs , C.R. , Davies , M.K. , et al. . ( 2000 ). ABC of heart failure: Pathophysiology . BMJ 320 , 167 – 170 . Google Scholar CrossRef Search ADS PubMed Kanehisa , M. , and Goto , S. ( 2000 ). KEGG: kyoto encyclopedia of genes and genomes . Nucleic Acids Res. 28 , 27 – 30 . Google Scholar CrossRef Search ADS PubMed Kang , J.H. , Lee , H.S. , Kang , Y.W. , et al. . ( 2016 ). Systems biological approaches to the cardiac signaling network . Brief. Bioinform. 17 , 419 – 428 . Google Scholar CrossRef Search ADS PubMed Kang , J.H. , Lee , H.S. , Park , D. , et al. . ( 2017 ). Context-independent essential regulatory interactions for apoptosis and hypertrophy in the cardiac signaling network . Sci. Rep. 7 , 34 . Google Scholar CrossRef Search ADS PubMed Kawasaki , K. , Watabe , T. , Sase , H. , et al. . ( 2008 ). Ras signaling directs endothelial specification of VEGFR2+ vascular progenitor cells . J. Cell Biol. 181 , 131 – 141 . Google Scholar CrossRef Search ADS PubMed Kehat , I. , and Molkentin , J.D. ( 2010 ). Molecular pathways underlying cardiac remodeling during pathophysiological stimulation . Circulation 122 , 2727 – 2735 . Google Scholar CrossRef Search ADS PubMed Kim , S. , Kim , J. , and Cho , K.H. ( 2007 ). Inferring gene regulatory networks from temporal expression profiles under time-delay and noise . Comput. Biol. Chem. 31 , 239 – 245 . Google Scholar CrossRef Search ADS PubMed Kim , J.R. , Kim , J. , Kwon , Y.K. , et al. . ( 2011 ). Reduction of complex signaling networks to a representative kernel . Sci. Signal. 4 , ra35 . Google Scholar CrossRef Search ADS PubMed Kim , J. , Park , S.M. , and Cho , K.H. ( 2013 ). Discovery of a kernel for controlling biomolecular regulatory networks . Sci. Rep. 3 , 2223 . Google Scholar CrossRef Search ADS PubMed Kim , J. , Vandamme , D. , Kim , J.R. , et al. . ( 2014 ). Robustness and evolvability of the human signaling network . PLoS Comput. Biol. 10 , e1003763 . Google Scholar CrossRef Search ADS PubMed Kitano , H. ( 2002 ). Computational systems biology . Nature 420 , 206 – 210 . Google Scholar CrossRef Search ADS PubMed Koitabashi , N. , Danner , T. , Zaiman , A.L. , et al. . ( 2011 ). Pivotal role of cardiomyocyte TGF-β signaling in the murine pathological response to sustained pressure overload . J. Clin. Invest. 121 , 2301 – 2312 . Google Scholar CrossRef Search ADS PubMed Kong , P. , Christia , P. , and Frangogiannis , N.G. ( 2014 ). The pathogenesis of cardiac fibrosis . Cell. Mol. Life Sci. 71 , 549 – 574 . Google Scholar CrossRef Search ADS PubMed Kuramochi , Y. , Guo , X. , and Sawyer , D.B. ( 2006 ). Neuregulin activates erbB2-dependent src/FAK signaling and cytoskeletal remodeling in isolated adult rat cardiac myocytes . J. Mol. Cell. Cardiol. 41 , 228 – 235 . Google Scholar CrossRef Search ADS PubMed Law , V. , Knox , C. , Djoumbou , Y. , et al. . ( 2014 ). DrugBank 4.0: shedding new light on drug metabolism . Nucleic Acids Res. 42 , D1091 – D1097 . Google Scholar CrossRef Search ADS PubMed Leask , A. ( 2015 ). Getting to the heart of the matter: new insights into cardiac fibrosis . Circ. Res. 116 , 1269 – 1276 . Google Scholar CrossRef Search ADS PubMed Lee , H.S. , Goh , M.J. , Kim , J. , et al. . ( 2015 ). A systems-biological study on the identification of safe and effective molecular targets for the reduction of ultraviolet B-induced skin pigmentation . Sci. Rep. 5 , 10305 . Google Scholar CrossRef Search ADS PubMed Li , C. ( 2017 ). Identifying the optimal anticancer targets from the landscape of a cancer-immunity interaction network . Phys. Chem. Chem. Phys. 19 , 7642 – 7651 . Google Scholar CrossRef Search ADS PubMed Li , C. , and Wang , J. ( 2013 ). Quantifying cell fate decisions for differentiation and reprogramming of a human stem cell network: landscape and biological paths . PLoS Comput. Biol. 9 , e1003165 . Google Scholar CrossRef Search ADS PubMed Li , C. , and Wang , J. ( 2014 a). Landscape and flux reveal a new global view and physical quantification of mammalian cell cycle . Proc. Natl Acad. Sci. USA 111 , 14130 – 14135 . Google Scholar CrossRef Search ADS Li , C. , and Wang , J. ( 2014 b). Quantifying the underlying landscape and paths of cancer . J. R. Soc. Interface 11 , 20140774 . Google Scholar CrossRef Search ADS PubMed Lim , J.Y. , Park , S.J. , Hwang , H.Y. , et al. . ( 2005 ). TGF-β1 induces cardiac hypertrophic responses via PKC-dependent ATF-2 activation . J. Mol. Cell. Cardiol. 39 , 627 – 636 . Google Scholar CrossRef Search ADS PubMed Ling , H. , Gray , C.B. , Zambon , A.C. , et al. . ( 2013 ). Ca2+/Calmodulin-dependent protein kinase II delta mediates myocardial ischemia/reperfusion injury through nuclear factor-κB . Circ. Res. 112 , 935 – 944 . Google Scholar CrossRef Search ADS PubMed Liu , S. , Li , Y. , Kim , S. , et al. . ( 2012 ). Phosphodiesterases coordinate cAMP propagation induced by two stimulatory G protein-coupled receptors in hearts . Proc. Natl Acad. Sci. USA 109 , 6578 – 6583 . Google Scholar CrossRef Search ADS Liu , Y.Y. , Slotine , J.J. , and Barabasi , A.L. ( 2011 ). Controllability of complex networks . Nature 473 , 167 – 173 . Google Scholar CrossRef Search ADS PubMed Loan Le , T.Y. , Mardini , M. , Howell , V.M. , et al. . ( 2012 ). Low-dose spironolactone prevents apoptosis repressor with caspase recruitment domain degradation during myocardial infarction . Hypertension 59 , 1164 – 1169 . Google Scholar CrossRef Search ADS PubMed Lu , J. , Zeng , H. , Liang , Z. , et al. . ( 2015 ). Network modelling reveals the mechanism underlying colitis-associated colon cancer and identifies novel combinatorial anti-cancer targets . Sci. Rep. 5 , 14739 . Google Scholar CrossRef Search ADS PubMed Macarthur , B.D. , Ma’ayan , A. , and Lemischka , I.R. ( 2009 ). Systems biology of stem cell fate and cellular reprogramming . Nat. Rev. Mol. Cell Biol. 10 , 672 – 681 . Google Scholar CrossRef Search ADS PubMed MacLellan , W.R. , and Schneider , M.D. ( 2000 ). Genetic dissection of cardiac growth control pathways . Annu. Rev. Physiol. 62 , 289 – 319 . Google Scholar CrossRef Search ADS PubMed Matsui , T. , and Rosenzweig , A. ( 2005 ). Convergent signal transduction pathways controlling cardiomyocyte survival and function: the role of PI 3-kinase and Akt . J. Mol. Cell. Cardiol. 38 , 63 – 71 . Google Scholar CrossRef Search ADS PubMed Murray , P.J. , Kang , J.W. , Mirams , G.R. , et al. . ( 2010 ). Modelling spatially regulated β-catenin dynamics and invasion in intestinal crypts . Biophys. J. 99 , 716 – 725 . Google Scholar CrossRef Search ADS PubMed Nadal-Ginard , B. , Kajstura , J. , Leri , A. , et al. . ( 2003 ). Myocyte death, growth, and regeneration in cardiac hypertrophy and failure . Circ. Res. 92 , 139 – 150 . Google Scholar CrossRef Search ADS PubMed Nakamura , T. , Mizuno , S. , Matsumoto , K. , et al. . ( 2000 ). Myocardial protection from ischemia/reperfusion injury by endogenous and exogenous HGF . J. Clin. Invest. 106 , 1511 – 1519 . Google Scholar CrossRef Search ADS PubMed Nakaoka , Y. , Nishida , K. , Fujio , Y. , et al. . ( 2003 ). Activation of gp130 transduces hypertrophic signal through interaction of scaffolding/docking protein Gab1 with tyrosine phosphatase SHP2 in cardiomyocytes . Circ. Res. 93 , 221 – 229 . Google Scholar CrossRef Search ADS PubMed Nakaoka , Y. , Nishida , K. , Narimatsu , M. , et al. . ( 2007 ). Gab family proteins are essential for postnatal maintenance of cardiac function via neuregulin-1/ErbB signaling . J. Clin. Invest. 117 , 1771 – 1781 . Google Scholar CrossRef Search ADS PubMed Narula , J. , Haider , N. , Virmani , R. , et al. . ( 1996 ). Apoptosis in myocytes in end-stage heart failure . N. Engl. J. Med. 335 , 1182 – 1189 . Google Scholar CrossRef Search ADS PubMed Pan , J. , Fukuda , K. , Saito , M. , et al. . ( 1999 ). Mechanical stretch activates the JAK/STAT pathway in rat cardiomyocytes . Circ. Res. 84 , 1127 – 1136 . Google Scholar CrossRef Search ADS PubMed Park , S.G. , Lee , T. , Kang , H.Y. , et al. . ( 2006 ). The influence of the signal dynamics of activated form of IKK on NF-κB and anti-apoptotic gene expressions: a systems biology approach . FEBS Lett. 580 , 822 – 830 . Google Scholar CrossRef Search ADS PubMed Pennanen , C. , Parra , V. , Lopez-Crisosto , C. , et al. . ( 2014 ). Mitochondrial fission is required for cardiomyocyte hypertrophy mediated by a Ca2+-calcineurin signaling pathway . J. Cell Sci. 127 , 2659 – 2671 . Google Scholar CrossRef Search ADS PubMed Phanstiel , D.H. , Brumbaugh , J. , Wenger , C.D. , et al. . ( 2011 ). Proteomic and phosphoproteomic comparison of human ES and iPS cells . Nat. Methods 8 , 821 – 827 . Google Scholar CrossRef Search ADS PubMed Ponikowski , P. , Voors , A.A. , Anker , S.D. , et al. . ( 2016 ). 2016 ESC Guidelines for the diagnosis and treatment of acute and chronic heart failure: the Task Force for the diagnosis and treatment of acute and chronic heart failure of the European Society of Cardiology (ESC) developed with the special contribution of the Heart Failure Association (HFA) of the ESC . Eur. Heart J. 37 , 2129 – 2200 . Google Scholar CrossRef Search ADS PubMed Ponnusamy , M. , Li , P.F. , and Wang , K. ( 2017 ). Understanding cardiomyocyte proliferation: an insight into cell cycle activity . Cell. Mol. Life Sci. 74 , 1019 – 1034 . Google Scholar CrossRef Search ADS PubMed Rehsia , N.S. , and Dhalla , N.S. ( 2010 ). Mechanisms of the beneficial effects of β-adrenoceptor antagonists in congestive heart failure . Exp. Clin. Cardiol. 15 , e86 – e95 . Google Scholar PubMed Reimand , J. , Arak , T. , Adler , P. , et al. . ( 2016 ). g:Profiler—a web server for functional interpretation of gene lists (2016 update) . Nucleic Acids Res. 44 , W83 – W89 . Google Scholar CrossRef Search ADS PubMed Remme , W.J. , and Swedberg , K. ( 2001 ). Guidelines for the diagnosis and treatment of chronic heart failure . Eur. Heart J. 22 , 1527 – 1560 . Google Scholar CrossRef Search ADS PubMed Roberge , S. , Roussel , J. , Andersson , D.C. , et al. . ( 2014 ). TNF-α-mediated caspase-8 activation induces ROS production and TRPM2 activation in adult ventricular myocytes . Cardiovasc. Res. 103 , 90 – 99 . Google Scholar CrossRef Search ADS PubMed Ronca , R. , Gualandi , L. , Crescini , E. , et al. . ( 2009 ). Fibroblast growth factor receptor-1 phosphorylation requirement for cardiomyocyte differentiation in murine embryonic stem cells . J. Cell. Mol. Med. 13 , 1489 – 1498 . Google Scholar CrossRef Search ADS PubMed Ryall , K.A. , Holland , D.O. , Delaney , K.A. , et al. . ( 2012 ). Network reconstruction and systems analysis of cardiac myocyte hypertrophy signaling . J. Biol. Chem. 287 , 42259 – 42268 . Google Scholar CrossRef Search ADS PubMed Saraste , A. , Pulkki , K. , Kallajoki , M. , et al. . ( 1997 ). Apoptosis in human acute myocardial infarction . Circulation 95 , 320 – 323 . Google Scholar CrossRef Search ADS PubMed Schaefer , C.F. , Anthony , K. , Krupa , S. , et al. . ( 2009 ). PID: the pathway interaction database . Nucleic Acids Res. 37 , D674 – D679 . Google Scholar CrossRef Search ADS PubMed Schaub , M.C. , and Hefti , M.A. ( 2007 ). The PGE2-Stat3 connection in cardiac hypertrophy . Cardiovasc. Res. 73 , 3 – 5 . Google Scholar CrossRef Search ADS PubMed Shin , S.Y. , Kim , T. , Lee , H.S. , et al. . ( 2014 ). The switching role of β-adrenergic receptor signalling in cell survival or death decision of cardiomyocytes . Nat. Commun. 5 , 5777 . Google Scholar CrossRef Search ADS PubMed Shin , S.Y. , Yang , H.W. , Kim , J.R. , et al. . ( 2011 ). A hidden incoherent switch regulates RCAN1 in the calcineurin-NFAT signaling network . J. Cell Sci. 124 , 82 – 90 . Google Scholar CrossRef Search ADS PubMed Song , H.K. , Hong , S.E. , Kim , T. , et al. . ( 2012 ). Deep RNA sequencing reveals novel cardiac transcriptomic signatures for physiological and pathological hypertrophy . PLoS One 7 , e35552 . Google Scholar CrossRef Search ADS PubMed Spiro , Z. , Kovacs , I.A. , and Csermely , P. ( 2008 ). Drug-therapy networks and the prediction of novel drug targets . J. Biol. 7 , 20 . Google Scholar CrossRef Search ADS PubMed Stewart , S. , MacIntyre , K. , Capewell , S. , et al. . ( 2003 ). Heart failure and the aging population: an increasing burden in the 21st century? Heart 89 , 49 – 53 . Google Scholar CrossRef Search ADS PubMed Thygesen , K. , Alpert , J.S. , Jaffe , A.S. , et al. . ( 2012 ). Third universal definition of myocardial infarction . Circulation 126 , 2020 – 2035 . Google Scholar CrossRef Search ADS PubMed Travers , J.G. , Kamal , F.A. , Robbins , J. , et al. . ( 2016 ). Cardiac fibrosis: the fibroblast awakens . Circ. Res. 118 , 1021 – 1040 . Google Scholar CrossRef Search ADS PubMed Ucar , A. , Gupta , S.K. , Fiedler , J. , et al. . ( 2012 ). The miRNA-212/132 family regulates both cardiac hypertrophy and cardiomyocyte autophagy . Nat. Commun. 3 , 1078 . Google Scholar CrossRef Search ADS PubMed Vega , R.B. , Harrison , B.C. , Meadows , E. , et al. . ( 2004 ). Protein kinases C and D mediate agonist-dependent cardiac hypertrophy through nuclear export of histone deacetylase 5 . Mol. Cell. Biol. 24 , 8374 – 8385 . Google Scholar CrossRef Search ADS PubMed Wang , Y. ( 2007 ). Mitogen-activated protein kinases in heart development and diseases . Circulation 116 , 1413 – 1423 . Google Scholar CrossRef Search ADS PubMed Wang , J. , Li , C. , and Wang , E. ( 2010 ). Potential and flux landscapes quantify the stability and robustness of budding yeast cell cycle network . Proc. Natl Acad. Sci. USA 107 , 8195 – 8200 . Google Scholar CrossRef Search ADS Wang , W. , Xu , S. , Yin , M. , et al. . ( 2015 ). Essential roles of Gab1 tyrosine phosphorylation in growth factor-mediated signaling and angiogenesis . Int. J. Cardiol. 181 , 180 – 184 . Google Scholar CrossRef Search ADS PubMed Watts , D.J. , and Strogatz , S.H. ( 1998 ). Collective dynamics of ‘small-world’ networks . Nature 393 , 440 – 442 . Google Scholar CrossRef Search ADS PubMed Yancy , C.W. , Jessup , M. , Bozkurt , B. , et al. . ( 2013 ). 2013 ACCF/AHA guideline for the management of heart failure: a report of the American College of Cardiology Foundation/American Heart Association Task Force on practice guidelines . Circulation 128 , e240 – e327 . Google Scholar CrossRef Search ADS PubMed Yang , S. , Hu , S. , Hsieh , Y.C. , et al. . ( 2006 ). Mechanism of IL-6-mediated cardiac dysfunction following trauma-hemorrhage . J. Mol. Cell. Cardiol. 40 , 570 – 579 . Google Scholar CrossRef Search ADS PubMed Yang , J.H. , and Saucerman , J.J. ( 2011 ). Computational models reduce complexity and accelerate insight into cardiac signaling networks . Circ. Res. 108 , 85 – 97 . Google Scholar CrossRef Search ADS PubMed Zarubin , T. , and Han , J. ( 2005 ). Activation and signaling of the p38 MAP kinase pathway . Cell Res. 15 , 11 – 18 . Google Scholar CrossRef Search ADS PubMed Zhao , L. , Cheng , G. , Jin , R. , et al. . ( 2016 ). Deletion of Interleukin-6 attenuates pressure overload-induced left ventricular hypertrophy and dysfunction . Circ. Res. 118 , 1918 – 1929 . Google Scholar CrossRef Search ADS PubMed © The Author(s) (2018). Published by Oxford University Press on behalf of Journal of Molecular Cell Biology, IBCB, SIBS, CAS. All rights reserved. This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/about_us/legal/notices) http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Journal of Molecular Cell Biology Oxford University Press

Attractor landscape analysis of the cardiac signaling network reveals mechanism-based therapeutic strategies for heart failure

Loading next page...
 
/lp/ou_press/attractor-landscape-analysis-of-the-cardiac-signaling-network-reveals-KqLl05eM0y
Publisher
Oxford University Press
Copyright
© The Author(s) (2018). Published by Oxford University Press on behalf of Journal of Molecular Cell Biology, IBCB, SIBS, CAS. All rights reserved.
ISSN
1674-2788
eISSN
1759-4685
D.O.I.
10.1093/jmcb/mjy019
Publisher site
See Article on Publisher Site

Abstract

Abstract Apoptosis and hypertrophy of cardiomyocytes are the primary causes of heart failure (HF), a global leading cause of death, and are regulated through the complicated intracellular signaling network, limiting the development of effective treatments due to its complexity. To identify effective therapeutic strategies for HF at a system level, we develop a large-scale comprehensive mathematical model of the cardiac signaling network by integrating all available experimental evidence. Attractor landscape analysis of the network model identifies distinct sets of control nodes that effectively suppress apoptosis and hypertrophy of cardiomyocytes under ischemic or pressure overload-induced HF, the two major types of HF. Intriguingly, our system-level analysis suggests that intervention of these control nodes may increase the efficacy of clinical drugs for HF and, of most importance, different combinations of control nodes are suggested as potentially effective candidate drug targets depending on the types of HF. Our study provides a systematic way of developing mechanism-based therapeutic strategies for HF. heart failure, mathematical modeling, cardiac signaling network, systems analysis, disease mechanism, mechanism-based therapeutic strategy, systems biology Introduction Heart failure (HF) is a clinical syndrome characterized by typical symptoms such as shortness of breathing and peripheral swelling resulting from a reduced cardiac output (Remme and Swedberg, 2001; Hilfiker-Klemer et al., 2006). The disease is caused by a structural and/or functional cardiac abnormality that arises from myocardial ischemia (e.g. angina, myocardial infarction) or sustained pressure overload (e.g. hypertension) (Remme and Swedberg, 2001; Hilfiker-Klemer et al., 2006). It is one of the global leading causes of death, and its incidence and prevalence are increasing over time (Stewart et al., 2003). To date, several drugs are available for HF treatment that are effective at an early stage in clinics (ACE Inhibitor Myocardial Infarction Collaborative Group, 1998), but yet they are reported to have limited efficacy in blocking the disease progression in the long-term treatment, which highlights the urgent need to develop effective therapeutic strategies. Apoptosis and hypertrophy of cardiomyocytes are the main causes of HF and they are the typically observed molecular phenotypes in cardiomyocytes during the disease progression (Remme and Swedberg, 2001). Cardiomyocyte apoptosis occurs in diverse pathological conditions, including ischemia or pressure overload, and causes the permanent loss of the functioning units in the myocardium, thereby contributing to the development of HF (Narula et al., 1996; Gonzalez et al., 2003). Meanwhile, cardiomyocyte hypertrophy occurs to compensate for the loss of myocytes following injury and to increase the cardiac pumping function temporally (Kang et al., 2016). However, over time, myocardial hypertrophy induces pathological changes in the heart that eventually lead to a functional deterioration, which also promotes HF (Nadal-Ginard et al., 2003). Apoptosis and hypertrophy of cardiomyocytes are known to be regulated by the complicated interactions in the intracellular cardiac signaling network (Ryall et al., 2012; Kang et al., 2016). Previous experimental researches were conducted to understand how the cardiac signal transduction pathways are regulated and how they are related to pathological phenotypes of cardiomyocytes (Saraste et al., 1997; Bernardo et al., 2010). However, there is a fundamental limitation in applying conventional experimental approaches to investigate the underlying mechanisms of HF because the cardiac signaling network involves complicated feedback regulation and crosstalk (Kitano, 2002), which makes it difficult to intuitively understand the cellular behavior and develop effective treatments for HF patients. To overcome such limitations and investigate effective therapeutic strategies for HF treatment, we have employed systems biological approaches based on mathematical modeling and simulation analysis that have been proven to facilitate the investigation of the complex underlying mechanisms of cellular signaling networks (Kitano, 2002; Kim et al., 2007; Murray et al., 2010). To this end, we have constructed a very large-scale comprehensive cardiac signaling network and developed its mathematical model by integrating all relevant experimental information that is closely associated with cardiomyocytes. Attractor landscape analysis, which aims to understand a cellular behavior by investigating a specific converging network state determined by the interaction of signaling molecules, has been widely used to explore the dynamics of a complex signaling network (Choi et al., 2012; Cho et al., 2016). From attractor landscape analysis of the cardiac signaling network model, we have identified distinct sets of control nodes that can effectively suppress apoptosis and hypertrophy of cardiomyocytes under ischemic or pressure overload-induced HF condition, the two major types of HF. Intriguingly, we found that intervention of these control nodes may increase the efficacy of currently available clinical drugs for HF patients. Most importantly, different combinations of control nodes are suggested as potentially effective candidate drug targets depending on the types of HF. Our study highlights the significance of network-level analyses for unraveling molecular regulatory mechanisms, and provides a systematic way of developing mechanism-based therapeutic strategies for HF. Results A mathematical model of the cardiac signaling network To understand the complicated dynamics of the underlying intracellular signaling process of cardiomyocytes, we first reconstructed a very large-scale comprehensive cardiac signaling network by integrating all relevant evidence on molecules and their interactions through an extensive search over literature and database such as Kyoto Encyclopedia of Genes and Genomes (KEGG) and Pathway Interaction Database (PID) (Kanehisa and Goto, 2000; Schaefer et al., 2009). The reconstructed cardiac signaling network, the most comprehensive cardiac signaling network to the best of our knowledge at present, consists of 141 nodes and 274 directed links, including 17 external input nodes and 15 output nodes (Figure 1 and Supplementary Dataset S1). The input nodes in the cardiac signaling network represent diverse external stimuli that activate downstream signaling pathways to transmit signals from the outer surface of cardiomyocytes to cell interior and therefore generate appropriate cellular responses (Figure 1 and Supplementary Dataset S1). The network involves various signaling pathways, including α- or β-adrenergic (AR) (Kanehisa and Goto, 2000), calcium (Ca) (Kehat and Molkentin, 2010), cyclic guanosine monophosphate (cGMP)-protein kinase G (PKG) (Francis, 2010), mitogen-activated protein kinase (MAPK) (Wang, 2007), phosphatidylinositol 3-kinase (PI3K)-Akt (Matsui and Rosenzweig, 2005), transforming growth factor beta (TGFβ) (Bujak and Frangogiannis, 2007), Src-focal adhesion kinase (FAK) (Kuramochi et al., 2006), Janus kinase (JAK)-signal transducer and activator of transcription (STAT) (Pan et al., 1999), tumor necrosis factor alpha (TNFα) (Condorelli et al., 2002), and nuclear factor kappa-light-chain-enhancer of activated B cells (NFκB) pathways (Ling et al., 2013), which play important roles in coordinating the cellular functions of cardiomyocytes. These signaling pathways relay and integrate biochemical signals, and transmit them to the nucleus to activate enzymes and transcription factors that regulate apoptosis or hypertrophy of cardiomyocytes. In the cardiac signaling network, caspase 3 (Casp3) was used as an indicator of apoptosis, while the remaining 14 output nodes were taken as indicators of hypertrophy on the basis of their crucial roles for the regulation of major cellular functions in cardiomyocytes (Figure 1 and Supplementary Table S1). Of note, several nodes that are involved multiple signaling pathways function to mediate crosstalks between different cascades and they are known to play a key role in coordinating cellular responses of cardiomyocytes under diverse pathophysiological conditions (Heineke and Molkentin, 2006; Kang et al., 2016) (Figure 1). We further found that the cardiac signaling network has a higher clustering coefficient and neighborhood connectivity than random networks (Supplementary Table S2) and that the in- and out-degree distributions of nodes in the network follow an approximate power-law distribution (Supplementary Figure S1), which are in accord with the general characteristics of biological networks (Watts and Strogatz, 1998; Barabasi and Oltvai, 2004; Barabasi et al., 2011; Ryall et al., 2012; Lu et al., 2015). Figure 1 View largeDownload slide A schematic diagram of the cardiac signaling network. The reconstructed cardiac signaling network consists of 141 nodes and 274 links; 212 links are activation links (red pointed arrows) and 62 links are inhibition links (blue blunted arrows). Of the 141 nodes, there are 17 input nodes and 15 output nodes. Nodes that are involved in two or more signaling pathways according to the KEGG classification are highlighted with multiple colors that represent the corresponding pathways. Detailed information for each link and assigned logic table for each node is shown in Supplementary Datasets S1 and S2. Figure 1 View largeDownload slide A schematic diagram of the cardiac signaling network. The reconstructed cardiac signaling network consists of 141 nodes and 274 links; 212 links are activation links (red pointed arrows) and 62 links are inhibition links (blue blunted arrows). Of the 141 nodes, there are 17 input nodes and 15 output nodes. Nodes that are involved in two or more signaling pathways according to the KEGG classification are highlighted with multiple colors that represent the corresponding pathways. Detailed information for each link and assigned logic table for each node is shown in Supplementary Datasets S1 and S2. To validate whether our network model components are all expressed in cardiomyocytes, we have further evaluated the gene or protein expression levels of the network components using mRNA-sequencing or biochemical experimental data obtained from previous studies (Nakamura et al., 2000; Cuenca et al., 2006; Yang et al., 2006; Loan Le et al., 2012; Song et al., 2012; Ucar et al., 2012; Zhao et al., 2016). Note that for the case of mRNA-sequencing data analysis, genes whose reads per kilobase per million mapped reads (RPKM) values are greater than 1 were determined to be expressed as in the previous studies (Phanstiel et al., 2011; Handler et al., 2013). As a result, we found that all the nodes (except input ligands that are derived from extracellular microenvironment of cardiomyocytes, and the chemical elements that are not transcriptionally or translationally expressed/activated such as Ca, cyclic adenosine monophosphate (cAMP), cGMP, and inositol trisphosphate (IP3)) are specifically expressed in cardiomyocytes (Supplementary Table S3), further supporting that the reconstructed network can comprehensively represent cardiomyocyte-specific signaling mechanisms. Due to the limited experimental information on detailed reaction kinetics of the interactions between biomolecules, continuous dynamic models (e.g. an ordinary differential equation-based model) have difficulties in estimating their parameter values (Park et al., 2006; Helikar et al., 2008). To avoid such a problem, parameter-free discrete Boolean network modeling approach has been widely used in modeling various biological systems, which is proven to be highly useful in investigating key dynamic properties of complex intracellular signaling networks (Helikar et al., 2008). This modeling method has been shown to successfully capture the well-known biological features and provides reliable predictions of the potential effects of drug treatments (Choi et al., 2012; Lee et al., 2015). Therefore, we have developed a large-scale Boolean model of the cardiac signaling network based on the qualitative properties of activation and/or inhibition of the signaling proteins in cardiomyocytes (Supplementary Dataset S2). In the cardiac signaling network model, the state value of a distinct node represents its activity, and is reflected as ‘1’ or ‘0’ for each active or inactive state, respectively. Logic tables for nodes were created based on experimental information from the relevant literature and database search (Supplementary Datasets S1 and S2). To investigate the optimal input setting that can reasonably represent the physiological condition of cardiomyocytes, we have performed extensive simulations for the analysis of input–output relationships of network components (see Materials and methods). To this end, 100 million sets of randomly sampled input intensities were applied to the network and the resulting cellular responses were monitored (see Materials and methods). As a result, the intensity of input nodes that can most closely reproduce the well-known qualitative activities of network components was determined as an optimal input setting for the physiological condition (Supplementary Table S4). The simulation results of the input–output relationship analysis are shown in Figure 2, Supplementary Figures S2 and S3, and Table S5, and these results indicate that the cardiac signaling network model can well replicate the cellular responses of cardiomyocytes that are experimentally observed (Aikawa et al., 1997; Condorelli et al., 2002; Chen et al., 2003; Nakaoka et al., 2003, 2007; Fu et al., 2004; Vega et al., 2004; Zarubin and Han, 2005; Schaub and Hefti, 2007; Colella et al., 2008; Fischer and Hilfiker-Kleiner, 2008; Kawasaki et al., 2008; Guo and Wang, 2009; Ronca et al., 2009; Backs et al., 2011; Liu et al., 2012; Pennanen et al., 2014; Roberge et al., 2014; Wang et al., 2015; Ponnusamy et al., 2017). Figure 2 View largeDownload slide Qualitative input–output relationships in the cardiac signaling network model using last 100 average steps for the activity calculation. (A) Positive relationship between TNFα and Casp8 activation (Roberge et al., 2014). (B) Positive relationship between TNFα and p38 activation (Chen et al., 2003). (C) Positive relationship between TNFα and ANP expression (Condorelli et al., 2002). (D) Positive relationship between TNFα and cJun activation (Condorelli et al., 2002). (E) Positive relationship between NE and Ca activation (Pennanen et al., 2014). (F) Positive relationship between NE and NFAT activation (Colella et al., 2008). (G) Positive relationship between NE and MEF2 activation (Backs et al., 2011). (H) Positive relationship between NE and Casp3 activation (Fu et al., 2004). (I) Positive relationship between NRG1 and SHP2 activation (Nakaoka et al., 2007). (J) Positive relationship between PGE2 and PKA activation (Liu et al., 2012). (K) Positive relationship between LIF and GAB1 activation (Nakaoka et al., 2003). (L) Positive relationship between LIF and STAT activation (Schaub and Hefti, 2007). (M) Positive relationship between TGFβ and Ras activation (Guo and Wang, 2009). (N) Positive relationship between TGFβ and p38 activation (Zarubin and Han, 2005). (O) Positive relationship between TGFβ and BNP expression (Vega et al., 2004). (P) Positive relationship between TGFβ and ANP expression (Vega et al., 2004). (Q) Positive relationship between FGF and Ras activation (Ronca et al., 2009). (R) Positive relationship between AngII and Src activation (Aikawa et al., 1997). (S) Positive relationship between CT1 and JAK activation (Fischer and Hilfiker-Kleiner, 2008). (T) Positive relationship between CT1 and Ras activation (Fischer and Hilfiker-Kleiner, 2008). (U) Positive relationship between HGF and SHP2 activation (Wang et al., 2015). (V) Positive relationship between VEGF and Ras activation (Kawasaki et al., 2008). (W) Negative relationship between FGF and FOXO activation (Ponnusamy et al., 2017). (X) Positive relationship between EGF and GAB1 activation (Nakaoka et al., 2007). The average ‘1’s over last 100 time steps were calculated and used for stable state activity on each node. The dose–response curve presented here is intended to illustrate how the cardiac signaling network model can qualitatively reproduce the referenced input–output relationship for a wide range of input signals. The simulation was repeated (n = 10) for the input at 2% noise (see Materials and methods). Figure 2 View largeDownload slide Qualitative input–output relationships in the cardiac signaling network model using last 100 average steps for the activity calculation. (A) Positive relationship between TNFα and Casp8 activation (Roberge et al., 2014). (B) Positive relationship between TNFα and p38 activation (Chen et al., 2003). (C) Positive relationship between TNFα and ANP expression (Condorelli et al., 2002). (D) Positive relationship between TNFα and cJun activation (Condorelli et al., 2002). (E) Positive relationship between NE and Ca activation (Pennanen et al., 2014). (F) Positive relationship between NE and NFAT activation (Colella et al., 2008). (G) Positive relationship between NE and MEF2 activation (Backs et al., 2011). (H) Positive relationship between NE and Casp3 activation (Fu et al., 2004). (I) Positive relationship between NRG1 and SHP2 activation (Nakaoka et al., 2007). (J) Positive relationship between PGE2 and PKA activation (Liu et al., 2012). (K) Positive relationship between LIF and GAB1 activation (Nakaoka et al., 2003). (L) Positive relationship between LIF and STAT activation (Schaub and Hefti, 2007). (M) Positive relationship between TGFβ and Ras activation (Guo and Wang, 2009). (N) Positive relationship between TGFβ and p38 activation (Zarubin and Han, 2005). (O) Positive relationship between TGFβ and BNP expression (Vega et al., 2004). (P) Positive relationship between TGFβ and ANP expression (Vega et al., 2004). (Q) Positive relationship between FGF and Ras activation (Ronca et al., 2009). (R) Positive relationship between AngII and Src activation (Aikawa et al., 1997). (S) Positive relationship between CT1 and JAK activation (Fischer and Hilfiker-Kleiner, 2008). (T) Positive relationship between CT1 and Ras activation (Fischer and Hilfiker-Kleiner, 2008). (U) Positive relationship between HGF and SHP2 activation (Wang et al., 2015). (V) Positive relationship between VEGF and Ras activation (Kawasaki et al., 2008). (W) Negative relationship between FGF and FOXO activation (Ponnusamy et al., 2017). (X) Positive relationship between EGF and GAB1 activation (Nakaoka et al., 2007). The average ‘1’s over last 100 time steps were calculated and used for stable state activity on each node. The dose–response curve presented here is intended to illustrate how the cardiac signaling network model can qualitatively reproduce the referenced input–output relationship for a wide range of input signals. The simulation was repeated (n = 10) for the input at 2% noise (see Materials and methods). Reflecting molecular characteristics under HF conditions to the cardiac signaling network model The aim of this study is to investigate most effective therapeutic strategies for HF depending on conditions. Although HF is caused by multiple pathologies (Hilfiker-Klemer et al., 2006), myocardial ischemia or sustained pressure overload is known to be its leading cause according to previous studies and clinical evidence (Blankesteijn et al., 2008; Thygesen et al., 2012). For this reason, we focused on investigating the signaling mechanisms of the two major types of HF (i.e. ischemic HF and pressure overload-induced HF). The heart function is regulated by various growth factors (i.e. insulin-like growth factor 1 (IGF1), epidermal growth factor (EGF), vascular endothelial growth factor (VEGF), TGFβ, and fibroblast growth factor (FGF)), catecholamines (i.e. norepinephrine (NE) and phenylephrine (PE)) or cytokines (i.e. TNFα, endothelin-1 (ET1), and interleukin (IL)) of which the secretion is modulated by autonomic nerve or endocrine system (Iwamoto et al., 2003; Florea and Cohn, 2014). During the development of HF due to the cardiac dysfunction, the sympathetic nervous system or the adrenal gland is activated as an early compensatory mechanism that supports inotropic action and maintains cardiac output (Jackson et al., 2000). This activated sympathetic nervous system or the adrenal gland changes the secretion profiles of various growth factors, catecholamines, or cytokines, which show different aspects from those of physiological condition (Kang et al., 2016). The changes in the secretion of those extracellular stimuli augment cardiac contractility and relaxation that increase the heart pumping function in the early phase of HF (Nadal-Ginard et al., 2003). However, in the long run, these stimuli activate various signaling pathways that play pivotal roles in the regulation of cardiomyocyte apoptosis and hypertrophy, which eventually contribute to the progression of HF (MacLellan and Schneider, 2000). Considering that the change in stimuli secretion profiles is one of the major causes of cardiomyocyte apoptosis and hypertrophy (MacLellan and Schneider, 2000), we have reflected ischemic or pressure overload-induced HF condition to the cardiac signaling network model by applying the distinct input intensities of individual HF condition. For this purpose, we integrated all relevant information on upregulated or downregulated molecular activities in each HF condition through an extensive search over literature and database (Supplementary Table S6). These collected observations include not only the changed activity profiles of the external ligands but also those of intracellular signaling molecules and transcription factors under HF conditions. Based on the fact that the changed stimuli secretion profiles in HF conditions cause the changes in the activity of downstream molecules (MacLellan and Schneider, 2000), we have performed extensive simulation analyses while randomly increasing (or decreasing) the intensity of input nodes of which the activity is enhanced (or suppressed) under distinct HF conditions in a wide range for 100 million times (Supplementary Table S6), and monitored the activity profiles of downstream signaling components in the cardiac signaling network. As a result, we found the optimal input setting that can reasonably represent distinct HF conditions of cardiomyocytes (Figure 3). To determine the robustness of the optimized input setting for HF conditions, additional computer simulations were performed while randomly perturbing the intensity of upregulated or downregulated inputs in the range of 5% or 20% for 100 million times. As a result, activity profiles of the molecules were well reproduced regardless of the perturbation range (Supplementary Figure S4). This result indicates that the mathematical model is adequate for studying the molecular mechanisms of the cardiac signaling network under HF conditions. Figure 3 View largeDownload slide The identified optimal input setting representing ischemic or pressure overload-induced HF condition of cardiomyocytes. The optimal input setting of cardiomyocytes was determined for ischemic or pressure overload-induced HF condition that reasonably reproduces qualitative characteristics of the activity profile for the network components under distinct HF conditions. Detailed information about the optimal input setting is shown in the main text. Blue or red nodes indicate differentially activated proteins (cutoffs: fold change >1.5, P < 0.05) under ischemic or pressure overload-induced HF, respectively. Green nodes indicate upregulated components in signaling pathways that contain the differentially activated proteins. Figure 3 View largeDownload slide The identified optimal input setting representing ischemic or pressure overload-induced HF condition of cardiomyocytes. The optimal input setting of cardiomyocytes was determined for ischemic or pressure overload-induced HF condition that reasonably reproduces qualitative characteristics of the activity profile for the network components under distinct HF conditions. Detailed information about the optimal input setting is shown in the main text. Blue or red nodes indicate differentially activated proteins (cutoffs: fold change >1.5, P < 0.05) under ischemic or pressure overload-induced HF, respectively. Green nodes indicate upregulated components in signaling pathways that contain the differentially activated proteins. To further investigate whether the signaling dynamics are differentially regulated under ischemic or pressure overload-induced HF condition, we performed pathway analysis of the signaling proteins that are significantly differentially activated under each HF condition using g:Profiler (Reimand et al., 2016). As a result, the differentially activated proteins (cutoffs: fold change >1.5, P < 0.05) between each individual HF type were found to be significantly enriched in different signaling pathways (Figure 3 and Table 1). These simulation results are consistent with previous experiments: for instance, adenylyl cyclase (AC)-cAMP signaling promotes the progression of ischemic HF (Appukuttan et al., 2012; Bravo et al., 2016), while TGFβ signaling is found to be important in the development of pressure overload-induced HF (Lim et al., 2005; Koitabashi et al., 2011). This further supports that the underlying signaling mechanisms may differ depending on the HF types and, therefore, different signaling pathways should be targeted to achieve optimal therapeutic responses, which can be useful for the development of mechanism-based treatment strategy for HF. Table 1 Pathway analysis of differentially activated proteins between ischemic and pressure overload-induced HF. Category Term Ischemic HF Pressure overload-induced HF Count Percent P-value Count Percent P-value KEGG Adrenergic signaling in cardiomyocytes 10 6.9 <0.001 0 0.0 N.S. Calcium signaling pathway 9 5.0 <0.001 1 0.6 N.S. cAMP signaling pathway 8 4.0 <0.001 1 0.5 N.S. cGMP–PKG signaling pathway 8 4.9 <0.001 2 1.2 N.S. Jak–STAT signaling pathway 5 3.2 0.020 0 0.0 N.S. Ras signaling pathway 6 2.6 0.013 3 1.3 N.S Renin–angiotensin system 3 13 0.006 0 0.0 N.S. Focal adhesion 2 1.0 N.S. 6 3.0 <0.001 FoxO signaling pathway 0 0.0 N.S. 5 3.8 <0.001 MAPK signaling pathway 5 2.0 N.S. 6 2.4 <0.001 TGF-β signaling pathway 1 1.2 N.S. 6 7.1 <0.001 Category Term Ischemic HF Pressure overload-induced HF Count Percent P-value Count Percent P-value KEGG Adrenergic signaling in cardiomyocytes 10 6.9 <0.001 0 0.0 N.S. Calcium signaling pathway 9 5.0 <0.001 1 0.6 N.S. cAMP signaling pathway 8 4.0 <0.001 1 0.5 N.S. cGMP–PKG signaling pathway 8 4.9 <0.001 2 1.2 N.S. Jak–STAT signaling pathway 5 3.2 0.020 0 0.0 N.S. Ras signaling pathway 6 2.6 0.013 3 1.3 N.S Renin–angiotensin system 3 13 0.006 0 0.0 N.S. Focal adhesion 2 1.0 N.S. 6 3.0 <0.001 FoxO signaling pathway 0 0.0 N.S. 5 3.8 <0.001 MAPK signaling pathway 5 2.0 N.S. 6 2.4 <0.001 TGF-β signaling pathway 1 1.2 N.S. 6 7.1 <0.001 Table showing the pathway analysis results of differentially activated proteins (cutoffs: fold change >1.5, P < 0.05) between ischemic and pressure overload-induced HF using g:Profiler (Reimand et al., 2016). N.S., no significance. Table 1 Pathway analysis of differentially activated proteins between ischemic and pressure overload-induced HF. Category Term Ischemic HF Pressure overload-induced HF Count Percent P-value Count Percent P-value KEGG Adrenergic signaling in cardiomyocytes 10 6.9 <0.001 0 0.0 N.S. Calcium signaling pathway 9 5.0 <0.001 1 0.6 N.S. cAMP signaling pathway 8 4.0 <0.001 1 0.5 N.S. cGMP–PKG signaling pathway 8 4.9 <0.001 2 1.2 N.S. Jak–STAT signaling pathway 5 3.2 0.020 0 0.0 N.S. Ras signaling pathway 6 2.6 0.013 3 1.3 N.S Renin–angiotensin system 3 13 0.006 0 0.0 N.S. Focal adhesion 2 1.0 N.S. 6 3.0 <0.001 FoxO signaling pathway 0 0.0 N.S. 5 3.8 <0.001 MAPK signaling pathway 5 2.0 N.S. 6 2.4 <0.001 TGF-β signaling pathway 1 1.2 N.S. 6 7.1 <0.001 Category Term Ischemic HF Pressure overload-induced HF Count Percent P-value Count Percent P-value KEGG Adrenergic signaling in cardiomyocytes 10 6.9 <0.001 0 0.0 N.S. Calcium signaling pathway 9 5.0 <0.001 1 0.6 N.S. cAMP signaling pathway 8 4.0 <0.001 1 0.5 N.S. cGMP–PKG signaling pathway 8 4.9 <0.001 2 1.2 N.S. Jak–STAT signaling pathway 5 3.2 0.020 0 0.0 N.S. Ras signaling pathway 6 2.6 0.013 3 1.3 N.S Renin–angiotensin system 3 13 0.006 0 0.0 N.S. Focal adhesion 2 1.0 N.S. 6 3.0 <0.001 FoxO signaling pathway 0 0.0 N.S. 5 3.8 <0.001 MAPK signaling pathway 5 2.0 N.S. 6 2.4 <0.001 TGF-β signaling pathway 1 1.2 N.S. 6 7.1 <0.001 Table showing the pathway analysis results of differentially activated proteins (cutoffs: fold change >1.5, P < 0.05) between ischemic and pressure overload-induced HF using g:Profiler (Reimand et al., 2016). N.S., no significance. Attractor landscape analysis of the cardiac signaling network model To investigate the dynamic behavior of the intracellular signaling process of cardiomyocytes and to identify effective molecular targets depending on the types of HF, we applied the state space analysis to the cardiac signaling network model, which has been proved to be effective for investigating the dynamic characteristics of a cellular system (Choi et al., 2012; Kim et al., 2014). It is known that complex signaling networks exhibit ordered (stable) dynamics, and that distinct cellular phenotypes can be described by attractor states (Choi et al., 2012; Kim et al., 2014). Attractors are stable states in which a dynamic system (i.e. the cardiac signaling network) tends to converge over time (Choi et al., 2012). The dynamic characteristics of a complex signaling network could be mapped to an attractor landscape (Lee et al., 2015). Every point in the attractor landscape corresponds to a specific network state which is defined by the activity state profile of every molecule in the network (Huang et al., 2009). The state space of the Boolean network with N nodes consists of 2N different states and the dynamics of the state transitions are characterized by the flow of sequential states which eventually converge towards attractors (Lee et al., 2015). The attractor is a particular subset of states, either a fixed point (a point attractor) which remains stable as a single network state, or a limited cycle of period p (a cyclic attractor), consisting of p states that are repeatedly visited by the network dynamics (Huang et al., 2005). The set of initial states with trajectories converging to a given attractor is the basin of attraction (Huang et al., 2005). Based on the state space analysis, distinct attractors of the cardiac signaling network model were classified into three groups characterized by particular cellular phenotypes. The apoptosis or hypertrophy phenotype was defined as an attractor in which any of output nodes used as indicators for corresponding phenotypes is ‘ON’ at least once in its cyclic state transition. The remaining attractors where all the output nodes maintain ‘OFF’ state were determined as a normal phenotype. During the attractor landscape analysis, physiological, ischemic HF, or pressure overload-induced HF condition was reflected to the network model by assigning the state value of particular input nodes to ‘1’ with an intensity determined in its corresponding input setting (Figure 3) while sampling a given number of initial states. From attractor landscape analysis of the cardiac signaling network model, we revealed that the basin sizes of attractors for normal phenotype were highly decreased whereas those for apoptosis or hypertrophy were dramatically increased under HF conditions compared to the case of physiological condition (Supplementary Figure S5 and Table S7). We note that distributions of the relative basin sizes of attractors are maintained regardless of the sampling size of the initial states (Supplementary Figures S5 and S6). In summary, the results suggest that apoptosis or hypertrophy of cardiomyocytes occurs through the coordinated activation of various pathways in the cardiac signaling network under HF conditions. Identification of distinct control nodes that effectively suppress apoptosis and hypertrophy of cardiomyocytes under HF conditions The cellular behavior is controlled by complicated regulatory interactions of various molecules that constitute the intracellular signaling network (Barabasi and Oltvai, 2004). Previous studies showed that intervention of only a small portion of the molecules, which we call control nodes, is needed to govern such complicated network dynamics and to control toward a desired cellular phenotype (Kim et al., 2011, 2013). To identify a set of nodes that can effectively suppress apoptosis and hypertrophy of cardiomyocytes under ischemic or pressure overload-induced HF condition, we employed genetic algorithm (GA) which is widely used to solve complex optimization problems (see Materials and methods for details). As a result, distinct sets of control nodes were found to be important for effectively suppressing apoptosis and hypertrophy of cardiomyocytes: G proteins (q) alpha subunit and G proteins (γ) alpha subunit (Gαq/11), TNF receptor associated factor (TRAF), PI3K, AC, and the growth factor receptor-bound protein 2 (Grb2) and Son of Sevenless (SOS) complex (G/S) under ischemic HF condition and TGFβ receptor (TGFR), Integrins, PI3K, Gαq/11, TRAF, Ras, and G proteins (s) alpha subunit (Gsα) under pressure overload-induced HF condition. These model predictions are supported by previous experimental findings: for instance, inhibition of AC, an identified control node under ischemic HF condition, was shown to suppress apoptosis of cardiomyocytes in response to ischemia injury (Appukuttan et al., 2012; Bravo et al., 2016). Likewise, myocyte knockdown of TGFR, an identified control node under pressure overload-induced HF condition, significantly reduced maladaptive hypertrophy of cardiomyocytes in transverse aortic constriction (an experimental model used for inducing pressure overload-induced cardiac hypertrophy and apoptosis in animals)-operated mice (Koitabashi et al., 2011). We further found that simultaneous inhibition of the control nodes highly decreases the basin sizes of attractors for apoptosis and hypertrophy under HF conditions, implying that the activities of these nodes should be negatively controlled to exhibit therapeutic effects (Figure 4A). To investigate the correlation between the number of controlled nodes and the basin sizes of attractors for normal phenotype, apoptosis, and hypertrophy, we have investigated the cellular behavior while controlling activities of different numbers of randomly sampled nodes. As a result, we found that the maximal basin sizes of attractors for normal phenotype were increased when a more number of sampled nodes were controlled, but controlling a few molecular targets was sufficient to regulate the cellular responses of cardiomyocytes (Supplementary Figure S7). Taken together, these results suggest that only a small fraction of network nodes is required to be controlled to obtain desired network states, which is consistent with the previous findings in other biological networks (Liu et al., 2011; Kim et al., 2013). Figure 4 View largeDownload slide Effect of inhibition of the control nodes on the cellular responses of cardiomyocytes under HF conditions. (A) Distribution of the basin sizes of attractors before and after inhibition of the control nodes under ischemic or pressure overload-induced HF condition. The colored boxes in the right panel describe the classified attractors representing different cellular phenotypes. White or black boxes represent cell phenotype as ‘off’ or ‘on’, respectively. Detailed information about the attractor classification is shown in the main text. (B−D) Distribution of the basin sizes of attractors for normal phenotype (B), apoptosis (C), or hypertrophy (D) after inhibiting the activities of a fixed number of sampled control nodes or random nodes under ischemic or pressure overload-induced HF condition. Simulation analysis was performed while controlling activities of possible combinations of control nodes or randomly sampled nodes and was repeated for 10 times using different random seeds of one million initial states following the same procedure. Data represent mean + SEM of all simulation results. *P < 0.05; **P < 0.01; ***P < 0.001; ****P < 0.0001; Student’s t-test. Figure 4 View largeDownload slide Effect of inhibition of the control nodes on the cellular responses of cardiomyocytes under HF conditions. (A) Distribution of the basin sizes of attractors before and after inhibition of the control nodes under ischemic or pressure overload-induced HF condition. The colored boxes in the right panel describe the classified attractors representing different cellular phenotypes. White or black boxes represent cell phenotype as ‘off’ or ‘on’, respectively. Detailed information about the attractor classification is shown in the main text. (B−D) Distribution of the basin sizes of attractors for normal phenotype (B), apoptosis (C), or hypertrophy (D) after inhibiting the activities of a fixed number of sampled control nodes or random nodes under ischemic or pressure overload-induced HF condition. Simulation analysis was performed while controlling activities of possible combinations of control nodes or randomly sampled nodes and was repeated for 10 times using different random seeds of one million initial states following the same procedure. Data represent mean + SEM of all simulation results. *P < 0.05; **P < 0.01; ***P < 0.001; ****P < 0.0001; Student’s t-test. To further evaluate whether the control nodes could effectively suppress apoptosis and hypertrophy of cardiomyocytes, we have investigated the cellular response while inhibiting all possible combinations of a fixed number of sampled control nodes under HF conditions. As a result, we found that the inhibition of identified control nodes significantly decreased the basin sizes of attractors for apoptosis and hypertrophy, but increased those for normal phenotype compared to the cases when the same number of nodes was randomly chosen (Figure 4B–D). Moreover, the efficacy for regulating the cellular behavior was improved when an increasing number of control nodes were inhibited (Figure 4B–D). These results were still consistent under the given model perturbations (Supplementary Figure S8, see Materials and methods for details). In summary, the result suggests that the identified sets of control nodes are important to effectively suppress apoptosis and hypertrophy of cardiomyocytes under HF conditions. Targeting control nodes increases the efficacy of currently available clinical drugs for HF treatment Previous studies suggested that nodes which are used as drug targets globally influence the cellular behavior (Spiro et al., 2008). Considering that the identified control nodes are found to be important for effectively suppressing apoptosis and hypertrophy of cardiomyocytes under HF conditions, we hypothesized that the control nodes might serve as drug targets. To this end, we have collected all approved or potential drug targets for HF from DrugBank database and literature search (Supplementary Table S8). Next, we compared the ratio of drug targets among the nodes of the control nodes with that of non-control nodes for the cardiac signaling network. Intriguingly, we found that the proportion of drug targets from control nodes was significantly higher than that from non-control nodes (Figure 5A). This suggests that the identified control nodes may be closely related to drug targets. Figure 5 View largeDownload slide Intervention of control nodes increases the efficacy of clinical drugs for HF treatment. (A) The fraction of drug targets in control nodes is significantly higher than that in non-control nodes under both ischemic and pressure overload-induced HF conditions. *P < 0.05; **P < 0.01; Fisher’s exact test. (B−D) Distribution of the increased (B) or decreased (C and D) basin sizes of attractors for normal phenotype (B), apoptosis (C), or hypertrophy (D) after inhibiting a fixed number of sampled control nodes–drug targets when combined with β-blocker treatment under ischemic or pressure overload-induced HF condition. Simulation analysis was performed while controlling activities of possible combinations of control nodes–drug targets or randomly sampled drug targets and was repeated for 10 times using different random seeds of one million initial states following the same procedure. Data represent mean + SEM of all simulation results. *P < 0.05; **P < 0.01; ***P < 0.001; ****P < 0.0001; Student’s t-test. See also Supplementary Table S9. Figure 5 View largeDownload slide Intervention of control nodes increases the efficacy of clinical drugs for HF treatment. (A) The fraction of drug targets in control nodes is significantly higher than that in non-control nodes under both ischemic and pressure overload-induced HF conditions. *P < 0.05; **P < 0.01; Fisher’s exact test. (B−D) Distribution of the increased (B) or decreased (C and D) basin sizes of attractors for normal phenotype (B), apoptosis (C), or hypertrophy (D) after inhibiting a fixed number of sampled control nodes–drug targets when combined with β-blocker treatment under ischemic or pressure overload-induced HF condition. Simulation analysis was performed while controlling activities of possible combinations of control nodes–drug targets or randomly sampled drug targets and was repeated for 10 times using different random seeds of one million initial states following the same procedure. Data represent mean + SEM of all simulation results. *P < 0.05; **P < 0.01; ***P < 0.001; ****P < 0.0001; Student’s t-test. See also Supplementary Table S9. Presently, several drugs (β-, α-, angiotensin II receptor (AT1R)-, and ET1 receptor (ET1R)-blockers) are clinically used for the HF patients, but their efficacy is often limited for blocking the disease progression (Rehsia and Dhalla, 2010). Since control nodes are largely overlapped with the drug targets, we further investigated whether the intervention of the control nodes that are drug targets (hereafter control nodes–drug targets) could improve the therapeutic effects of currently available clinical drugs. For this purpose, we have investigated the cellular response while inhibiting all possible combinations of a fixed number of sampled control nodes–drug targets when combined with each clinically available drug under HF conditions. Note that during the mathematical simulation, the effect of the clinical drugs (i.e. β-, α-, AT1R-, and ET1R-blockers) was reflected by pinning the state of their corresponding target nodes, β1/β2AR, αAR, AT1R, and ET1R, to ‘OFF’, respectively. As a result, inhibition of sampled control nodes–drug targets decreased the basin sizes of attractors for apoptosis and hypertrophy, but increased those for normal phenotype when combined with each clinical drug (Figure 5B–D and Supplementary Table S9). These results were still consistent under the given model perturbations (Supplementary Figures S9–S12, see Materials and methods for details). Intriguingly, optimal combinations of control nodes–drug targets that are most effective in increasing the basin sizes of attractors for normal phenotype were different depending on the types of HF or clinical drugs, suggesting the importance of a system-level development of mechanism-based therapeutic strategies for HF (Table 2). Table 2 The optimal combinations of control nodes for increasing the basin sizes of attractors for normal phenotype under each distinct HF condition. The optimal combinations of control nodes–drug targets that are most effective in increasing the basin sizes of attractors for normal phenotype when combined with each clinical drug under distinct HF condition. Black boxes indicate corresponding optimal combinations of control nodes. Table 2 The optimal combinations of control nodes for increasing the basin sizes of attractors for normal phenotype under each distinct HF condition. The optimal combinations of control nodes–drug targets that are most effective in increasing the basin sizes of attractors for normal phenotype when combined with each clinical drug under distinct HF condition. Black boxes indicate corresponding optimal combinations of control nodes. To further evaluate whether G/S or Gsα, the identified control node that is not included in drug targets, could also improve the therapeutic effects of currently available clinical drugs, we have investigated the cellular behavior while perturbing activities of G/S and Gsα when combined with each clinical drug under ischemic and pressure overload-induced HF conditions, respectively. As a result, we found that the inhibition of G/S or Gsα together with combinations of sampled control nodes–drug targets decreased the basin sizes of attractors for apoptosis and hypertrophy, but increased those for normal phenotype when combined with each clinical drug compared to the cases when the same number of drug targets was randomly chosen (Supplementary Table S10). This result suggests that G/S or Gsα may have a potential to serve as a therapeutic target for HF treatment. Discussion HF is a pathophysiological state in which the heart is unable to pump out a sufficient amount of blood throughout the body (Hilfiker-Klemer et al., 2006). Apoptosis and hypertrophy of cardiomyocytes are the leading causes of HF and they are the typically observed pathological phenotypes in cardiomyocytes during the progression of HF (Remme and Swedberg, 2001). A number of previous studies have used in silico simulation analysis to understand the complicated regulatory mechanisms in cardiomyocytes and to develop the effective therapeutic strategies for HF (Yang and Saucerman, 2011). However, those studies still exhibited limitations because their network models were established focusing only on either apoptosis or hypertrophy separately, and only a few relevant signaling pathways were partially considered in the models (Shin et al., 2011, 2014; Yang and Saucerman, 2011; Ryall et al., 2012; Kang et al., 2017). To overcome such limitations, we have developed a novel mathematical model of the cardiac signaling network, the most comprehensive network among currently existing models to the best of our knowledge, by collecting all available relevant experimental information, and investigated effective therapeutic strategies for HF treatment. Our novel findings obtained from these efforts led to the following conclusions: (i) we have identified distinct sets of control nodes that can effectively suppress apoptosis and hypertrophy of cardiomyocytes under ischemic or pressure overload-induced HF condition. (ii) Intervention of the control nodes may increase the efficacy of clinical drugs for HF. (iii) Different combinations of control nodes are suggested to serve as potentially effective candidate drug targets depending on the types of HF. In general, HF is not a single disease entity but can be considered as a wide disease spectrum including various pathological phenotypes (Ponikowski et al., 2016). In this regard, the underlying molecular pathophysiology can be different between HF patients although they have the same clinical features (De Keulenaer and Brutsaert, 2007). In other words, individual HF patient follows a unique disease trajectory depending on his/her clinical presentations, comorbidities, as well as the molecular pathophysiology including complex interplay of signaling dynamics at the molecular level (De Keulenaer and Brutsaert, 2007). However, the current guideline for HF treatment is largely based on the clinical features of the disease such as symptoms and signs rather than the molecular level understanding (Yancy et al., 2013; Ponikowski et al., 2016), and, therefore, the efficacy of currently available drugs considerably varies between HF patients and the drug treatment may even cause unwanted adverse effects (Yancy et al., 2013). Here, intriguingly, distinct combinations of control nodes were suggested as potentially effective candidate drug targets depending on HF types. This finding implies that individualized drug treatment strategies based on the molecular mechanism of HF can be used to achieve optimal therapeutic responses, underscoring the urgent need to develop mechanism-based therapeutic strategies for HF. We also note that the drugs that are widely used for HF patients in the clinics are mostly receptor antagonists (Law et al., 2014). Considering that the identified control nodes comprise not only receptors but also intracellular kinases, the results suggest a rationale for the development of the downstream kinase inhibitors for HF treatment (Law et al., 2014). Collectively, our findings may contribute to the improvement of the current guideline for HF treatment and the development of effective therapeutic strategies. Epigenetic landscape approach has been widely employed to study the complex regulatory dynamics of various biological phenomena (Wang et al., 2010; Li and Wang, 2013, 2014a, b; Li, 2017). This method attempts to uncover the functional mechanism of transition between distinct cellular states by analyzing the landscape topography and comparing the height of barriers that separate different basins of attraction in biological network systems (Wang et al., 2010; Li and Wang, 2013, 2014a, b; Li, 2017). Previous studies have proven the usefulness of landscape topography analysis in investigating biological processes and cellular functions since it provides a quantitative way of understanding the dynamic behavior of complex intracellular signaling networks (Wang et al., 2010; Li and Wang, 2013, 2014a, b; Li, 2017). Based on the epigenetic landscape concept and topography analysis of continuous mathematical models of specific gene regulatory systems, including stem cell differentiation, cancer, and cancer immunity (Li and Wang, 2013, 2014b; Li, 2017), they revealed key factors that are critically responsible for determining the barrier height between different attractors and identified optimal intervention strategies that may produce desired outcomes by controlling the barrier height (Li and Wang, 2013, 2014b; Li, 2017). In this study, we have developed a large-scale discrete Boolean model of cardiac signaling network and focused on the investigation of distinct sets of control nodes that can effectively maximize the basin sizes of attractors for normal phenotype under HF conditions, rather than analyzing the landscape topography of the network model. A future study employing the epigenetic landscape approach based on topography analysis of our cardiac signaling network model will further advance the understanding of regulatory dynamics of HF and may provide a rationale for a novel therapeutic strategy for HF treatment. In this study, we have developed a mathematical model of a comprehensive cardiac signaling network by integrating available cardiomyocyte-specific experimental evidence through an extensive search over various literatures and databases. However, there were some difficulties in establishing a large-scale network model by only using mechanical information obtained from such literatures and databases due to the intrinsic biological noise and false positive results of biochemical experiments, and the lack of sufficient available information about the cooperativity among the incoming regulatory signals of biomolecules (Macarthur et al., 2009). To overcome such challenge, we used GA, an optimization algorithm (Goldberg, 1989), in assigning the most appropriate logic tables to the nodes that have multiple positive and/or negative upstream regulators as such that the mathematical model can well reproduce the cellular responses of cardiomyocytes (see Materials and methods). Sufficient accumulation of biological data that can provide accurate mechanistic information on functional regulation of genes and proteins based on advanced experimental techniques would facilitate the development of more reliable mathematical model. Moreover, we mainly focused on apoptosis and hypertrophy of cardiomyocytes as key molecular phenotypic features of HF and investigated optimal therapeutic strategies that can effectively suppress such maladaptive cellular responses. However, myocardial fibrosis developed by complex heterocellular interactions between cardiomyocytes and other various cell types, including cardiac fibroblasts, is also known to play an important role in the HF progression (Kong et al., 2014; Leask, 2015; Gourdie et al., 2016; Travers et al., 2016). Development and analysis of mathematical models that can comprehensively encompass both intra- and intercellular signaling mechanisms of pathophysiology of HF may provide more improved therapeutic approaches. In summary, from attractor landscape analysis of the cardiac signaling network model, we have identified distinct sets of control nodes that can effectively suppress apoptosis and hypertrophy of cardiomyocytes under ischemic or pressure overload-induced HF condition. We further found that intervention of these control nodes may increase the efficacy of clinical drugs for HF treatment and that different combinations of control nodes are suggested as potentially effective candidate drug targets depending on the types of HF. Our study provides a systematic way of developing effective mechanism-based therapeutic strategies for HF. Materials and methods Development of the cardiac signaling network model The cardiac signaling network was reconstructed by incorporating all relevant information on the complicated interactions between individual proteins in cardiomyocytes through an extensive search over literature and database, and it comprises 141 nodes and 274 directed links. The 274 directed links comprising the cardiac signaling network are shown in Supplementary Dataset S1. The column ‘Source’ indicates the molecules that regulate downstream ‘Target’ molecules; and the column ‘Target’ indicates the molecules being activated or inhibited by the upstream ‘Source’ molecules. In the column ‘Interaction’, ‘+’ and ‘−’ represent activation and inhibition, respectively. The activity states of intracellular signaling nodes except for the input nodes were determined by their assigned logic tables as shown in Supplementary Dataset S2. The logic table for each node was determined based on the biological information on activation and/or inhibition of distinct molecules obtained from the related literature and database search (Supplementary Dataset S2). For the cases of nodes that have two or more positive and/or negative upstream regulators, however, there were some difficulties in generating their logic tables based only on the literature mining due to the lack of sufficient available information about the cooperativity among the incoming regulatory signals. Therefore, we used the weighted sum Boolean function to assign the most appropriate logic tables to the nodes with at least two positive and/or negative input nodes as in the previous study (Fumia and Martins, 2013). In this weighted sum modeling method, each node i has only two states, Si = 1 and Si = 0, representing the active and the inactive state of a node, respectively. The state values of nodes in the next time step (t + 1) are determined by the state values of their input nodes in the present time step (t) according to the following rule: Si(t+1)=sgn(∑j=1kin(i)JjiSj(t)−θi) Here, Jji is the interaction weight from input j on node i, where activating and inhibitory interactions are positive and negative numbers, respectively. kin(i) indicates the number of input nodes that regulate node i and determine its activation state. The threshold function sgn(x) is a unitary step function (sgn(x) = 1 if x > 0, sgn(x) = 0 if x ≤ 0). θi is the activation threshold of node i. Therefore, if a weighted sum of input signals of node i exceeds its activation threshold θi, the node is activated or stays active if it was already active; otherwise, it becomes inactive or stays inactive. We employed GA (Goldberg, 1989), which is widely used to solve complex optimization problems, to identify the optimal values of interaction weight Jji and activation threshold θi that can reasonably recapitulate the dynamic activity profiles of network components as reported from previous experiments (Figure 2). All the Boolean simulation analyses of the cardiac signaling network model were implemented using Matlab® 2014b in a Window Cluster system composed of 288 CPUs in parallel. In silico simulation analysis of input–output relationships of network components While the state space of the cardiac signaling network is huge (2141 states), we found that each network state eventually converges to a relatively small periodic cycle within 300 time steps from numerous iterative simulations. Thus, we generated Boolean function updates for 1000 time steps to measure the stable state activities of nodes using various levels of input stimulations. To represent the input intensity of each simulation run, we placed the input node in the cycle that yields the desired percentage of ON (Lee et al., 2015). For instance, if an input is set to 50% ON, it would be put on a cycle of ‘1010’ or ‘0101’, and this condition is fixed during the mathematical simulation. A small percentage of background noise (5%, described below) was added to the inputs to produce biologically reliable characteristics. For instance, if an input is set to 40% ON for a given run under 5% noise, the input changes randomly from 35% to 45% ON. The initial condition of each node (except input nodes) was randomized in every run (50% ‘0’, 50% ‘1’). The average ‘1’s over the last 50, 100, or 200 time steps were calculated and used to obtain the steady-state activity of each node. The results showed that the input–output relationships of the candidate nodes are similar regardless of the number of last average steps (Figure 2; Supplementary Figures S2 and S3). Note that, during the extensive simulation for the input–output relationship analysis of the cardiac signaling network model, the input intensity of a particular input node of interest was set to 0%−100% ON with 1% increments while intensities of remaining input nodes were fixed to the optimal input settings, as previously described (Helikar et al., 2008) (Supplementary Table S4). Attractor landscape analysis of the cardiac signaling network model To understand the cellular characteristics of the cardiac signaling network model, we performed the attractor landscape analysis (Choi et al., 2012). The attractor landscape is composed of attractors, the stable states of the network (Cho et al., 2016). The set of network states that converge to the corresponding attractors is called the basin of attraction (Huang et al., 2005). For a given input condition, we obtained the attractor landscape of the cardiac signaling network model using one million randomly sampled initial states until the state trajectory of an initial state reaches a point attractor with a single network state or a cyclic attractor with repetitive states (Choi et al., 2012; Cho et al., 2016). Note that the simulation result was maintained regardless of the sampling size of initial states (Supplementary Figures S5 and S6). During the attractor landscape analysis, physiological, ischemic HF, or pressure overload-induced HF condition was reflected to the network model by assigning the state value of particular input nodes to ‘1’ with an intensity determined in its corresponding input setting (Figure 3) while sampling a given number of initial states. To examine the robustness of the simulation results, we have performed additional computer simulation experiments while randomly perturbing the network structure and parameter settings to confirm the robustness of the dynamic analysis. Of note, two forms of network model perturbations were given in a combined manner during the attractor landscape analysis: (i) we have randomly rewired indicated percentage of the total number of links comprising the cardiac signaling network and (ii) logic tables of the indicated percentage of all internal regulatory nodes in the network were randomly perturbed. Identification of control nodes To identify the control nodes for minimizing the basin sizes of attractors for apoptosis and hypertrophy under ischemic or pressure overload-induced HF condition, we adopted the previously developed algorithm which identifies a set of nodes needed to be regulated to control the network states to converge to a desired attractor (Cho et al., 2016). During the mathematical simulation, the state value of each candidate control node was pinned to either 0 or 1. We employed GA that is used for optimization and search problems by artificially mutating, evolving, and selecting the most promising chromosomes containing all possible solutions (Kim et al., 2013). The fitness function was adopted from Cho et al. (2016), and is as follows: Fitness={B3×(n−∑i=1nXi)2×2,ifB=1B3×(n−∑i=1nXi)2,otherwise where B is the basin sizes of attractors for a desired phenotype, n is the number of nodes in the cardiac signaling network, and Xi is either 1 or 0 when node i is included or not in the chromosome X, respectively (Cho et al., 2016). Specifically, the basin size B is defined as follows: the basin of attraction converging to the attractors for normal phenotype when the candidate control nodes were perturbed under ischemic or pressure overload-induced HF condition. Here, attractors where all the output nodes used as indicators of hypertrophy or apoptosis of cardiomyocytes maintain ‘OFF’ state were determined as a normal phenotype. To measure the basin size B, we performed the attractor analysis of the cardiac signaling network model using one million randomly sampled initial states until the state trajectory of an initial state reaches a point attractor with a single network state or a cyclic attractor with repetitive states. Then, the fraction of initial states that converge to the attractors for normal phenotype was determined as the basin size B. The node set with the highest fitness value (i.e. the node set that can maximize the basin size B under distinct HF conditions) in the chromosome is chosen as a set of control nodes (Cho et al., 2016). Supplementary material Supplementary material is available at Journal of Molecular Cell Biology online. Funding This work was supported by the National Research Foundation of Korea (NRF) grants funded by the Korea Government, the Ministry of Science and ICT (2017R1A2A1A17069642, 2015M3A9A7067220, and 2013M3A9A7046303). Conflict of interest none declared. References ACE Inhibitor Myocardial Infarction Collaborative Group . ( 1998 ). Indications for ACE inhibitors in the early treatment of acute myocardial infarction: systematic overview of individual data from 100000 patients in randomized trials . Circulation 97 , 2202 – 2212 . CrossRef Search ADS PubMed Aikawa , R. , Komuro , I. , Yamazaki , T. , et al. . ( 1997 ). Oxidative stress activates extracellular signal-regulated kinases through Src and Ras in cultured cardiac myocytes of neonatal rats . J. Clin. Invest. 100 , 1813 – 1821 . Google Scholar CrossRef Search ADS PubMed Appukuttan , A. , Kasseckert , S.A. , Micoogullari , M. , et al. . ( 2012 ). Type 10 adenylyl cyclase mediates mitochondrial Bax translocation and apoptosis of adult rat cardiomyocytes under simulated ischaemia/reperfusion . Cardiovasc. Res. 93 , 340 – 349 . Google Scholar CrossRef Search ADS PubMed Backs , J. , Worst , B.C. , Lehmann , L.H. , et al. . ( 2011 ). Selective repression of MEF2 activity by PKA-dependent proteolysis of HDAC4 . J. Cell Biol. 195 , 403 – 415 . Google Scholar CrossRef Search ADS PubMed Barabasi , A.L. , Gulbahce , N. , and Loscalzo , J. ( 2011 ). Network medicine: a network-based approach to human disease . Nat. Rev. Genet. 12 , 56 – 68 . Google Scholar CrossRef Search ADS PubMed Barabasi , A.L. , and Oltvai , Z.N. ( 2004 ). Network biology: understanding the cell’s functional organization . Nat. Rev. Genet. 5 , 101 – 113 . Google Scholar CrossRef Search ADS PubMed Bernardo , B.C. , Weeks , K.L. , Pretorius , L. , et al. . ( 2010 ). Molecular distinction between physiological and pathological cardiac hypertrophy: experimental findings and therapeutic strategies . Pharmacol. Ther. 128 , 191 – 227 . Google Scholar CrossRef Search ADS PubMed Blankesteijn , W.M. , van de Schans , V.A. , ter Horst , P. , et al. . ( 2008 ). The Wnt/frizzled/GSK-3β pathway: a novel therapeutic target for cardiac hypertrophy . Trends Pharmacol. Sci. 29 , 175 – 180 . Google Scholar CrossRef Search ADS PubMed Bravo , C.A. , Vatner , D.E. , Pachon , R. , et al. . ( 2016 ). A food and drug administration-approved antiviral agent that inhibits adenylyl cyclase type 5 protects the ischemic heart even when administered after reperfusion . J. Pharmacol. Exp. Ther. 357 , 331 – 336 . Google Scholar CrossRef Search ADS PubMed Bujak , M. , and Frangogiannis , N.G. ( 2007 ). The role of TGF-β signaling in myocardial infarction and cardiac remodeling . Cardiovasc. Res. 74 , 184 – 195 . Google Scholar CrossRef Search ADS PubMed Chen , Y. , Ke , Q. , Yang , Y. , et al. . ( 2003 ). Cardiomyocytes overexpressing TNF-α attract migration of embryonic stem cells via activation of p38 and c-Jun amino-terminal kinase . FASEB J. 17 , 2231 – 2239 . Google Scholar CrossRef Search ADS PubMed Cho , S.H. , Park , S.M. , Lee , H.S. , et al. . ( 2016 ). Attractor landscape analysis of colorectal tumorigenesis and its reversion . BMC Syst. Biol. 10 , 96 . Google Scholar CrossRef Search ADS PubMed Choi , M. , Shi , J. , Jung , S.H. , et al. . ( 2012 ). Attractor landscape analysis reveals feedback loops in the p53 network that control the cellular response to DNA damage . Sci. Signal. 5 , ra83 . Google Scholar CrossRef Search ADS PubMed Colella , M. , Grisan , F. , Robert , V. , et al. . ( 2008 ). Ca2+ oscillation frequency decoding in cardiac cell hypertrophy: role of calcineurin/NFAT as Ca2+ signal integrators . Proc. Natl Acad. Sci. USA 105 , 2859 – 2864 . Google Scholar CrossRef Search ADS Condorelli , G. , Morisco , C. , Latronico , M.V. , et al. . ( 2002 ). TNF-α signal transduction in rat neonatal cardiac myocytes: definition of pathways generating from the TNF-α receptor . FASEB J. 16 , 1732 – 1737 . Google Scholar CrossRef Search ADS PubMed Cuenca , J. , Martin-Sanz , P. , Alvarez-Barrientos , A.M. , et al. . ( 2006 ). Infiltration of inflammatory cells plays an important role in matrix metalloproteinase expression and activation in the heart during sepsis . Am. J. Pathol. 169 , 1567 – 1576 . Google Scholar CrossRef Search ADS PubMed De Keulenaer , G.W. , and Brutsaert , D.L. ( 2007 ). Systolic and diastolic heart failure: different phenotypes of the same disease? Eur. J. Heart Fail . 9 , 136 – 143 . Google Scholar CrossRef Search ADS PubMed Fischer , P. , and Hilfiker-Kleiner , D. ( 2008 ). Role of gp130-mediated signalling pathways in the heart and its impact on potential therapeutic aspects . Br. J. Pharmacol. 153(Suppl 1) , S414 – S427 . Google Scholar PubMed Florea , V.G. , and Cohn , J.N. ( 2014 ). The autonomic nervous system and heart failure . Circ. Res. 114 , 1815 – 1826 . Google Scholar CrossRef Search ADS PubMed Francis , S.H. ( 2010 ). The role of cGMP-dependent protein kinase in controlling cardiomyocyte cGMP . Circ. Res. 107 , 1164 – 1166 . Google Scholar CrossRef Search ADS PubMed Fu , Y.C. , Chi , C.S. , Yin , S.C. , et al. . ( 2004 ). Norepinephrine induces apoptosis in neonatal rat cardiomyocytes through a reactive oxygen species-TNFα-caspase signaling pathway . Cardiovasc. Res. 62 , 558 – 567 . Google Scholar CrossRef Search ADS PubMed Fumia , H.F. , and Martins , M.L. ( 2013 ). Boolean network model for cancer pathways: predicting carcinogenesis and targeted therapy outcomes . PLoS One 8 , e69008 . Google Scholar CrossRef Search ADS PubMed Goldberg , D.E. ( 1989 ). Genetic Algorithms in Search, Optimization, and Machine Learning (1st edn) . Reading, MA : Addison-Wesley Professional . Gonzalez , A. , Fortuno , M.A. , Querejeta , R. , et al. . ( 2003 ). Cardiomyocyte apoptosis in hypertensive cardiomyopathy . Cardiovasc. Res. 59 , 549 – 562 . Google Scholar CrossRef Search ADS PubMed Gourdie , R.G. , Dimmeler , S. , and Kohl , P. ( 2016 ). Novel therapeutic strategies targeting fibroblasts and fibrosis in heart disease . Nat. Rev. Drug Discov. 15 , 620 – 638 . Google Scholar CrossRef Search ADS PubMed Guo , X. , and Wang , X.F. ( 2009 ). Signaling cross-talk between TGF-β/BMP and other pathways . Cell Res. 19 , 71 – 88 . Google Scholar CrossRef Search ADS PubMed Handler , D. , Meixner , K. , Pizka , M. , et al. . ( 2013 ). The genetic makeup of the Drosophila piRNA pathway . Mol. Cell 50 , 762 – 777 . Google Scholar CrossRef Search ADS PubMed Heineke , J. , and Molkentin , J.D. ( 2006 ). Regulation of cardiac hypertrophy by intracellular signalling pathways . Nat. Rev. Mol. Cell Biol. 7 , 589 – 600 . Google Scholar CrossRef Search ADS PubMed Helikar , T. , Konvalina , J. , Heidel , J. , et al. . ( 2008 ). Emergent decision-making in biological signal transduction networks . Proc. Natl Acad. Sci. USA 105 , 1913 – 1918 . Google Scholar CrossRef Search ADS Hilfiker-Klemer , D. , Landmesser , U. , and Drexler , H. ( 2006 ). Molecular mechanisms in heart failure—focus on cardiac hypertrophy, inflammation, angiogenesis, and apoptosis . J. Am. Coll. Cardiol. 48 , A56 – A66 . Google Scholar CrossRef Search ADS Huang , S. , Eichler , G. , Bar-Yam , Y. , et al. . ( 2005 ). Cell fates as high-dimensional attractor states of a complex gene regulatory network . Phys. Rev. Lett. 94 , 128701 . Google Scholar CrossRef Search ADS PubMed Huang , S. , Ernberg , I. , and Kauffman , S. ( 2009 ). Cancer attractors: a systems view of tumors from a gene network dynamics and developmental perspective . Semin. Cell Dev. Biol. 20 , 869 – 876 . Google Scholar CrossRef Search ADS PubMed Iwamoto , R. , Yamazaki , S. , Asakura , M. , et al. . ( 2003 ). Heparin-binding EGF-like growth factor and ErbB signaling is essential for heart function . Proc. Natl Acad. Sci. USA 100 , 3221 – 3226 . Google Scholar CrossRef Search ADS Jackson , G. , Gibbs , C.R. , Davies , M.K. , et al. . ( 2000 ). ABC of heart failure: Pathophysiology . BMJ 320 , 167 – 170 . Google Scholar CrossRef Search ADS PubMed Kanehisa , M. , and Goto , S. ( 2000 ). KEGG: kyoto encyclopedia of genes and genomes . Nucleic Acids Res. 28 , 27 – 30 . Google Scholar CrossRef Search ADS PubMed Kang , J.H. , Lee , H.S. , Kang , Y.W. , et al. . ( 2016 ). Systems biological approaches to the cardiac signaling network . Brief. Bioinform. 17 , 419 – 428 . Google Scholar CrossRef Search ADS PubMed Kang , J.H. , Lee , H.S. , Park , D. , et al. . ( 2017 ). Context-independent essential regulatory interactions for apoptosis and hypertrophy in the cardiac signaling network . Sci. Rep. 7 , 34 . Google Scholar CrossRef Search ADS PubMed Kawasaki , K. , Watabe , T. , Sase , H. , et al. . ( 2008 ). Ras signaling directs endothelial specification of VEGFR2+ vascular progenitor cells . J. Cell Biol. 181 , 131 – 141 . Google Scholar CrossRef Search ADS PubMed Kehat , I. , and Molkentin , J.D. ( 2010 ). Molecular pathways underlying cardiac remodeling during pathophysiological stimulation . Circulation 122 , 2727 – 2735 . Google Scholar CrossRef Search ADS PubMed Kim , S. , Kim , J. , and Cho , K.H. ( 2007 ). Inferring gene regulatory networks from temporal expression profiles under time-delay and noise . Comput. Biol. Chem. 31 , 239 – 245 . Google Scholar CrossRef Search ADS PubMed Kim , J.R. , Kim , J. , Kwon , Y.K. , et al. . ( 2011 ). Reduction of complex signaling networks to a representative kernel . Sci. Signal. 4 , ra35 . Google Scholar CrossRef Search ADS PubMed Kim , J. , Park , S.M. , and Cho , K.H. ( 2013 ). Discovery of a kernel for controlling biomolecular regulatory networks . Sci. Rep. 3 , 2223 . Google Scholar CrossRef Search ADS PubMed Kim , J. , Vandamme , D. , Kim , J.R. , et al. . ( 2014 ). Robustness and evolvability of the human signaling network . PLoS Comput. Biol. 10 , e1003763 . Google Scholar CrossRef Search ADS PubMed Kitano , H. ( 2002 ). Computational systems biology . Nature 420 , 206 – 210 . Google Scholar CrossRef Search ADS PubMed Koitabashi , N. , Danner , T. , Zaiman , A.L. , et al. . ( 2011 ). Pivotal role of cardiomyocyte TGF-β signaling in the murine pathological response to sustained pressure overload . J. Clin. Invest. 121 , 2301 – 2312 . Google Scholar CrossRef Search ADS PubMed Kong , P. , Christia , P. , and Frangogiannis , N.G. ( 2014 ). The pathogenesis of cardiac fibrosis . Cell. Mol. Life Sci. 71 , 549 – 574 . Google Scholar CrossRef Search ADS PubMed Kuramochi , Y. , Guo , X. , and Sawyer , D.B. ( 2006 ). Neuregulin activates erbB2-dependent src/FAK signaling and cytoskeletal remodeling in isolated adult rat cardiac myocytes . J. Mol. Cell. Cardiol. 41 , 228 – 235 . Google Scholar CrossRef Search ADS PubMed Law , V. , Knox , C. , Djoumbou , Y. , et al. . ( 2014 ). DrugBank 4.0: shedding new light on drug metabolism . Nucleic Acids Res. 42 , D1091 – D1097 . Google Scholar CrossRef Search ADS PubMed Leask , A. ( 2015 ). Getting to the heart of the matter: new insights into cardiac fibrosis . Circ. Res. 116 , 1269 – 1276 . Google Scholar CrossRef Search ADS PubMed Lee , H.S. , Goh , M.J. , Kim , J. , et al. . ( 2015 ). A systems-biological study on the identification of safe and effective molecular targets for the reduction of ultraviolet B-induced skin pigmentation . Sci. Rep. 5 , 10305 . Google Scholar CrossRef Search ADS PubMed Li , C. ( 2017 ). Identifying the optimal anticancer targets from the landscape of a cancer-immunity interaction network . Phys. Chem. Chem. Phys. 19 , 7642 – 7651 . Google Scholar CrossRef Search ADS PubMed Li , C. , and Wang , J. ( 2013 ). Quantifying cell fate decisions for differentiation and reprogramming of a human stem cell network: landscape and biological paths . PLoS Comput. Biol. 9 , e1003165 . Google Scholar CrossRef Search ADS PubMed Li , C. , and Wang , J. ( 2014 a). Landscape and flux reveal a new global view and physical quantification of mammalian cell cycle . Proc. Natl Acad. Sci. USA 111 , 14130 – 14135 . Google Scholar CrossRef Search ADS Li , C. , and Wang , J. ( 2014 b). Quantifying the underlying landscape and paths of cancer . J. R. Soc. Interface 11 , 20140774 . Google Scholar CrossRef Search ADS PubMed Lim , J.Y. , Park , S.J. , Hwang , H.Y. , et al. . ( 2005 ). TGF-β1 induces cardiac hypertrophic responses via PKC-dependent ATF-2 activation . J. Mol. Cell. Cardiol. 39 , 627 – 636 . Google Scholar CrossRef Search ADS PubMed Ling , H. , Gray , C.B. , Zambon , A.C. , et al. . ( 2013 ). Ca2+/Calmodulin-dependent protein kinase II delta mediates myocardial ischemia/reperfusion injury through nuclear factor-κB . Circ. Res. 112 , 935 – 944 . Google Scholar CrossRef Search ADS PubMed Liu , S. , Li , Y. , Kim , S. , et al. . ( 2012 ). Phosphodiesterases coordinate cAMP propagation induced by two stimulatory G protein-coupled receptors in hearts . Proc. Natl Acad. Sci. USA 109 , 6578 – 6583 . Google Scholar CrossRef Search ADS Liu , Y.Y. , Slotine , J.J. , and Barabasi , A.L. ( 2011 ). Controllability of complex networks . Nature 473 , 167 – 173 . Google Scholar CrossRef Search ADS PubMed Loan Le , T.Y. , Mardini , M. , Howell , V.M. , et al. . ( 2012 ). Low-dose spironolactone prevents apoptosis repressor with caspase recruitment domain degradation during myocardial infarction . Hypertension 59 , 1164 – 1169 . Google Scholar CrossRef Search ADS PubMed Lu , J. , Zeng , H. , Liang , Z. , et al. . ( 2015 ). Network modelling reveals the mechanism underlying colitis-associated colon cancer and identifies novel combinatorial anti-cancer targets . Sci. Rep. 5 , 14739 . Google Scholar CrossRef Search ADS PubMed Macarthur , B.D. , Ma’ayan , A. , and Lemischka , I.R. ( 2009 ). Systems biology of stem cell fate and cellular reprogramming . Nat. Rev. Mol. Cell Biol. 10 , 672 – 681 . Google Scholar CrossRef Search ADS PubMed MacLellan , W.R. , and Schneider , M.D. ( 2000 ). Genetic dissection of cardiac growth control pathways . Annu. Rev. Physiol. 62 , 289 – 319 . Google Scholar CrossRef Search ADS PubMed Matsui , T. , and Rosenzweig , A. ( 2005 ). Convergent signal transduction pathways controlling cardiomyocyte survival and function: the role of PI 3-kinase and Akt . J. Mol. Cell. Cardiol. 38 , 63 – 71 . Google Scholar CrossRef Search ADS PubMed Murray , P.J. , Kang , J.W. , Mirams , G.R. , et al. . ( 2010 ). Modelling spatially regulated β-catenin dynamics and invasion in intestinal crypts . Biophys. J. 99 , 716 – 725 . Google Scholar CrossRef Search ADS PubMed Nadal-Ginard , B. , Kajstura , J. , Leri , A. , et al. . ( 2003 ). Myocyte death, growth, and regeneration in cardiac hypertrophy and failure . Circ. Res. 92 , 139 – 150 . Google Scholar CrossRef Search ADS PubMed Nakamura , T. , Mizuno , S. , Matsumoto , K. , et al. . ( 2000 ). Myocardial protection from ischemia/reperfusion injury by endogenous and exogenous HGF . J. Clin. Invest. 106 , 1511 – 1519 . Google Scholar CrossRef Search ADS PubMed Nakaoka , Y. , Nishida , K. , Fujio , Y. , et al. . ( 2003 ). Activation of gp130 transduces hypertrophic signal through interaction of scaffolding/docking protein Gab1 with tyrosine phosphatase SHP2 in cardiomyocytes . Circ. Res. 93 , 221 – 229 . Google Scholar CrossRef Search ADS PubMed Nakaoka , Y. , Nishida , K. , Narimatsu , M. , et al. . ( 2007 ). Gab family proteins are essential for postnatal maintenance of cardiac function via neuregulin-1/ErbB signaling . J. Clin. Invest. 117 , 1771 – 1781 . Google Scholar CrossRef Search ADS PubMed Narula , J. , Haider , N. , Virmani , R. , et al. . ( 1996 ). Apoptosis in myocytes in end-stage heart failure . N. Engl. J. Med. 335 , 1182 – 1189 . Google Scholar CrossRef Search ADS PubMed Pan , J. , Fukuda , K. , Saito , M. , et al. . ( 1999 ). Mechanical stretch activates the JAK/STAT pathway in rat cardiomyocytes . Circ. Res. 84 , 1127 – 1136 . Google Scholar CrossRef Search ADS PubMed Park , S.G. , Lee , T. , Kang , H.Y. , et al. . ( 2006 ). The influence of the signal dynamics of activated form of IKK on NF-κB and anti-apoptotic gene expressions: a systems biology approach . FEBS Lett. 580 , 822 – 830 . Google Scholar CrossRef Search ADS PubMed Pennanen , C. , Parra , V. , Lopez-Crisosto , C. , et al. . ( 2014 ). Mitochondrial fission is required for cardiomyocyte hypertrophy mediated by a Ca2+-calcineurin signaling pathway . J. Cell Sci. 127 , 2659 – 2671 . Google Scholar CrossRef Search ADS PubMed Phanstiel , D.H. , Brumbaugh , J. , Wenger , C.D. , et al. . ( 2011 ). Proteomic and phosphoproteomic comparison of human ES and iPS cells . Nat. Methods 8 , 821 – 827 . Google Scholar CrossRef Search ADS PubMed Ponikowski , P. , Voors , A.A. , Anker , S.D. , et al. . ( 2016 ). 2016 ESC Guidelines for the diagnosis and treatment of acute and chronic heart failure: the Task Force for the diagnosis and treatment of acute and chronic heart failure of the European Society of Cardiology (ESC) developed with the special contribution of the Heart Failure Association (HFA) of the ESC . Eur. Heart J. 37 , 2129 – 2200 . Google Scholar CrossRef Search ADS PubMed Ponnusamy , M. , Li , P.F. , and Wang , K. ( 2017 ). Understanding cardiomyocyte proliferation: an insight into cell cycle activity . Cell. Mol. Life Sci. 74 , 1019 – 1034 . Google Scholar CrossRef Search ADS PubMed Rehsia , N.S. , and Dhalla , N.S. ( 2010 ). Mechanisms of the beneficial effects of β-adrenoceptor antagonists in congestive heart failure . Exp. Clin. Cardiol. 15 , e86 – e95 . Google Scholar PubMed Reimand , J. , Arak , T. , Adler , P. , et al. . ( 2016 ). g:Profiler—a web server for functional interpretation of gene lists (2016 update) . Nucleic Acids Res. 44 , W83 – W89 . Google Scholar CrossRef Search ADS PubMed Remme , W.J. , and Swedberg , K. ( 2001 ). Guidelines for the diagnosis and treatment of chronic heart failure . Eur. Heart J. 22 , 1527 – 1560 . Google Scholar CrossRef Search ADS PubMed Roberge , S. , Roussel , J. , Andersson , D.C. , et al. . ( 2014 ). TNF-α-mediated caspase-8 activation induces ROS production and TRPM2 activation in adult ventricular myocytes . Cardiovasc. Res. 103 , 90 – 99 . Google Scholar CrossRef Search ADS PubMed Ronca , R. , Gualandi , L. , Crescini , E. , et al. . ( 2009 ). Fibroblast growth factor receptor-1 phosphorylation requirement for cardiomyocyte differentiation in murine embryonic stem cells . J. Cell. Mol. Med. 13 , 1489 – 1498 . Google Scholar CrossRef Search ADS PubMed Ryall , K.A. , Holland , D.O. , Delaney , K.A. , et al. . ( 2012 ). Network reconstruction and systems analysis of cardiac myocyte hypertrophy signaling . J. Biol. Chem. 287 , 42259 – 42268 . Google Scholar CrossRef Search ADS PubMed Saraste , A. , Pulkki , K. , Kallajoki , M. , et al. . ( 1997 ). Apoptosis in human acute myocardial infarction . Circulation 95 , 320 – 323 . Google Scholar CrossRef Search ADS PubMed Schaefer , C.F. , Anthony , K. , Krupa , S. , et al. . ( 2009 ). PID: the pathway interaction database . Nucleic Acids Res. 37 , D674 – D679 . Google Scholar CrossRef Search ADS PubMed Schaub , M.C. , and Hefti , M.A. ( 2007 ). The PGE2-Stat3 connection in cardiac hypertrophy . Cardiovasc. Res. 73 , 3 – 5 . Google Scholar CrossRef Search ADS PubMed Shin , S.Y. , Kim , T. , Lee , H.S. , et al. . ( 2014 ). The switching role of β-adrenergic receptor signalling in cell survival or death decision of cardiomyocytes . Nat. Commun. 5 , 5777 . Google Scholar CrossRef Search ADS PubMed Shin , S.Y. , Yang , H.W. , Kim , J.R. , et al. . ( 2011 ). A hidden incoherent switch regulates RCAN1 in the calcineurin-NFAT signaling network . J. Cell Sci. 124 , 82 – 90 . Google Scholar CrossRef Search ADS PubMed Song , H.K. , Hong , S.E. , Kim , T. , et al. . ( 2012 ). Deep RNA sequencing reveals novel cardiac transcriptomic signatures for physiological and pathological hypertrophy . PLoS One 7 , e35552 . Google Scholar CrossRef Search ADS PubMed Spiro , Z. , Kovacs , I.A. , and Csermely , P. ( 2008 ). Drug-therapy networks and the prediction of novel drug targets . J. Biol. 7 , 20 . Google Scholar CrossRef Search ADS PubMed Stewart , S. , MacIntyre , K. , Capewell , S. , et al. . ( 2003 ). Heart failure and the aging population: an increasing burden in the 21st century? Heart 89 , 49 – 53 . Google Scholar CrossRef Search ADS PubMed Thygesen , K. , Alpert , J.S. , Jaffe , A.S. , et al. . ( 2012 ). Third universal definition of myocardial infarction . Circulation 126 , 2020 – 2035 . Google Scholar CrossRef Search ADS PubMed Travers , J.G. , Kamal , F.A. , Robbins , J. , et al. . ( 2016 ). Cardiac fibrosis: the fibroblast awakens . Circ. Res. 118 , 1021 – 1040 . Google Scholar CrossRef Search ADS PubMed Ucar , A. , Gupta , S.K. , Fiedler , J. , et al. . ( 2012 ). The miRNA-212/132 family regulates both cardiac hypertrophy and cardiomyocyte autophagy . Nat. Commun. 3 , 1078 . Google Scholar CrossRef Search ADS PubMed Vega , R.B. , Harrison , B.C. , Meadows , E. , et al. . ( 2004 ). Protein kinases C and D mediate agonist-dependent cardiac hypertrophy through nuclear export of histone deacetylase 5 . Mol. Cell. Biol. 24 , 8374 – 8385 . Google Scholar CrossRef Search ADS PubMed Wang , Y. ( 2007 ). Mitogen-activated protein kinases in heart development and diseases . Circulation 116 , 1413 – 1423 . Google Scholar CrossRef Search ADS PubMed Wang , J. , Li , C. , and Wang , E. ( 2010 ). Potential and flux landscapes quantify the stability and robustness of budding yeast cell cycle network . Proc. Natl Acad. Sci. USA 107 , 8195 – 8200 . Google Scholar CrossRef Search ADS Wang , W. , Xu , S. , Yin , M. , et al. . ( 2015 ). Essential roles of Gab1 tyrosine phosphorylation in growth factor-mediated signaling and angiogenesis . Int. J. Cardiol. 181 , 180 – 184 . Google Scholar CrossRef Search ADS PubMed Watts , D.J. , and Strogatz , S.H. ( 1998 ). Collective dynamics of ‘small-world’ networks . Nature 393 , 440 – 442 . Google Scholar CrossRef Search ADS PubMed Yancy , C.W. , Jessup , M. , Bozkurt , B. , et al. . ( 2013 ). 2013 ACCF/AHA guideline for the management of heart failure: a report of the American College of Cardiology Foundation/American Heart Association Task Force on practice guidelines . Circulation 128 , e240 – e327 . Google Scholar CrossRef Search ADS PubMed Yang , S. , Hu , S. , Hsieh , Y.C. , et al. . ( 2006 ). Mechanism of IL-6-mediated cardiac dysfunction following trauma-hemorrhage . J. Mol. Cell. Cardiol. 40 , 570 – 579 . Google Scholar CrossRef Search ADS PubMed Yang , J.H. , and Saucerman , J.J. ( 2011 ). Computational models reduce complexity and accelerate insight into cardiac signaling networks . Circ. Res. 108 , 85 – 97 . Google Scholar CrossRef Search ADS PubMed Zarubin , T. , and Han , J. ( 2005 ). Activation and signaling of the p38 MAP kinase pathway . Cell Res. 15 , 11 – 18 . Google Scholar CrossRef Search ADS PubMed Zhao , L. , Cheng , G. , Jin , R. , et al. . ( 2016 ). Deletion of Interleukin-6 attenuates pressure overload-induced left ventricular hypertrophy and dysfunction . Circ. Res. 118 , 1918 – 1929 . Google Scholar CrossRef Search ADS PubMed © The Author(s) (2018). Published by Oxford University Press on behalf of Journal of Molecular Cell Biology, IBCB, SIBS, CAS. All rights reserved. This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/about_us/legal/notices)

Journal

Journal of Molecular Cell BiologyOxford University Press

Published: Apr 24, 2018

There are no references for this article.

You’re reading a free preview. Subscribe to read the entire article.


DeepDyve is your
personal research library

It’s your single place to instantly
discover and read the research
that matters to you.

Enjoy affordable access to
over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month

Explore the DeepDyve Library

Search

Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly

Organize

Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place.

Access

Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals.

Your journals are on DeepDyve

Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more.

All the latest content is available, no embargo periods.

See the journals in your area

DeepDyve

Freelancer

DeepDyve

Pro

Price

FREE

$49/month
$360/year

Save searches from
Google Scholar,
PubMed

Create lists to
organize your research

Export lists, citations

Read DeepDyve articles

Abstract access only

Unlimited access to over
18 million full-text articles

Print

20 pages / month

PDF Discount

20% off