Skip to main content

Advertisement

We’d like to understand how you use our websites in order to improve them. Register your interest.
Download PDF

Modeling cardiac β-adrenergic signaling with normalized-Hill differential equations: comparison with a biochemical model

BMC Systems Biology volume 4, Article number: 157 (2010) Cite this article

Abstract

Background

New approaches are needed for large-scale predictive modeling of cellular signaling networks. While mass action and enzyme kinetic approaches require extensive biochemical data, current logic-based approaches are used primarily for qualitative predictions and have lacked direct quantitative comparison with biochemical models.

Results

We developed a logic-based differential equation modeling approach for cell signaling networks based on normalized Hill activation/inhibition functions controlled by logical AND and OR operators to characterize signaling crosstalk. Using this approach, we modeled the cardiac β1-adrenergic signaling network, including 36 reactions and 25 species. Direct comparison of this model to an extensively characterized and validated biochemical model of the same network revealed that the new model gave reasonably accurate predictions of key network properties, even with default parameters. Normalized Hill functions improved quantitative predictions of global functional relationships compared with prior logic-based approaches. Comprehensive sensitivity analysis revealed the significant role of PKA negative feedback on upstream signaling and the importance of phosphodiesterases as key negative regulators of the network. The model was then extended to incorporate recently identified protein interaction data involving integrin-mediated mechanotransduction.

Conclusions

The normalized-Hill differential equation modeling approach allows quantitative prediction of network functional relationships and dynamics, even in systems with limited biochemical data.

Background

The β-adrenergic signaling pathway plays a key role in the regulation of normal heart function and the development of heart failure [1–5]. Systems analysis of β-adrenergic signaling in the heart may provide important new insights into the mechanisms of heart failure and reveal new therapeutic targets. Previous mathematical models of cardiac β-adrenergic signaling have characterized how biochemical mechanisms of this pathway determine its coordinated regulation of cell contractility in health and disease [6–8]. However, this work relied on extensive biochemical data from the literature that may not be available for more recently discovered pathways. Therefore, more scalable modeling approaches are needed.

As an alternative to generating biochemically detailed kinetic models, several modeling approaches that are more closely based on network topology have been developed including Boolean modeling [9], fuzzy logic modeling [10] and extreme pathways analysis [11]. These approaches require few or no parameters and facilitate large-scale analysis of systems properties, such as feedback loops and feasible solution spaces. But these approaches have a variety of limitations. While extreme pathways analysis predicts the entire feasible steady-state solution space of a network, its ability to predict dynamic time-courses for given experiments is limited [12]. Simulations from discrete-level models (e.g. Boolean) can be difficult to interpret due to sensitivity of model predictions to temporal updating schemes [13], assignment of discrete activity-levels to continuous-valued variables like concentration [14], and the limited ability to describe realistic timescales [15]. The tradeoffs inherent in many of these logic-based modeling approaches has recently been reviewed [16]. In addition, these modeling approaches are generally not compatible with the wealth of systems analysis tools for differential equations from control theory and dynamical systems. Piecewise-linear differential equation models overcome some of these limitations by making both species values and time continuous, but steady-state species activities are still binary [9, 15, 17]. Others have modeled signaling networks with continuous approximations of Boolean functions [18] that are implemented to minimize steady-state differences between Boolean and continuous models.

To address these limitations, we developed a normalized-Hill differential equation modeling approach that combines advantages of both biochemical and Boolean models. This approach uses normalized Hill functions and logical AND and OR operators to describe network crosstalk. We used this approach to model the cardiac β-adrenergic signaling pathway and performed a direct comparison with a previously validated biochemical model of the same network [6, 7]. We then used this model to gain insight into the roles of feedback and feed-forward loops in the β-adrenergic pathway and examined potential crosstalk with integrin-mediated mechanotransduction. The analysis presented here demonstrates that the normalized-Hill differential equation modeling approach can provide reasonably accurate predictions of signaling properties, even when little parameter data is available.

Results

Toy signaling network

For demonstration, we created a toy signaling network using our normalized-Hill differential equation approach. This simple network consists of two input ligands ("A" and "B") that activate receptors "C" and "D", respectively. A positive feedback loop exists between "C" and "E" that is inhibited when "D" is activated (see Figure 1A). The state variables represent the "fractional activation" of the signaling species, which is normalized to the maximal possible activity. Fractional activation varies continuously with time and can take on any value between 0 and 1, inclusive. For example, fractional activation for a substrate that is active only when phosphorylated is equivalent to the ratio of phosphorylated to total protein.

Figure 1
figure1

Normalized-Hill toy network model. A) Schematic of the 5-species toy network, including two inputs, an AND reaction, and a positive feedback loop. B) Characteristics of sample normalized-Hill functions (n = 4 for both curves). C) Simulated signaling dynamics in response to transient inputs A, B, and A + B. Output E remains active after removal of input A due to a bistable positive feedback loop, though this memory is erased with addition of B. Addition of both stimuli at the same time transiently activates output E. Parameters used are: τ = 1, W = 1, n = 1.4, and EC50 = 0.5.

Full size image

Model equations for this toy network are provided in "Methods" and Additional File 1, Supplemental Methods, though properties of our modeling approach are discussed here. Interactions between species are modeled using normalized Hill functions with 3 reaction parameters: the reaction weight "W", half-maximal effective concentration "EC50", and Hill coefficient "n". The reaction weight determines how much a given interaction activates or inhibits an output species and can take on any value between 0 and 1, inclusive. EC50 is the fractional activation of an input species required to induce half-maximal activation of an output species. Lastly, the Hill coefficient determines sensitivity to changes in inputs. The normalized Hill functions are constrained to f(0) = 0, f(1) = 1 and f(EC50) = 0.5. Additionally, species activities are controlled by the parameters τ and YMAX, which are the reaction time constant and species maximal fractional activation, respectively. While YMAX = 1 is typical, this value can be altered to reflect a change in protein expression relative to a reference condition. Typical default reaction and species parameter values are W = 1, EC50 = 0.5, n = 1.4, τ = 1, and YMAX = 1; choice of default values for EC50 and n are examined in more detail below. Crosstalk between species is modeled using continuous functions analogous to Boolean AND and OR operations (see Methods).

While default parameters can provide qualitatively reasonable results, the parameters involved in these Hill functions can be directly measured in cellular experiments to quantitatively refine model predictions. For example, in the toy model, WD, KB and nB could all be determined directly from a single steady-state concentration response where B is varied and D is measured. Likewise, τD could be determined experimentally by measuring the dynamic response of D in response to a step in B. If D is a kinase substrate, typical experiments could include quantitative Western blots, cellular immunofluorescence, or live-cell FRET biosensors similar to those used previously for measuring active PKA dynamics in cardiac myocytes [19].

Figure 1C shows a sample simulation of toy network signaling dynamics in response to a transient exposure to input "A", input "B", followed by both inputs simultaneously. Note that "C" and "E" retain full activity even after input "A" has been removed due to the presence of a bistable positive feedback loop. However, this memory is erased once input "B" is activated. When "A" and "B" are activated simultaneously, "E" is transiently activated but switches off once precursors "C" and "D" become highly active.

