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.
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.
β-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].
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].
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).
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.
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).
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].