The switching role of β-adrenergic receptor signalling in cell survival or death decision of cardiomyocytes

How cell fate (survival or death) is determined and whether such determination depends on the strength of stimulation has remained unclear. In this study, we discover that the cell fate of cardiomyocytes switches from survival to death with the increase of β-adrenergic receptor (β-AR) stimulation. Mathematical simulations combined with biochemical experimentation of β-AR signalling pathways show that the gradual increment of isoproterenol (a non-selective β1/β2-AR agonist) induces the switching response of Bcl-2 expression from the initial increase followed by a decrease below its basal level. The ERK1/2 and ICER-mediated feed-forward loop is the hidden design principle underlying such cell fate switching characteristics. Moreover, we find that β1-blocker treatment increases the survival effect of β-AR stimuli through the regulation of Bcl-2 expression leading to the resistance to cell death, providing new insight into the mechanism of therapeutic effects. Our systems analysis further suggests a novel potential therapeutic strategy for heart disease.

Supplementary Fig. 5. The temporal profiles of cAMP-related apoptotic molecules in response to ISO stimulation. Cardiomyocytes were stimulated with 1 μM ISO, and the mRNA level of cAMPrelated apoptotic molecules were measured over time by qRT-PCR. The expression levels of Xiap, cIAP2, and Bim did not show any significant change in response to ISO. cIAP1 was transiently increased for the initial 12 hours and then decreased to its basal level. The data represent means ± SEM, n≥3 biological and technical replicates (independent culture preparations). *, p < 0.05; **, p < 0.01; ***, p < 0.001 versus control group; Student's t-test.
Supplementary Fig. 6. The effect of the duration of β-AR stimuli on the cellular decision for survival or death. To examine the effect of the duration of the stimuli on the Bcl-2 expression level and cell viability, experiments and simulation were performed. (a) A schematic diagram showing the protocol of investigating the effect of the duration of β-AR stimuli on the cellular decision for survival or death. Isolated adult cardiomyocytes were incubated with indicated concentrations of ISO for 1, 3, 6, 9 or 12 hours. At the end of stimulation, ISO was washed out. Cardiomyocytes were lysed for western blotting analysis or imaged for calculating the survival rate at 12 hours from the start of ISO treatment. (b-c) Representative Bcl-2 and α-tubulin (loading control) immunoblots and the semiquantified data showed that the effect of ISO on the cellular decision for survival or death at the low or high concentration was enhanced along with the washout time, respectively (the longer the stimulation, the stronger the effect). The data represent means ± SEM, n ≥ 4 biological and technical replicates (independent culture preparations). (d) Both simulation (d) and experimental (c) results showed the linear correlation between the duration of stimulation and the enhancement of the effect on the expression level of Bcl-2. On the other hand, the cell fate determination does not change by the duration of the stimuli. (e-f) The survival rate was assessed by live-cell imaging. Data represent means ± SEM, n = 3 biological and technical replicates (independent culture preparations), n=20 for simulation (d). *, p < 0.05; **, p<0.01; ***, p<0.001 compared to control group; Student's t-test. (a) Representative immunoblots for Bax, Bcl-2, and α-tubulin (loading control) expressed with 12 hours of incubation at the indicated concentrations of ISO. The Bcl-2 expression level obviously increased at the concentration range of 10 -10 -10 -9 M, but further decreased to below its basal level at the concentration range of 10 -7 -10 -6 M, but the Bax expression level was not changed significantly. (b) The ratio of Bax to Bcl-2 quantified from the immunoblots shown in (a).
Supplementary Fig. 8. PKA plays an important role in increasing the level of Bcl-2 expression at the low concentration range of ISO through ERK1/2 activities. Simulation results show that the PKA inhibitor (RP-cAMPS) significantly suppressed ERK1/2 activity (left panel) and the switching profile of Bcl-2 disappears after the PKA inhibition (right panel). The data represent mean ± SEM for the repetitive simulations (n = 20) over up to 20% random variation of parameter values.
Supplementary Fig. 9. The efficiency of siRNA transfection in adult cardiomyocytes. To evaluate the efficiency of siRNA transfection into adult cardiomyocytes, siGLO Transfection Indicators (Thermo Scientific Inc.) was used. 24 hours after transfection, cardiomyocytes were imaged using a confocal microscope. Most of the rod-shape and multi-nuclei cardiomyocytes contain siGLO (green) signal, which is enriched in the nuclear region (arrow heads). Scale bar = 100 μm. Supplementary Fig. 10. Validation of the role of ERK1/2-mediated feed-forward loop using a RSK inhibitor. Cardiomyocytes were stimulated with the indicated concentrations of ISO in the presence of BRD7389. Representative immunoblots and plots of semi-quantification showing Bcl-2 expression at the indicated ISO concentration determined in the presence or absence of BRD7389. Data represent means ± SEM, n = 4 biological and technical replicates (independent culture preparations). The control data were taken from Fig. 6i. *, p < 0.05; Student's t-test. Supplementary Fig. 11. Potential therapeutic targets were predicted from the sensitivity analysis of the kinetic parameters and the total concentrations. (A) Before performing the sensitivity analysis, we had screened all the kinetic parameters and identified 21 parameters that could be the potential therapeutic targets by surveying the previous experimental data (see Supplementary Table 1  The simulation results show that Bcl-2 level was increased during the initial period of treatment with a cAMP analogue (8-CPT-cAMP), but then it was decreased for a longer treatment. The data represent mean ± SEM for the repetitive simulations (n = 20) over up to 30% random variation of parameter values.

