Context-independent essential regulatory interactions for apoptosis and hypertrophy in the cardiac signaling network

Apoptosis and hypertrophy of cardiomyocytes are the primary causes of heart failure and are known to be regulated by complex interactions in the underlying intracellular signaling network. Previous experimental studies were successful in identifying some key signaling components, but most of the findings were confined to particular experimental conditions corresponding to specific cellular contexts. A question then arises as to whether there might be essential regulatory interactions that prevail across diverse cellular contexts. To address this question, we have constructed a large-scale cardiac signaling network by integrating previous experimental results and developed a mathematical model using normalized ordinary differential equations. Specific cellular contexts were reflected to different kinetic parameters sampled from random distributions. Through extensive computer simulations with various parameter distributions, we revealed the five most essential context-independent regulatory interactions (between: (1) αAR and Gαq, (2) IP3 and calcium, (3) epac and CaMK, (4) JNK and NFAT, and (5) p38 and NFAT) for hypertrophy and apoptosis that were consistently found over all our perturbation analyses. These essential interactions are expected to be the most promising therapeutic targets across a broad spectrum of individual conditions of heart failure patients.


II. Supplementary Tables
Note that the mark ' → ' represents influence that is either activation or inhibition. Effect of catecholamine on the distribution of the 16 network components evaluated using each function type separately or using different types of response functions in a combined manner.

Supplementray
* Results of the mathematical simulation analysis performed using the combination of the four different response functions Supplementray Table S6. Results of the one-distribution perturbation analyses in the cardiac signaling network.
1) Result of one-distribution perturbation analysis for apoptosis or hypertrophy when each response function was separately applied Results of one-distribution perturbation analysis for apoptosis and hypertrophy when the threshold for determining the marginal distributions was set to top 5% or 20%. The effect is represented as the ratio between the degree of appearance of phenotypes in the one-distribution perturbation analysis and that in the control distributions. Parameters of which the marginal distributions significantly (p<0.05) increased apoptosis or hypertrophy in all response function types are shown. The mathematical analysis was all repeated for 10 times using different random seeds of 1 million parameter sets for each case. P-values were determined by comparison with the control distributions using Student's t test. See Supplementary Data Sets for full data. 33 2) Result of reverse one-distribution perturbation analysis for apoptosis or hypertrophy when each response function was separately applied Results of reverse one-distribution perturbation analysis for apoptosis and hypertrophy when the threshold for determining the marginal distributions was set to top 5% or 20%. The effect is represented as the ratio between the degree of appearance of phenotypes in the onedistribution perturbation analysis and that in the control distributions. Parameters of which the marginal distributions significantly (p<0.05) decreased apoptosis or hypertrophy in all response function types are shown. The mathematical analysis was all repeated for 10 times using different random seeds of 1 million parameter sets for each case. P-values were determined by comparison with the control distributions using Student's t test. See Supplementary Data Sets for full data.
3) Result of one-distribution perturbation analysis when four different response function types were applied in a combined manner Results of one-distribution perturbation analysis for apoptosis and hypertrophy when four different response function types were applied in a combined manner. The effect is represented as the ratio between the degree of appearance of phenotypes in one-distribution perturbation analysis and that in the control distributions. Threshold for determining the marginal distributions was set to top 5%, 10%, or 20%. Parameters of which the marginal distributions significantly (p<0.05) increased apoptosis or hypertrophy are shown. The mathematical analysis was all repeated for 10 times using different random seeds of 1 million parameter sets for each case. P-values were determined by comparison with the control distributions using Student's t test. See Supplementary Data Sets for full data. Results of two-distribution perturbation analysis for apoptosis and hypertrophy when four different response function types were applied in a combined manner. The synergistic effect was calculated as the difference between the effect of simultaneous perturbation of marginal distributions of two parameters on the phenotype and the sum of that obtained from perturbing either individual marginal distribution. Higher values indicate stronger synergistic effect. Parameter pairs exhibiting synergistic effect for apoptosis or hypertrophy with significance (p<0.05) are shown. The mathematical analysis was all repeated for 10 times using different random seeds of 1, 10, or 100 million parameter sets for each case. P-values were determined using Student's t test. See Supplementary Data Sets for full data.      Synergistic effect of GAB1 and SHP2 pm rasgap12 sum rasgap =pm rasgap1 × pm rasgap1 + pm rasgap2 × pm rasgap3 + pm rasgap3 × pm rasgap1 pm synrasgp123 = pm synrasgap ×pm synrasgap pm rasgap12 =(pm rasgap1 × pm rasgap2 / sum rasgap )×(pm synrasgap -pm synrasgap123 ) pm rasgap23 =(pm rasgap2 × pm rasgap3 / sum rasgap )×(pm synrasgap -pm synrasgap123 ) pm rasgap31 =(pm rasgap3 × pm rasgap1 / sum rasgap )×(pm synrasgap -pm synrasgap123 ) Pulsation is the ability to generate pulse at the moment of input signal is given or is removed.