Because this modeling approach uses ordinary differential equations, well-established systems approaches, such as quantitative sensitivity analysis, can be readily applied to study network relationships. The sensitivity plot shown in Figure 2 quantifies the global functional relationships between all species in the system. Each column simulates an individual experiment in which the maximum activity of one species is perturbed (e.g. "AMAX") and the subsequent changes in all network species' activities are quantified. The normalized steady-state sensitivities of these species to parameter perturbations are quantified according to S = (ΔY/ΔP)(Po/Yo) (see Methods for details). Graded positive and negative sensitivities are represented with shades of red and green, respectively. For example, the "B" column indicates that increased fractional activation of "B" causes increased activity of itself and "D", decreased activity of "E" and (to a lesser extent) "C", and no change in the activity of "A". This type of quantitative sensitivity analysis is not possible with many other logic-based or topological approaches.

Figure 2
figure2

Sensitivity analysis of the toy network: Each column of the sensitivity matrix represents a numerical experiment in which one species was perturbed by -25%. The steady-state sensitivities of all model species were computed according to S = (ΔY/ΔP)(Po/Yo). The model predicts that E is positively regulated by A and C but negatively regulated by B and D.

Full size image

β-adrenergic signaling network

We next applied our normalized-Hill modeling approach to the cardiac β-adrenergic signaling network, the major pathway that regulates contractility, metabolism, and gene regulation of the heart [20]. Following stimulation by β agonists (e.g. norepinephrine or NE), the β-adrenergic receptor couples with G proteins, which subsequently activate adenylyl cyclase (AC). AC synthesizes the second messenger cAMP, which causes dissociation of the regulatory and catalytic subunits of protein kinase A (PKAR and PKAC, respectively). PKAC phosphorylates several substrates, including phospholamban (PLB), the ryanodine receptor (RyR), troponin I (TnI), and inhibitor-1 (Inhib1), resulting in enhanced cardiac contractility. Figure 3 shows a schematic of this model, which includes 25 species and 36 reactions. All interactions were implemented using the same default parameters used in the toy signaling network. Importantly, the network topology for this normalized-Hill model was designed to be very similar to an extensively characterized and validated biochemical model of this same network [7].

Figure 3
figure3

Schematic of normalized-Hill β-adrenergic signaling network. The normalized-Hill model consists of 25 species and 36 reactions. While developed using a normalized-Hill differential equation approach, the network topology closely mimics a well-validated biochemical model [6] to allow for direct comparisons.

Full size image

To test the normalized-hill β-adrenergic model, we simulated a simple experiment in which NE is added to the system and subsequently removed. The addition of NE indirectly activates G proteins (Gsα), cAMP, and PLB (Figure 4A). The response delay of cAMP and PLB to NE (as compared to Gsα) occurs because these species are farther downstream in the signaling cascade. An adaptive response to the sustained NE input is observed in these activation profiles, as seen experimentally for cAMP [21, 22], which may be attributed to the negative feedback loops consisting of PKAC and GRK phosphorylation of β1ARs (see below). Responses for individual species activities during simulations with varying levels of NE were also evaluated; a sample dose response is shown in Figure 4B, where increased levels of NE produced concomitant increases in the fractional activities of Gsα, cAMP, and PLB, as expected from experiments [23, 24].

Figure 4
figure4

Predicted dynamics and dose responses of the normalized-Hill β-adrenergic model. A) Predicted dynamics of G proteins (Gsα), cyclic AMP (cAMP), and phospholamban (PLB) phosphorylation in response to a transient norepinephrine (NE) exposure (NE = 1). B) Steady-state concentration responses of Gsα, cAMP and PLB to varying levels of NE.

Full size image

Role of feedback and forward loops in shaping network dynamics

Feedback and feed-forward loops are key network motifs that can enhance information processing by altering signaling dynamics or providing adaptation [25]. In this network, there are two negative feedback loops, both involving receptor phosphorylation by either PKA or G-protein receptor kinase (GRK), which desensitize the β-adrenergic receptor to ligand inputs. Receptor desensitization leads to an adaptive response to the ligand NE. For 4 PKA substrates (PLB, IKs, ICa, RyR) there are coherent positive feed-forward loops formed when PKA enhances substrate phosphorylation both directly (PKA phosphorylates the substrate) and indirectly by blocking inhibitor-1 (Inhib1), which phosphorylates protein phosphatase 1 (PP1), inhibiting substrate phosphorylation. We hypothesized that these feedback and feed-forward loops contributed to the predicted adaptive dynamics seen in our default simulations of PLB activity (see Figure 4A). To test this, we performed simulations in which each of these feedback or feed-forward loops were disrupted (i.e. setting their corresponding reactions weights, WGRK-B1ARPG, WPKAC-B1ARPA, or WPKAC-Inhib1, to zero) (Figure 5).

Figure 5
figure5

Dynamic roles for feedback and feed-forward loops in the β-adrenergic network. Time courses of PLB activation in response to transient NE exposures (NE = 1) with several network perturbations. The thick solid line is the response of the default model, showing moderate PLB activation and damped oscillations. Also shown are simulations in which either the inhibitor-1 coherent positive feed-forward loop (thick dashed line; WPKAC-Inhib1 = 0), the GRK negative feedback loop (thin dashed line; WGRK-B1ARPG = 0) or the PKAC negative feedback loop (thin solid line; WPKAC-B1ARPA = 0) were eliminated. Inhibitor-1 amplified PLB activation without affecting PLB dynamics. Both GRK and PKAC feedback loops decreased PLB phosphorylation via β1-adrenergic receptor desensitization, though PKAC controlled steady-state adaptation while GRK feedback dampened the PLB response.

Full size image

Surprisingly, disrupting individual loops had dramatically different consequences. Inhibition of the GRK negative feedback loop resulted in sustained oscillations, indicating that GRK negative feedback contributed to damping PLB phosphorylation. In contrast, inhibition of the PKA negative feedback loop raised the steady-state PLB activity and disrupted oscillations, showing that this feedback loop controlled steady-state PLB adaptation. Lastly, blocking the inhibitor-1 feed-forward loop reduced PLB activity without qualitatively affecting the timecourse, suggesting that this loop amplifies PLB signaling. Previous experimental and modeling studies comparing GRK and PKA feedback loops studied the β2-AR isoform, where receptor desensitization was driven primarily by GRK [26, 27]. However, steady-state measurements in cells expressing β1-AR [28], the receptor isoform considered here, are consistent with the current model predictions that both GRK and PKA feedbacks contribute significantly to PLB responses, via β1-adrenergic receptor desensitization. Thus, the normalized-Hill modeling approach can be used to assess how various network architectures drive signaling dynamics, though results can be refined when experimental data is available (discussed further in subsequent sections).

Quantitative sensitivity analysis of the β-adrenergic network