II. Supplementary Tables
Supplementary Table 1   We assumed that the volume of cardiomyocytes is 26.1 pL (where the radius of cardiomyocyte is 8 μm and the length is 130 μm) and the protein amount in one cardiomyocyte is 2.8 ng (where the number of cells in a 35-mm culture dish is 95,000 and the total amount of protein from one 35mm culture dish is 270 μg). Based on these parameter values, the total concentration of the protein was estimated from the RNA-seq data 1 using the following equation: where RPKM is the Reads Per Kilobase per Million mapped reads, Conc is the concentration of a protein, M is the total number of isoforms of the protein, N is the number of genes, MW is the molecular weight of the protein, A is the protein amount of the cardiomyocyte (2.8 ng), and V is the volume of the cardiomyocyte (26.1 pL). We also assumed that the amount of a protein is linearly correlated with the amount of mRNA 2 .

Supplementary Note 1. Mathematical model
The computational model of β-AR signalling network is composed of four major modules: cAMP-PKA signalling module, central feedback regulation module, extracellular signalregulated kinases 1 and 2 (ERK1/2) signalling module and Ca 2+ regulatory module. The cAMP-PKA signalling module includes two types of β-AR (β 1 and β 2 ) which sequentially activate G s protein and adenylyl cyclase (AC) 36,37 . The active AC promotes the accumulation of cytosolic cAMP and subsequently leads to PKA activation 38 . The active PKA phosphorylates β 1 -and β 2 -AR simultaneously; however, the functional consequences of these two receptors are quite different 37,39,40 : The phosphorylated β 1 -AR is internalized and thus desensitized to the stimulation of agonist, while the phosphorylated β 2 -AR is coupled from G s to G i 41 . Similar to PKA, GRK2 also phosphorylates both types of β-ARs but the difference is that GRK2 phosphorylates only the ligand bound receptors 42,43 . Thus, the beta-AR signalling pathway is regulated by two feedback loops that are mediated by PKA and GRK2. In this study, for the simplicity of mathematical modeling, we explicitly included only the feedback loops of PKA in the model, since we assumed that the feedback mechanism mediated through PKA and GRK2 are functionally redundant in terms of the cell fate determination. Note, however, that the effect of GRK2 was still reflected in our model because the model parameters were fitted to the time course data, where the functions of GRK2 were already reflected.
The central feedback regulation module consists of CREB, ICER and PDE3. The active PKA is translocated to the nucleus and phosphorylates CREB, which leads to the transcription of inducible cAMP early repressor (ICER) and Bcl-2 44 . The induced ICER protein is further stabilized by cAMP-dependent signalling 45 . ICER protein contains DNA binding and leucine zipper domains but not the N-terminal transactivation domain, which makes it function as endogenous inhibitors of gene transcription driven by its cognates such as CREB, including its own expression by competing with CREB in binding CRE sequence 44 . In addition, ICER transcriptionally represses the expression of PDE3A, which subsequently contributes to the acceleration of the cAMP accumulation 44 . Thus, in modeling, ICER mediates two feedback loops: one is constructed by its own negative autoregulation and the other one is PDE3mediated feedback loop. Note that although a number of experimental results showed that PDE4 may play an essential role in cellular processes such as muscle contraction 46 , for the simplicity of mathematical modeling, we explicitly included only PDE3 in the model since PDE4 does not significantly affect the apoptosis of cardiomyocytes, even though its inhibition increases the cAMP level 47 . However, the functional effect of PDE4 was implicitly reflected in our model by fitting the model parameters to the time-course data of cAMP and PKA measured in the presence of PDE4. Regarding the compartmentalization of cAMP signal, we also simplified the model by assuming that the cAMP signal in cytosol is homogeneous, since the coarse-grained network that simplifies intermediate signalling steps would not significantly affect the steadystate properties of the system 48 .
As mentioned above, PKA switches the coupling of β 2 -AR for G s to G i . This switching releases G i _α_GTP from G i _β/γ subunit. G i _α_GTP inhibits the AC activity while promoting the binding of SOS and Grb2 via β-arrestins, which subsequently activates the ERK1/2 signalling module 41 . ERK1/2 regulates the central feedback module in two ways: the active ERK1/2 destabilizes the ICER protein 45 and activates CREB via the phosphorylation of RSK 49,50 . The Ca 2+ regulatory module is quite complex because it includes many channels and pumps that regulate the inward and outward Ca 2+ flux through the plasma membrane and intracellular Ca 2+ store 51 . However, we do not intend to investigate all of the detailed mechanism in this study. Instead, we construct a basic and minimal module. As a result, this module includes the Ca 2+ regulatory machinery, calcineurin and Ca 2+ /calmodulin-dependent protein kinase II, which regulates the β-AR signalling network through feedforward and feedback regulations.
The model was developed using ordinary differential equations (ODEs) on the basis of Michaelis-Mention-type functions and mass action law. In particular, the reaction steps of β 1and β 2 -ARs in the model are completed in a closed system, where all components are linked by reversible reactions and thereby the reactions consist of a closed-loop. Thus, for these reaction steps we applied detailed balance principle as previously reported 52,53 . In other words, the product of rate constants in one direction around the closed loop is equal to the product of rate constants in the opposite direction around the loop. As explained in the above, the simplification of some signalling components might reduce the accuracy of simulated temporal dynamics. However, that would not much significantly affect the long-term response in terms of cell fate determination 48 .

Supplementary Note 3. Sensitivity analysis
We carried out the global sensitivity analysis by perturbing the kinetic parameter values of all the regulatory links (39 links in total) in the range of 10-fold, and examined their influences on the Bcl-2 switching response profile. The results of this global sensitivity analysis showed that the parameters (vm18, kc19a, ki39, kc48, kc21b, vm23, vm25, vm27, vm34 and ki34) associated with eight essential regulatory links have significant effects on shaping the Bcl-2 switching response profile ( Supplementary Fig. 12a), consistent with the results obtained from our coarsegraining network analysis (Fig. 5). In the next, we further investigated the hidden mechanism underlying the cellular decision for survival or death response of Bcl-2. For this purpose, the parameter sets were clustered into two groups leading to either 'survival-only' or 'death-only' response which was defined on the basis of Bcl-2 response profile. By comparing the clustered parameter values, we found that most of the parameters (e.g., kc19b, kc19a, kc28b, ka39b, vm23, ki23, ki29 and ki34) representing the eight essential regulatory links (Fig. 5b) showed statistically significant differences, implying that they play a central role in determining the Bcl-2 switching response property (Supplementary Fig. 12b).