Supplementray
Consistency measures similarity between the values of Z at each time point between those evaluated at different parameter sets.

II. Normalized equation modeling of the 16 network motifs
Differential equations of each model are made according to the principles described in METHODS (differential equations for all models are provided in Table S8). Each model has two parameters: p1 represents synergistic effect of two links while p2 represents dominant effect of one link over that of the other. p2 means dominant effect from Z to X over Input to X in negative feedback and positive feedback motifs while p2 means dominant effect from Y to Z over X to Z in coherent feedforward and incoherent feedforward motifs. Parameter p1 and p2 were sampled from uniform distribution between 0 and 1 with constraint that the sum of p1 and p2 should be less than or equal to 1. Total 1,000 parameter sets were generated. For each model and for each form of the response function, for each parameter set, numerical simulation was performed using ode15s function in the MATLAB R2009a (i.e. total 64,000 (16×4×1,000) simulations were conducted). Input signal is given from t=10 to t=100.
Dynamical features are evaluated for every simulation.

III. Observed dynamical characteristics in the 16 network motifs
Adaptation was mainly observed in negative feedback (motif 1, 2, 3, and 4) and incoherent feedforward motifs (motif 13, 14, 15, and 16). The result was consistent regardless of the form of the response function. Memory was only observed in positive feedback motifs (motif 5, 6, 7, and 8) except for motif 5 with Hill, motif 7 with Hill, and motif 8 with Hill.
Consistency was observed in all coherent feedforward motifs (motif 9, 10, 11, and 12) regardless of the form of the response function. Pulsation was observed in incoherent feedforward motifs (mainly with the function 'Hill'). From the results, we can conclude that the specific dynamical features are associated with the specific network motifs. These features may be determined by network structure of the motifs rather than the detailed arrangement of the links inside the motifs or form of the response functions. Therefore, dynamics of the network was successfully represented by normalized equation modeling. In Table S9, we summarized the association of network motifs and dynamical features. These results were consistent with previous studies (Table S10) (1, 2).

IV. Calculation of marginal distributions of the parameters
Then, which link is associated with the dynamical features? For example, which link is related with the property of adaption in type 1 negative feedback motif? To investigate the question, we observed the shape of the marginal distributions for the parameters in four such associations: association 1 -adaptation with negative feedback motifs, association 2adaptation with incoherent feedforward motifs, association 3 -memory with positive feedback motifs, association 4 -pulsation with incoherent feedforward motifs. Marginal distribution was calculated for the parameters that show upper 10% of the dynamics of interest for each association (Fig. S3A).

IV-1. Association 1: adaptation in negative feedback motifs
High p1 and low p2 were associated with adaptation in all types of negative feedback motifs. That is, the higher the synergistic effect of Input and Z, the lower the lower the dominant effect of Z to X in inhibition of X, adaptation can be observed more frequently.
Therefore, coordination of two links in regulation of X may be important to generate the dynamical feature, "adaptation". The distributions were similar regardless of the form of the response functions.