Quantitative sensitivity analysis provides an approach to systematically survey the functional relationships within a signaling network, often revealing unexpected systems properties. We performed sensitivity analysis on our normalized-Hill β-adrenergic model using the same approach as that implemented for our toy network (Figure 6A). For this analysis, model inputs were set at moderate activity levels to avoid saturation of the network's dynamic range (NE = 0.5, Fsk = 0.2, IBMX = 0.2). Species are generally ordered from upstream (top/left) to downstream (bottom/right) in the sensitivity matrix. From this analysis, several global properties of the network can be identified. The diagonal represents self-activation, which is not always prominent due to the negative feedback loops discussed earlier. While upstream species in the left portion of the matrix (NE through PKAC) affect many others, species further downstream in the pathway (PP2A through PP1) affect fewer. This indicates that the downstream species are in "branches" of the pathway and do not significantly feedback on upstream components.

Figure 6
figure6

Sensitivity analysis of the β-adrenergic network identifies key system regulators. A) Sensitivity analysis of the normalized-Hill β-adrenergic model. Note the bright green squares prevalent in the phosphodiesterase (PDE), protein kinase inhibitor (PKI), and phosphatase 2A (PP2A) columns, implicating the importance of these three species as key negative regulators of the system. B) Sensitivity analysis of a detailed biochemical model of the same pathway, allowing for direct comparison of the two modeling approaches. A Pearson correlation coefficient of 0.75 was computed between the individual sensitivities of the two models, demonstrating substantial quantitative agreement between the functional outputs of the two modeled networks.

Full size image

Importantly, this analysis reveals quantitative relationships that are not apparent from the network topology alone. Several inhibiting (green) relationships, mainly involving PDE, PKI, and PP2A, are quantitatively prominent. These proteins appear to be key negative regulators of this network that may serve as potential therapeutic targets. PDE strongly inhibits (green) PKAC, TnI and B1ARPA, while it modestly activates B1AR and AC. These effects are not explained completely by path length, as PDE has a smaller effect on its direct target (cAMP) than species which are three steps away (e.g. TnI, B1ARPA). Indeed, the PDE inhibitor milrinone has been used for patients with congestive heart failure [29]. However, milrinone actually worsens mortality by elevating occurrences of ventricular arrhythmias [30, 31], perhaps due to the large number of PDE-sensitive species suggested by our modeling results. PP2A is a strong inhibitor of TnI with less strong inhibition of IKs and PLB. In contrast, other perturbations such as Fsk and IBMX exhibit rather modest effects on the network despite affecting many other species in a qualitative sense.

Examining a particular row in this matrix allows one to identify perturbations that are more or less likely to affect a given output. For example, the calcium channel ICa appears very sensitive to activation by cAMP or inhibition by PDE, but is less sensitive to Fsk activation (same path length as PDE) or direct PP1 inhibition. Some species have low quantitative sensitivity even though they are within the same pathway. For example, direct substrates of PKAC (e.g. β1ARPA, ICa, PLB) are highly sensitive to perturbations in AC, cAMP and PKAC, while other species in the same pathway, such as β1AR or Gsα, are much less sensitive. Many of these quantitative predictions cannot be achieved by qualitative graph analysis, and may be used to prioritize future experiments.

Direct comparison of normalized-Hill and biochemical models

In order to assess the predictive accuracy of the normalized-Hill β-adrenergic model, the model's sensitivity matrix was compared to a similar matrix generated from a detailed biochemical model of the same signaling network [7] (Figure 6B, additional details in Methods). The difference between these two matrices is shown in Figure 6C. Overall, there are a number of similarities between the structures of the two sensitivity matrices. Note the clear divisions between upstream and downstream components, indicating similar predictions of global functional relationships across the two models. To quantify these similarities, we computed the Pearson correlation coefficient between the corresponding sensitivities of the two models and obtained a value of 0.75. This analysis indicates that there is substantial quantitative agreement between the two models, considering the significant differences in their formulation and the use of default parameters in the normalized-Hill model.

To test the appropriateness of our default selections for Hill coefficients and EC50 values, we examined correlation coefficients between the biochemical model and normalized-Hill models while varying default "n" and EC50. We found that the predictions were insensitive to the choice of default Hill coefficient (though n = 1.4 was optimal), but highly sensitive to the intuitive EC50 value of 0.5. Higher sensitivity to EC-50 arises because a linear pathway tends to systematically amplify or diminish signals when default EC-50 deviates from 0.5 (see Additional File 2, Figure S1).

Despite striking similarities between predictions by the normalized-Hill and biochemical models, there are also several notable discrepancies that may provide further insight. To highlight these, we re-classified the predicted sensitivities from both models as either "activating" (S = 1), "inhibiting" (S = -1), or "neutral" (S = 0) and produced qualitative sensitivity matrices using only these values (see Additional File 3, Figure S2, additional details in Methods). Globally, the two models showed good agreement in terms of individual sensitivity types, with 457 out of 484 (94%) individual sensitivities qualitatively matching. Of the 27 mismatches, 3 of these were in opposite directions (0.62% of the total). This was seen, for example, when Gsα/βγ (GsaBg) was perturbed: the normalized-Hill model predicted that reduction of Gsα/βγ increases activation of B1AR and B1ARPG, while the opposite result was obtained from biochemical model predictions. In the normalized-Hill model, reduced PKA feedback via reduced GsaBg enhanced fractional activation of B1AR and, thus, B1ARPG. These results are expected given that all normalized-Hill model interactions are unidirectional. On the other hand, the biochemical model uses detailed mass action kinetics to describe these interactions, where a reduction in total Gsα/βγ can actually pull additional free receptors to a bound form, resulting in reduced B1AR and B1ARPG. Thus, there are competing mechanisms in this portion of the network (PKA feedback versus G-protein activation) that are represented differently between the two model types. Related to this issue, most species in the biochemical network are sensitive to perturbations in Gβγ, though to a very small extent quantitatively (compare Gbg columns in Figures 6B and S2). This subtle difference arises because Gβγ is a terminal node in the normalized-Hill model, whereas the additional details of the biochemical model allow Gβγ perturbations to very modestly influence downstream signaling. Though quantitatively insignificant in most cases, these results highlight subtle limitations in the normalized-Hill β-adrenergic signaling model that can be addressed with additional reactions (though not done here).

Other discrepancies are attributed to differences in the extent to which spatial compartmentation was incorporated in the two models. For example, all downstream PKA substrates are sensitive to inhibitor-1 perturbation in the normalized-Hill model since this species inhibits global PP1 activity. In the biochemical model, however, inhibitor-1 is only responsible for local inhibition of PP1 near PLB (but not other PKA substrates; see Inhib1 columns in Figure 6). Thus, many of the discrepancies can be attributed to subtle differences in network connectivity rather than the modeling approaches themselves.

Other logic-based modeling approaches have also used logical AND/OR interactions and differential equations to describe biochemical networks [9, 15, 18]. While these logic-based approaches have not previously been directly compared to a biochemical model, we extended the sensitivity analysis (as in Figure 6) to examine these approaches as well. The β-adrenergic model was re-implemented using piece-wise linear, Hill, or linear activation functions. As shown in Additional File 4, Figure S3 and Additional File 5, Figure S4, the normalized-Hill approach exhibits substantially better performance compared with the piece-wise linear approach or Hill equations with previously-used default parameters (n = 3, K = 0.3) [18]. The poor performance of the piece-wise linear approach is largely due to the fact that its steady-state values are restricted, hindering predictions of sensitivity to a quantitative perturbation. The Hill approach performed poorly as well, but its performance could be improved somewhat by optimizing parameters "n" and "K". Linear activation functions worked fairly well for this steady-state sensitivity analysis, but linear activation functions are not able to predict nonlinear phenomena such as the bistability shown in Figure 1.

To further compare predictions from the normalized-Hill and biochemical models, we examined the predicted dynamics of Gsα and PLB in response to a transient NE exposure. Though global functional relationships are strikingly similar between the two models (Figure 6), the normalized-Hill model outputs using all default parameters contains damped oscillations that are not prominent in the biochemical model (compare Figure 4A with the left panel in Figure 7). To test whether the normalized-Hill model predictions can be refined further, we fit several parameters in the normalized-Hill model to time-course data from the biochemical model using a nonlinear least squares optimization algorithm (lsqnonlin in Matlab). Eleven parameters (4 τ's, 3 weights, 3 EC50's, and an additional basal receptor activity term) were adjusted to fit model predictions (see Additional file 1, Supplemental Methods).

Figure 7
figure7