IV-2. Association 2: adaptation in incoherent feedforward motifs
For adaptation, different distributions of parameters were observed depending on the type of incoherent feedforward motifs (in motif 13 and 16, high p1 and low p2 was associated with adaptation while low p1 and high p2 was associated with adaptation in motif 14 and motif 15). However, the marginal distributions were not different according to the form of response functions.

IV-3. Association 3: memory in positive feedback motifs
Calculation of marginal distributions for type 1 positive feedback with Hill, type 3 positive feedback with Hill, and type 4 positive feedback with Hill was impossible because of zero or small number of cases. Except for these, low p1 and high p2 are associated with memory.
That is, the lower the synergistic effect of Input and Z, the higher the dominant effect of Z to X in activation of X, adaptation can be observed more frequently. Therefore, positive feedback link may take a role to generate the dynamical feature, "memory".

IV-4. Association 4: pulsation in incoherent feedforward motifs
Low p1 and medium p2 were associated with pulsation in all types of incoherent feedforward motifs. That is, pulsation was observed more frequently when the synergistic effect of Input and Z was low and the effect of Y on Z (indirect effect for the activation of Z) was neither high nor low. The latter is an interesting point. Therefore, coordination of direct effect (i.e. X on Z) and indirect effect (i.e. Y on Z) may be required to generate the dynamical feature, "pulsation". The distributions were similar regardless of the form of the response functions.

V. Distribution perturbation analysis
Although distinct distributions of parameters were observed for specific dynamics in all previously investigated network-dynamics associations, it does not always imply that these distributions generate the dynamics. To further investigate this issue, we performed "distribution perturbation analysis". After the one or two parameters were sampled from marginal distribution while the others were sampled from uniform distribution, numerical simulation was performed and change of phenotypes (or dynamics) was evaluated. If the results showed certain dynamics more frequently, it will suggest that the distribution generates the dynamics. For one distribution perturbation analysis, we sampled the value of one parameter (p1 or p2) from marginal distribution whereas the other (p2 or p1) was sampled from uniform distribution. For two distribution perturbation analysis, we sampled both parameters (p1 and p2) from marginal distribution. One constraint required to be satisfied during the process of random sampling is that the sum of p1 and p2 is less 1. The results of distribution perturbation analyses were demonstrated Fig. S4 and Table S11.

V-1. Association 1: adaptation in negative feedback motifs
For all types of negative feedback motifs and all forms of response functions, perturbation of any one distribution and two distributions generated more degree of adaptation. The perturbation effect of the two distributions was always larger than that of one marginal distribution. This results implies that high p1 and low p2 (i.e., strong synergistic effect with dominant role of input) actively generates the dynamical feature, adaptation.

V-2. Association 2: adaptation in incoherent feedforward motifs
Distribution perturbation has more impacts in type 1 and type 4 than in type 2 and type 3.
That is, marginal distributions of p1 and p2 in motif 14 and motif 15 do not have significant causal relationship with adaptation.

V-3. Association 3: memory in positive feedback motifs
For all types of positive feedback motifs and all forms of response functions, distribution perturbation of p2 generates memory more effectively. In addition, differences are very slight between the results of distribution perturbtion of p2 and those of distribution perturbation of both p1 and p2. These results imply that high p2 have more significant causal relationship with memory than low p1.

V-4. Association 4: pulsation in incoherent feedforward motifs
For all types and all forms of response functions, perturbation of any one distribution and two distribution generated more degree of adaptation. The perturbation effect of two distributions was always larger than that of one distribution. This results implies that low p1 and medium p2 (i.e., strong synergistic effect with dominant role of input) actively generates the dynamical feature, pulsation.