Comparing predicted dynamics of the normalized-Hill model with a detailed biochemical model. The left panel shows normalized dynamic responses of Gsα and PLB phosphorylation predicted by the biochemical model during a transient 10 nM NE exposure. Normalized-Hill model predictions of the same outputs during a similar NE exposure are shown in the right panel. Activity levels and dynamics of the normalized-Hill model traces were fit to those of the biochemical model via least-squares parameter estimation. In all, 11 parameters were adjusted (4 τ's, 3 W's, 3 EC50's, and an additional basal receptor activity term) to mimic biochemical model outputs using the normalized-Hill model.

Full size image

These parameter adjustments allowed for more similar signaling dynamics and comparable peak fractional activities of Gsα and PLB compared with the biochemical model (Figure 7), with more gradual adaptation rather than damped oscillations. We further probed whether conclusions drawn regarding the feedforward and feedback loops from Figure 5 would be maintained in the adjusted model. As shown in Additional File 6, Figure S5, while the PLB response exhibits gradual adaptation rather than damped oscillations, the role of the feedforward and feedback loops are similar: Inhibitor 1 amplifies PLB, B1ARPG attenuates PLB, and B1ARPA dominates the degree of steady-state adaptation. These analyses demonstrate that the normalized-Hill model largely captures key features of the more detailed biochemical model, that specific model discrepancies can be quantitatively explained by parameter differences, and that dynamic predictions can be refined by fitting relevant parameters to available data.

Model extension to incorporate integrin-mediated mechanotransduction

Due to the wealth of data needed to parameterize a biochemical model, it can be difficult to extend the model with new experimental findings. However, the normalized-Hill modeling framework may be more easily extended. We searched the Pathway Commons database http://www.pathwaycommons.org/pc/ for additional PKA substrates, finding 46 pathways and 63 catalysis reactions. The first listed interaction was PKA-mediated phosphorylation of β-integrins, which regulates cellular adhesion [32]. Alenghat et al. recently showed that shear stresses applied to β1 integrins using RGD-coated magnetic microbeads stimulated cAMP via Gsα [33]. By integrating this new data into our model, we predicted that mechanical stresses may activate the cAMP/PKA pathway in cardiac myocytes in a manner independent of β-adrenergic receptors. Additional File 7, Figure S6 shows how these mechanisms were incorporated into the normalized-Hill model. We ran a simulation in which only mechanical stress was applied to the network (i.e. NE = Fsk = IBMX = 0), and saw activation of Gsα, cAMP, and integrin β-subunit phosphorylation (Itgbp) (Figure 8). Note that the response times of these outputs to the stress stimulus are dependent on their relative positions in the signaling cascade. As an independent validation of these predictions, prior work has shown that stretch of intact heart induces a rise in cAMP, which may contribute to stretch-induced increases in cardiac contractility [34].

Figure 8
figure8

Model extension to include integrin activation and phosphorylation in response to mechanical stress. Reactions for stress-induced integrin activation and subsequent phosphorylation by PKA were identified from recent protein interaction data http://www.pathwaycommons.org/pc/. Shown are predicted dynamics during a transient mechanical stress stimulus for Gsα (which is stimulated by the stress input), cAMP, and integrin phosphorylation (Itgbp) by PKA. The mechanisms for this additional pathway are schematized in Additional File 7, Figure S6.

Full size image

Discussion and conclusions

We developed a normalized-Hill differential equation modeling approach that combines advantages of both biochemical and Boolean models. Even when parameters are not available, this approach allows for predictions of signaling dynamics and is compatible with a wide range of existing systems analyses, including quantitative sensitivity analysis. Furthermore, these models can be iteratively refined by either tuning parameter values (as done here) or adding additional reactions to better reflect quantitative features of experimental data. This approach was evaluated by direct quantitative comparison with a well characterized and experimentally validated biochemical model of the cardiac β-adrenergic signaling network [7].

The analysis revealed several new insights into relationships between β-adrenergic network topology and dynamics. Dissection of the feedforward and feedback loops showed that each loop could play a unique role in regulating PLB phosphorylation dynamics. The Inhibitor-1-mediated coherent positive feedforward loop was predicted to amplify signaling, which is consistent with reports of decreased PLB phosphorylation in mice with Inhibitor-1 ablation [35]. In contrast, PKA- and GRK-mediated negative feedback loops (B1ARPA and B1ARPG, respectively) reduced PLB phosphorylation either by enhancing steady-state adaptation (PKA) or attenuating the overall PLB signal (GRK). While the kinetics and damped oscillations were sensitive to the rate constants for these feedback loops, the overall roles of the feedback loops were maintained using both default and adjusted parameters in the normalized-Hill model. Quantitative sensitivity analysis revealed the global structure of functional relationships in this network, which was not clear from previous analyses of the biochemical model [6]. This suggests that therapeutically targeting species within the PKA negative feedback loop may be less specific than targeting individual branches involving PKA substrates. Sensitivity analysis highlighted key hubs (such as PDE's) as major inhibitors of this pathway, consistent with drug development in this area [36]. Sensitivity analysis also revealed quantitative relationships that are not seen from the topology alone, such as lower sensitivity of Gsα and AC to perturbations in PKAC. Finally, model extension allowed analysis of crosstalk with integrin-mediated mechanotransduction, helping to explain experimentally-observed cAMP synthesis and increased cardiac contractility during stretch [34].

The normalized-Hill differential equation approach is related to logic-based approaches that also attempt to simulate biological networks with limited parameter data. The most well-characterized qualitative approach is Boolean network modeling, where species states are binary and the network is simulated over discrete time [37]. The binary state limitation can be partially overcome by allowing multiple discrete states [38], but this requires additional information to determine the number of states and to what state a given reaction activates a species. Discrete-time models can also require complex updating schemes to avoid artifacts [15], resulting in stochastic simulation data that increases the complexity of subsequent analysis. To address these issues, piecewise linear differential equations have been developed, which allow continuous-valued species states and continuous time [9, 15, 17]. Like our approach, piecewise linear differential equations use logical AND/OR operations and represent time-dependence with differential equations. However, this approach still uses discrete thresholds for activation of reactions that may not be consistent with experimental data, and the steady-state predictions of these models are binary [39]. As a result of limited possible steady-state values, we found that piecewise linear differential equations could not accurately predict quantitative model sensitivities. However, an advantage of the piecewise linear approach is that it requires fewer parameters: time constants and activation threshold parameters.

As opposed to discrete thresholds, a wide variety of biological networks exhibit smoothly saturating activation profiles that are well-approximated by Hill functions. Hill functions have been used to model a wide range of phenomena, including hemoglobin binding [40], dynamics of synthetic gene networks [41] and signaling pathways [42]. Because our state variables are generally limited to values between 0 and 1, we developed a normalization scheme that constrained activating Hill functions such that f(0) = 0, f(1) = 1 and f(EC50) = 0.5. We found that normalizing with just 2 of these 3 constraints, as done previously [18], caused undesirable shifts in EC50's or maximal activities, resulting in either diminished or enhanced signaling down a linear pathway. These artifacts limited the dynamic range of the pathway and hindered quantitative sensitivity analysis as compared with the biochemical model. In addition, we have included parameters that allow for further quantitative refinement of model predictions: reaction weights (W) and maximum activity values (YMAX). Reaction weights allow certain reactions to exert more influence than others. This was important for refinement of our β-adrenergic model because only partial receptor desensitization is seen experimentally [22]. Ability to perturb YMAX values was critical for performing quantitative sensitivity analysis in the present work, but this also allows for future incorporation of changes in protein expression that are independent of the level of fractional activation.

Direct quantitative comparison of the normalized-Hill β-adrenergic model with a validated biochemical model allowed characterization of the strengths and weaknesses of this approach in a controlled environment. To our knowledge, such a direct comparison has not been previously made for other logic-based approaches. Through sensitivity analysis and time course comparisons, we found that default parameter values were sufficient to predict global functional relationships with relatively high quantitative agreement, especially in terms of responses to systematic perturbations (correlation coefficient of 0.75). Indeed, the vast majority of experimental data on signaling networks provides relative rather than absolute quantification [43]. The normalized-Hill activation function outperformed piecewise linear, Hill and linear activation functions. Absolute quantitative predictions of particular time constants or signal magnitudes can be achieved with either prior parameter knowledge or optimization-based parameter estimation. Because only a small number of steps in a pathway are highly sensitive or rate-limiting, high quantitative agreement was obtained when fitting 11 parameters. By comparison, the biochemical model contains 88 parameters.

Because the model is composed of differential equations, this modeling framework is compatible with the wealth of analysis tools currently available from dynamical systems including sensitivity analysis (shown here), bifurcation analysis [44] and parameter estimation (performed in specific cases here) [45]. These models may also be directly integrated with multi-scale models of electrical-mechanical coupling or other aspects of physiology [46].

A number of differences from biochemical models must be considered regarding model structure. First, species states are quantified in terms of fractional activation rather than absolute quantities such as concentration. While this may be appropriate for proteins where the data is often normalized and the focus is on post-translational modifications, small molecules such as cyclic AMP are always in an "active" form. For such cases, fractional activation may be considered to be relative to the amount measured under a highly-stimulated condition. Second, not all relationships are well-represented by activating or inhibiting Hill functions. Other activation functions can be used as appropriate, and the use of differential equations allows a detailed biochemical module to be easily embedded within a larger normalized-Hill model or vice versa. Third, this normalized-Hill approach does not explicitly incorporate competitive inhibition, leading to a small number of incorrect predictions when Gβγ was perturbed. Finally, multiple post-translational modifications on a single protein may be required to be represented as separate species, as we have done for the GRK and PKA phosphorylation sites on the β-adrenergic receptor (B1ARPG and B1ARPA, respectively).

Despite these limitations, the normalized-Hill differential equation modeling approach predicts quantitative systems properties that typically require detailed biochemical models. Using default parameters, we obtained accurate simulations of the β-adrenergic signaling network that revealed unexpected functional relationships, generated experimentally testable predictions and readily incorporated new parameter and reaction data. This suggests that this approach may facilitate larger scale reconstruction and analysis of signaling networks, particularly for those where biochemical characterization is limited.

Methods

Normalized-Hill differential equation modeling approach

Species dynamics are predicted using ordinary differential equations, where the state variables represent fractional activation of each species. Species interactions were defined with normalized activating or inhibiting Hill functions (fact or finhib) which are described below. Pathway crosstalk was implemented using logical AND and OR operations: "f(x)f(y)" and "f(x)+f(y)-f(x)f(y)", respectively. Using this approach, the toy model described in Figure 1 was implemented with the following differential equations:

d D d t = 1 τ D ( W B D f a c t ( B ) D M A X − D ) d E d t = 1 τ E ( W C D E f a c t ( C ) f i n h i b ( D ) E M A X − E ) d C d t = 1 τ C [ ( W A C f a c t ( A ) + W E C f a c t ( E ) − W A C f a c t ( A ) W E C f a c t ( E ) ) C M A X − C ]
(1.1)

where τ is the time constant for a given species, W is the reaction weight (constrained to 0≤W≤1), and YMAX is the maximal fractional activation allowing simulations of knock-down (YMAX < 1) or overexpression (YMAX > 1). The normalized activating or inhibiting Hill functions, depicted graphically in Figure 1B, have the following form:

f a c t ( X ) = B X n K n + X n ; f i n h i b ( X ) = 1 − B X n K n + X n
(1.2)

where B and K are constrained such that fact(0) = 0, fact(EC50) = 0.5 and fact(1) = 1. From these constraints, we derived that:

B = E C 50 n − 1 2 E C 50 n − 1 ; K = ( B − 1 ) 1 / n
(1.3)

We further constrained fact(X) = 1 for X≥1 to ensure that species activities are limited to YMAX. As default parameters, we used W = 1, EC50 = 0.5, n = 1.4, τ = 1, and YMAX = 1. Sensitivity of β-adrenergic model predictions to the choice of default values for n and EC50 were examined (see Additional File 2, Figure S1). The model was insensitive to the default Hill coefficient, with Pearson correlation coefficients ≥ 0.93 for 1.01 ≤ n ≤1.7 compared to n = 1.4. Sensitivity to default EC50 was higher because a value other than 0.5 causes a linear pathway to systematically amplify or diminish signaling at subsequent steps. All differential equations were implemented in MATLAB (MathWorks, Natwick MA) and solved numerically with the ode23 function.

Sensitivity analysis

Two-dimensional sensitivity matrices were generated from normalized-Hill model simulations according to the equation Sij = (ΔYi/ΔPj)(Po,j/Yo,i), where S ij is the sensitivity of species "i" to perturbation of species "j", ΔYi is the change in steady-state output of species "i", and ΔPj is the change in parameter "j". These sensitivities are normalized to the original output (Yo,i) and parameter value (Po,j) to facilitate comparison between parameters and species. For a given numerical experiment (i.e. a column in a sensitivity matrix), a single species' YMAX was perturbed by ΔP (in general, -25%). The model was run to steady state and the normalized sensitivities in all model species (rows of the sensitivity matrix) were computed. Thus, the sensitivity matrix represents a total of ni*nj predictions from nj numerical experiments. We found that different perturbation magnitudes (-10 to -75%) produced very similar sensitivity matrices (correlation coefficient > 0.98).

Comparison of normalized-Hill model with biochemical model

Normalized-Hill model predictions were compared to our previously described biochemical model [7] in terms of parameter sensitivities and activation dynamics. Sensitivity analysis was computed as described for the normalized-Hill model, except that the biochemical model does not incorporate YMAX parameters. Therefore, the most analogous parameters were perturbed: total species concentrations (where explicitly defined) or catalytic rate constants associated with production of a particular output. Biochemical model outputs were also selected to be consistent with the normalized-Hill model species as closely as possible. For quantitative comparison of normalized-Hill and biochemical model sensitivity matrices, a Pearson correlation coefficient was calculated from the individual elements of the sensitivity matrices. For a more qualitative comparison of the two models, we classified individual elements of the sensitivity matrix as either "activating" (S = 1), "inhibiting" (S = -1) or "neutral" (S = 0), in which threshold sensitivities of ±0.003 (determined by visual inspection of the histograms of individual sensitivities from the two models) were used to reassign individual values. Comparison of normalized-Hill and biochemical models was performed by counting the number of qualitatively matching individual sensitivities and dividing by the total number of matrix elements.

Biochemical model activation time courses were performed in similar fashion to normalized-Hill simulations. Reported values are normalized to the maximal activation level for that particular output (e.g. concentration of phosphorylated PLB subunits over total PLB). Ligand concentration was set to 10 nM where 'NE' is labeled in figures depicting biochemical model outputs.

In some simulations (Figure 7), normalized-Hill model parameters were adjusted to match outputs from the biochemical model. A list of all adjusted parameters is provided in Additional file 1, Supplemental Methods, although the fitting procedure is described here. Time courses for Gsα and PLB phosphorylation were obtained from the biochemical model during a transient NE exposure and normalized to their maximum values (i.e. phosphorylated PLB over total PLB protein). Parameters in the normalized-Hill model were then optimized using nonlinear least-squares minimization (lsqnonlin in MATLAB). This function adjusts selected parameter values to minimize the sum of square errors between data points; in this case, errors between normalized-Hill and the normalized biochemical model predictions were minimized. Several rounds of fitting with randomized initial parameter guesses, along with some manual adjustments, produced the final set of parameters used to generate the normalized-Hill model curves in Figure 7. See Additional File 1, Supplemental Methods for further details.

Acknowledgements

The authors thank Dr. Jason Papin for critical review of the manuscript. This work was supported by National Institutes of Health grant HL094476 and the American Heart Association grant 0830470N.

References

  1. 1.

    Kiuchi K, Sato N, Shannon RP, Vatner DE, Morgan K, Vatner SF: Depressed beta-adrenergic receptor- and endothelium-mediated vasodilation in conscious dogs with heart failure. Circ Res. 1993, 73 (6): 1013-1023.

  2. 2.

    Brodde OE: Beta-adrenoceptors in cardiac disease. Pharmacol Ther. 1993, 60 (3): 405-430. 10.1016/0163-7258(93)90030-H

  3. 3.

    Bristow MR, Ginsburg R, Minobe W, Cubicciotti RS, Sageman WS, Lurie K, Billingham ME, Harrison DC, Stinson EB: Decreased catecholamine sensitivity and beta-adrenergic-receptor density in failing human hearts. N Engl J Med. 1982, 307 (4): 205-211. 10.1056/NEJM198207223070401

  4. 4.

    Brodde OE, Michel MC: Adrenergic and muscarinic receptors in the human heart. Pharmacol Rev. 1999, 51 (4): 651-690.

  5. 5.

    Engelhardt S, Bohm M, Erdmann E, Lohse MJ: Analysis of beta-adrenergic receptor mRNA levels in human ventricular biopsy specimens by quantitative polymerase chain reactions: progressive reduction of beta 1-adrenergic receptor mRNA in heart failure. J Am Coll Cardiol. 1996, 27 (1): 146-154. 10.1016/0735-1097(95)00425-4

  6. 6.

    Saucerman JJ, Brunton LL, Michailova AP, McCulloch AD: Modeling beta-adrenergic control of cardiac myocyte contractility in silico. J Biol Chem. 2003, 278 (48): 47997-48003. 10.1074/jbc.M308362200

  7. 7.

    Saucerman JJ, Healy SN, Belik ME, Puglisi JL, McCulloch AD: Proarrhythmic consequences of a KCNQ1 AKAP-binding domain mutation: computational models of whole cells and heterogeneous tissue. Circ Res. 2004, 95 (12): 1216-1224. 10.1161/01.RES.0000150055.06226.4e

  8. 8.

    Saucerman JJ, McCulloch AD: Cardiac beta-adrenergic signaling: from subcellular microdomains to heart failure. Ann N Y Acad Sci. 2006, 1080: 348-361. 10.1196/annals.1380.026

  9. 9.

    Albert I, Thakar J, Li S, Zhang R, Albert R: Boolean network simulations for life scientists. Source Code Biol Med. 2008, 3: 16- 10.1186/1751-0473-3-16

  10. 10.

    Aldridge BB, Saez-Rodriguez J, Muhlich JL, Sorger PK, Lauffenburger DA: Fuzzy logic analysis of kinase pathway crosstalk in TNF/EGF/insulin-induced signaling. PLoS Comput Biol. 2009, 5 (4): e1000340- 10.1371/journal.pcbi.1000340

  11. 11.

    Papin JA, Price ND, Palsson BO: Extreme pathway lengths and reaction participation in genome-scale metabolic networks. Genome Res. 2002, 12 (12): 1889-1900. 10.1101/gr.327702

  12. 12.

    Klamt S, Stelling J: Two approaches for metabolic pathway analysis?. Trends Biotechnol. 2003, 21 (2): 64-69. 10.1016/S0167-7799(02)00034-3

  13. 13.

    Klemm K, Bornholdt S: Stable and unstable attractors in Boolean networks. Phys Rev E Stat Nonlin Soft Matter Phys. 2005, 72 (5 Pt 2): 055101-

  14. 14.

    Ruths D, Muller M, Tseng JT, Nakhleh L, Ram PT: The signaling petri net-based simulator: a non-parametric strategy for characterizing the dynamics of cell-specific signaling networks. PLoS Comput Biol. 2008, 4 (2): e1000005- 10.1371/journal.pcbi.1000005

  15. 15.

    Chaves M, Sontag ED, Albert R: Methods of robustness analysis for Boolean models of gene control networks. Syst Biol (Stevenage). 2006, 153 (4): 154-167.

  16. 16.

    Morris MK, Saez-Rodriguez J, Sorger PK, Lauffenburger DA: Logic-based models for the analysis of cell signaling networks. Biochemistry. 2010, 49 (15): 3216-3224. 10.1021/bi902202q

  17. 17.

    Glass L, Kauffman SA: The logical analysis of continuous, non-linear biochemical control networks. J Theor Biol. 1973, 39 (1): 103-129. 10.1016/0022-5193(73)90208-7

  18. 18.

    Wittmann DM, Krumsiek J, Saez-Rodriguez J, Lauffenburger DA, Klamt S, Theis FJ: Transforming Boolean models to continuous models: methodology and application to T-cell receptor signaling. BMC Syst Biol. 2009, 3 (1): 98- 10.1186/1752-0509-3-98

  19. 19.

    Saucerman JJ, Zhang J, Martin JC, Peng LX, Stenbit AE, Tsien RY, McCulloch AD: Systems analysis of PKA-mediated phosphorylation gradients in live cardiac myocytes. Proc Natl Acad Sci USA. 2006, 103 (34): 12923-12928. 10.1073/pnas.0600137103

  20. 20.

    Bers DM: Cardiac excitation-contraction coupling. Nature. 2002, 415 (6868): 198-205. 10.1038/415198a

  21. 21.

    Vila Petroff MG, Egan JM, Wang X, Sollott SJ: Glucagon-like peptide-1 increases cAMP but fails to augment contraction in adult rat cardiac myocytes. Circ Res. 2001, 89 (5): 445-452. 10.1161/hh1701.095716

  22. 22.

    Hausdorff WP, Caron MG, Lefkowitz RJ: Turning off the signal: desensitization of beta-adrenergic receptor function. FASEB J. 1990, 4 (11): 2881-2889.

  23. 23.

    Hagemann D, Kuschel M, Kuramochi T, Zhu W, Cheng H, Xiao RP: Frequency-encoding Thr17 phospholamban phosphorylation is independent of Ser16 phosphorylation in cardiac myocytes. J Biol Chem. 2000, 275 (29): 22532-22536. 10.1074/jbc.C000253200

  24. 24.

    Laflamme MA, Becker PL: Do beta 2-adrenergic receptors modulate Ca2+ in adult rat ventricular myocytes?. Am J Physiol. 1998, 274 (4 Pt 2): H1308-1314.

  25. 25.

    Yakoby N, Lembong J, Schupbach T, Shvartsman SY: Drosophila eggshell is patterned by sequential action of feedforward and feedback loops. Development. 2008, 135 (2): 343-351. 10.1242/dev.008920

  26. 26.

    Xin W, Tran TM, Richter W, Clark RB, Rich TC: Roles of GRK and PDE4 activities in the regulation of beta2 adrenergic signaling. J Gen Physiol. 2008, 131 (4): 349-364. 10.1085/jgp.200709881

  27. 27.

    Violin JD, DiPilato LM, Yildirim N, Elston TC, Zhang J, Lefkowitz RJ: beta2-adrenergic receptor signaling and desensitization elucidated by quantitative modeling of real time cAMP dynamics. J Biol Chem. 2008, 283 (5): 2949-2961. 10.1074/jbc.M707009200

  28. 28.

    Rapacciuolo A, Suvarna S, Barki-Harrington L, Luttrell LM, Cong M, Lefkowitz RJ, Rockman HA: Protein kinase A and G protein-coupled receptor kinase phosphorylation mediates beta-1 adrenergic receptor endocytosis through different pathways. J Biol Chem. 2003, 278 (37): 35403-35411. 10.1074/jbc.M305675200

  29. 29.

    Packer M: Effect of phosphodiesterase inhibitors on survival of patients with chronic congestive heart failure. Am J Cardiol. 1989, 63 (2): 41A-45A. 10.1016/0002-9149(89)90392-5

  30. 30.

    Endoh M, Hori M: Acute heart failure: inotropic agents and their clinical uses. Expert Opin Pharmacother. 2006, 7 (16): 2179-2202. 10.1517/14656566.7.16.2179

  31. 31.

    Packer M, Carver JR, Rodeheffer RJ, Ivanhoe RJ, DiBianco R, Zeldis SM, Hendrix GH, Bommer WJ, Elkayam U, Kukin ML, et al.: Effect of oral milrinone on mortality in severe chronic heart failure. The PROMISE Study Research Group. N Engl J Med. 1991, 325 (21): 1468-1475. 10.1056/NEJM199111213252103

  32. 32.

    Jones SL: Protein kinase A regulates beta2 integrin avidity in neutrophils. J Leukoc Biol. 2002, 71 (6): 1042-1048.

  33. 33.

    Alenghat FJ, Tytell JD, Thodeti CK, Derrien A, Ingber DE: Mechanical control of cAMP signaling through integrins is mediated by the heterotrimeric Galphas protein. J Cell Biochem. 2009, 106 (4): 529-538. 10.1002/jcb.22001

  34. 34.

    Todaka K, Ogino K, Gu A, Burkhoff D: Effect of ventricular stretch on contractile strength, calcium transient, and cAMP in intact canine hearts. Am J Physiol. 1998, 274 (3 Pt 2): H990-1000.

  35. 35.

    Carr AN, Schmidt AG, Suzuki Y, del Monte F, Sato Y, Lanner C, Breeden K, Jing SL, Allen PB, Greengard P, et al.: Type 1 phosphatase, a negative regulator of cardiac function. Mol Cell Biol. 2002, 22 (12): 4124-4135. 10.1128/MCB.22.12.4124-4135.2002

  36. 36.

    Van Tassell BW, Radwanski P, Movsesian M, Munger MA: Combination therapy with beta-adrenergic receptor antagonists and phosphodiesterase inhibitors for chronic heart failure. Pharmacotherapy. 2008, 28 (12): 1523-1530. 10.1592/phco.28.12.1523

  37. 37.

    Pal R, Datta A, Bittner ML, Dougherty ER: Intervention in context-sensitive probabilistic Boolean networks. Bioinformatics. 2005, 21 (7): 1211-1218. 10.1093/bioinformatics/bti131

  38. 38.

    Pomerance A, Ott E, Girvan M, Losert W: The effect of network topology on the stability of discrete state models of genetic control. Proc Natl Acad Sci USA. 2009, 106 (20): 8209-8214. 10.1073/pnas.0900142106

  39. 39.

    Albert R, Wang RS: Discrete dynamic modeling of cellular signaling networks. Methods in enzymology. 2009, 467: 281-306. full_text

  40. 40.

    Saroff HA, Yap WT: The action of hemoglobin. Cooperative effects in tetrameric proteins. Biopolymers. 1972, 11 (5): 957-971. 10.1002/bip.1972.360110504

  41. 41.

    Rosenfeld N, Young JW, Alon U, Swain PS, Elowitz MB: Accurate prediction of gene feedback circuit behavior from component properties. Mol Syst Biol. 2007, 3: 143- 10.1038/msb4100185

  42. 42.

    Aldridge BB, Burke JM, Lauffenburger DA, Sorger PK: Physicochemical modelling of cell signalling pathways. Nat Cell Biol. 2006, 8 (11): 1195-1203. 10.1038/ncb1497

  43. 43.

    Malmstrom J, Beck M, Schmidt A, Lange V, Deutsch EW, Aebersold R: Proteome-wide cellular protein concentrations of the human pathogen Leptospira interrogans. Nature. 2009, 460 (7256): 762-765. 10.1038/nature08184

  44. 44.

    Higham CF: Bifurcation analysis informs Bayesian inference in the Hes1 feedback loop. BMC Syst Biol. 2009, 3: 12- 10.1186/1752-0509-3-12

  45. 45.

    Bandara S, Schloder JP, Eils R, Bock HG, Meyer T: Optimal experimental design for parameter estimation of a cell signaling model. PLoS Comput Biol. 2009, 5 (11): e1000558- 10.1371/journal.pcbi.1000558

  46. 46.

    Bassingthwaighte J, Hunter P, Noble D: The Cardiac Physiome: perspectives for the future. Exp Physiol. 2009, 94 (5): 597-605. 10.1113/expphysiol.2008.044099

Download references

Author information

Affiliations

  1. Department of Biomedical Engineering, University of Virginia, Charlottesville, VA, 22908, USA
    • Matthew J Kraeutler
    • , Anthony R Soltis
    •  & Jeffrey J Saucerman
Authors
  1. Matthew J Kraeutler
    View author publications
    You can also search for this author in
  2. Anthony R Soltis
    View author publications
    You can also search for this author in
  3. Jeffrey J Saucerman
    View author publications
    You can also search for this author in

Corresponding author

Correspondence to Jeffrey J Saucerman.

Additional information

Authors' contributions

JS planned the project; MK and AS performed the analysis; MK, AS and JS wrote the manuscript. All authors read and approved the final version of the manuscript.

Electronic supplementary material

Additional file 1: Supplemental Methods. Model equations and description of parameter estimation methods. (PDF 185 KB)

Additional file 2: Figure S1: Dependence of model accuracy on choices for default parameters. A) Varying the choice of default Hill coefficient has little effect on the predictive capability of the normalized-Hill model compared to the detailed biochemical model, though n = 1.4 is optimal. B) Strength of normalized-Hill predictions are very sensitive to the choice of default EC50, where model predictions dramatically worsened as default EC50 deviated from the intuitive value of 0.5. (PDF 24 KB)

Additional file 3: Figure S2: Qualitative comparison of sensitivity matrices from normalized-Hill and biochemical β-adrenergic models. A) Individual elements of each sensitivity matrix were re-classified as "activating" (red), "inhibiting" (green) or "neutral" (black), in which a threshold sensitivity of +/- 0.003 was determined by visual inspection of a histogram of sensitivities for each model. B) Difference plot showing regions of qualitative discrepancy between normalized-Hill model and the biochemical model. Globally, the two models showed good agreement in terms of individual sensitivity types, with 457 out of 484 (94%) individual sensitivities qualitatively matching. Of the 27 mismatches, only 3 sensitivities were in opposite directions (difference of -2; 0.62% of the total, see main text for further details). (PDF 38 KB)

Additional file 4: Figure S3: Comparison of biochemical model with normalized-Hill model and alternative modeling implementations. Differential equation models of the same β-adrenergic network were generated using linear, piece-wise linear, or traditional Hill activation functions. All other parts of the model were kept constant, including time constants, reaction weights, AND/OR logic, and network topology. A) Pearson correlation coefficients were computed by comparing the sensitivity matrix of each model against the sensitivity matrix of the biochemical model. The normalized-Hill approach produces the highest level of agreement. B, left) Sensitivity matrix using Hill activation functions of n = 3 and EC50 = 0.3, as used previously for default parameters a for T-cell signaling network [17]. B, right) Difference between the sensitivity matrices from the traditional Hill model (n = 3, EC50 = 0.3) and the biochemical model, indicating areas of discrepancy. Sensitivity matrices for remaining modeling approaches are shown in Additional File 5, Figure S4. (PDF 41 KB)