VI. Summary
Normalized equation modeling revealed well-known dynamical features of network motifs: negative feedback and incoherent feedforward motifs showed adaptation; positive feedback motifs showed memory; coherent feedforward motif showed consistency.
Furthermore, the conditions for the dynamic features can be identified through distribution perturbation analysis. Most of the results were consistent regardless of the form of the response function which strengthen the reliability of the developed method.

I. Construction of the EGFR network
Network structure is the simplified form of the previously constructed EGFR network (Fig.   S5A).

II. Normalized equation modeling of the EGFR network.
Differential equations of each model are made according to the principles described in METHODS (differential equations for all models are provided in Table S12). The model requires 14 parameters (i.e. pmsyngs, pmgs1, pmgs2, pmsynpitk, pmpitk1, pmpitk2, pmras, pmsynras, pmsyngab, pmgab1, pmgab2, pmsynrasgap, pmrasgap1, pmrasgap2) which are sampled from standard uniform distribution. The remaining parameter can be calculated from the sampled parameters. In this manner, 100,000 parameter sets were randomly generated. The detailed information of the parameters is provided in Table S13.

III. Measurement of the MEK inhibitor resistance
Numerical simulation was performed using ode15s function in the MATLAB R2009a.
MEK inhibitor is treated at t=100 (Fig. S5B). The resistance of the MEK inhibitor was calculated as the increased amount of AKT after MEK inhibitor is given.

IV. Marginal distribution of parameters
To find out which parameter is associated with MEK inhibitor resistance, we investigated the marginal distributions of all 14 parameters. Marginal distribution was calculated for the parameter sets that showed top 10% of the MEK inhibitor resistance. Among 14 parameters, six parameters (pmgs2, pmpitk1, pmsynras, pmras, pmsynrasgap, pmrasgap) represented non-uniform, coherent marginal distributions (Fig. S5C).
Among them, the importance of GAB1 regulation was supported by published experimental results (3).

VI. Summary
We investigated the underlying mechanism of MEK inhibitor resistance in the EGFR network using normalized equation modeling and one-distribution perturbation analysis. The regulation of Grb2Sos, Ras, and GAB1 were identified as essential regulatory processes to generate MEK inhibitor resistance, among which GAB1 regulation had been also identified as a key process in previously published results.
Characterization of IGF-1-induced activation of the JAK/STAT Pathway in Rat Cardiomyocytes.
The inhibitory effect of ANP on ET-1-stimulated IP3 production was abolished in cells overexpressing RGS4DN ( Figure 2A) 18443239 Gbg Beta and gamma subunit of G protein Gbg is generated simultaneously with activation of Gi or  5) GSK3beta phosphorylates conserved serines necessary for nuclear export, promotes nuclear exit.

7)
In the presence of Ang2 or PE, NFAT is up-regulated. This up-regulation was completely abolished in the presence of CsA or FK506, supporting the conclusion that AngII and PE activate NFAT through a calcineurin-dependent signal transduction pathway ( Figure 4H).
Here, biological meaning of 'NFATact' is the transcriptional activity of NFAT.
1) MEK1-ERK1/2 signaling augments NFAT transcriptional activity independent of calcineurin, independent of changes in NFAT nuclear localization, and independent of alterations in NFAT transactivation potential. In contrast, MEK1-ERK1/2 signaling enhances NFAT-dependent gene expression through an indirect mechanism involving induction of cardiac AP-1 activity, which functions as a necessary NFAT interacting partner.

5) GSK3B inhibits the DNA binding activity of NFATc
3) Overexpression of wild-type or constitutively active mutant c-Src, but not kinase inactive mutant c-Src, lead to increased tyrosine kinase activity in Shc immunoprecipitates.
2) Inhibitory feedback phosphorylation of SOS by ERK provides a mechanism for the inhibition of Ras signaling.
3) SHP2 that binds to Gab1 was reported to be a positive regulator of the MAPK pathway. Grb2&Gab1 growth factor receptor-bound protein 2(Grb2)-Grb2-associated binding protein 1(Gab1) complex