Additional file 5: Figure S4: Sensitivity matrices and differences from the biochemical model for β-adrenergic models implemented with "piece-wise linear", Hill, and linear activation functions. A) Traditional Hill activation with n = 2 and EC50 = 0.5. B) Activation was implemented as fact(x < 0.5) = 0, fact(x≥0.5) = 1 as described for "piece-wise linear differential equations" in [9, 15]. White indicates values with no Real solution (NaN), which were excluded from computation of the correlation coefficient. The poor agreement is related to the limited number of possible steady-state values using this activation function. C) Linear activation functions, implemented as fact(x) = x. (PDF 69 KB)

Additional file 6: Figure S5: Dynamic roles for feedback and feed-forward loops in the β-adrenergic network in the adjusted normalized-Hill Model. Using the adjusted model described in Figure 7, perturbations to feedback and feed-forward loops were performed as in Figure 5. The thick solid line is the response of the adjusted model to transient NE exposure, showing moderate PLB activation and partial adaptation to constant input. Also shown are simulations in which either the inhibitor-1 coherent positive feed-forward loop (thick dashed line; WPKAC-Inhib1 = 0), the GRK negative feedback loop (thin dashed line; WGRK-B1ARPG = 0) or the PKAC negative feedback loop (thin solid line; WPKAC-B1ARPA = 0) were eliminated. Inhibitor-1 amplified PLB activation without affecting PLB dynamics. Both GRK and PKAC feedback loops decreased PLB phosphorylation via β1-adrenergic receptor desensitization, though PKAC controlled steady-state adaptation while GRK feedback attenuated the PLB response. While the adjusted model does not exhibit the damped PLB oscillations of the default model (Figure 5), the roles of these feed-forward and feedback loops were largely maintained. (PDF 27 KB)

12918_2010_568_MOESM7_ESM.PDF

Additional file 7: Figure S6: Schematic of the β-adrenergic network with addition of integrin signaling. Based on recent protein interaction data http://www.pathwaycommons.org/pc/, reactions involving mechanical stress, RGD-coupled microbeads, integrins (Itg), and PKA-mediated integrin phosphorylation (Itgbp) were added to the model to demonstrate the extensibility of the normalized-Hill modeling approach. (PDF 102 KB)

Authors’ original submitted files for images

Rights and permissions

Reprints and Permissions

About this article

Cite this article

Kraeutler, M.J., Soltis, A.R. & Saucerman, J.J. Modeling cardiac β-adrenergic signaling with normalized-Hill differential equations: comparison with a biochemical model. BMC Syst Biol 4, 157 (2010). https://doi.org/10.1186/1752-0509-4-157

Download citation

Keywords

  • Hill Function
  • Biochemical Model
  • Linear Activation Function
  • Reaction Weight
  • Piecewise Linear Approach

BMC Systems Biology

ISSN: 1752-0509

Contact us