Abstract
Background
New approaches are needed for largescale predictive modeling of cellular signaling networks. While mass action and enzyme kinetic approaches require extensive biochemical data, current logicbased approaches are used primarily for qualitative predictions and have lacked direct quantitative comparison with biochemical models.
Results
We developed a logicbased 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 logicbased 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 integrinmediated mechanotransduction.
Conclusions
The normalizedHill 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 [15]. 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 [68]. 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 largescale 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 steadystate solution space of a network, its ability to predict dynamic timecourses for given experiments is limited [12]. Simulations from discretelevel models (e.g. Boolean) can be difficult to interpret due to sensitivity of model predictions to temporal updating schemes [13], assignment of discrete activitylevels to continuousvalued variables like concentration [14], and the limited ability to describe realistic timescales [15]. The tradeoffs inherent in many of these logicbased 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. Piecewiselinear differential equation models overcome some of these limitations by making both species values and time continuous, but steadystate 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 steadystate differences between Boolean and continuous models.
To address these limitations, we developed a normalizedHill 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 feedforward loops in the βadrenergic pathway and examined potential crosstalk with integrinmediated mechanotransduction. The analysis presented here demonstrates that the normalizedHill 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 normalizedHill 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. NormalizedHill toy network model. A) Schematic of the 5species toy network, including two inputs, an AND reaction, and a positive feedback loop. B) Characteristics of sample normalizedHill 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 EC_{50 }= 0.5.
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", halfmaximal effective concentration "EC_{50}", 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. EC_{50 }is the fractional activation of an input species required to induce halfmaximal 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(EC_{50}) = 0.5. Additionally, species activities are controlled by the parameters τ and Y_{MAX}, which are the reaction time constant and species maximal fractional activation, respectively. While Y_{MAX }= 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, EC_{50 }= 0.5, n = 1.4, τ = 1, and Y_{MAX }= 1; choice of default values for EC_{50 }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).
Additional file 1. Supplemental Methods. Model equations and description of parameter estimation methods.
Format: PDF Size: 185KB Download file
This file can be viewed with: Adobe Acrobat Reader
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, W_{D}, K_{B }and n_{B }could all be determined directly from a single steadystate 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 livecell 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, wellestablished 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. "A_{MAX}") and the subsequent changes in all network species' activities are quantified. The normalized steadystate sensitivities of these species to parameter perturbations are quantified according to S = (ΔY/ΔP)(P_{o}/Y_{o}) (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 logicbased or topological approaches.
Figure 2. 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 steadystate sensitivities of all model species were computed according to S = (ΔY/ΔP)(P_{o}/Y_{o}). The model predicts that E is positively regulated by A and C but negatively regulated by B and D.
βadrenergic signaling network
We next applied our normalizedHill 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 inhibitor1 (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 normalizedHill model was designed to be very similar to an extensively characterized and validated biochemical model of this same network [7].
Figure 3. Schematic of normalizedHill βadrenergic signaling network. The normalizedHill model consists of 25 species and 36 reactions. While developed using a normalizedHill differential equation approach, the network topology closely mimics a wellvalidated biochemical model [6] to allow for direct comparisons.
To test the normalizedhill β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 β_{1}ARs (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. Predicted dynamics and dose responses of the normalizedHill β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) Steadystate concentration responses of Gsα, cAMP and PLB to varying levels of NE.
Role of feedback and forward loops in shaping network dynamics
Feedback and feedforward 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 Gprotein 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 feedforward loops formed when PKA enhances substrate phosphorylation both directly (PKA phosphorylates the substrate) and indirectly by blocking inhibitor1 (Inhib1), which phosphorylates protein phosphatase 1 (PP1), inhibiting substrate phosphorylation. We hypothesized that these feedback and feedforward 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 feedforward loops were disrupted (i.e. setting their corresponding reactions weights, W_{GRKB1ARPG}, W_{PKACB1ARPA}, or W_{PKACInhib1}, to zero) (Figure 5).
Figure 5. Dynamic roles for feedback and feedforward 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 inhibitor1 coherent positive feedforward loop (thick dashed line; W_{PKACInhib1 }= 0), the GRK negative feedback loop (thin dashed line; W_{GRKB1ARPG }= 0) or the PKAC negative feedback loop (thin solid line; W_{PKACB1ARPA }= 0) were eliminated. Inhibitor1 amplified PLB activation without affecting PLB dynamics. Both GRK and PKAC feedback loops decreased PLB phosphorylation via β_{1}adrenergic receptor desensitization, though PKAC controlled steadystate adaptation while GRK feedback dampened the PLB response.
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 steadystate PLB activity and disrupted oscillations, showing that this feedback loop controlled steadystate PLB adaptation. Lastly, blocking the inhibitor1 feedforward 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, steadystate 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 normalizedHill 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 normalizedHill β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 selfactivation, 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. Sensitivity analysis of the βadrenergic network identifies key system regulators. A) Sensitivity analysis of the normalizedHill β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.
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 PDEsensitive 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. β_{1}ARPA, ICa, PLB) are highly sensitive to perturbations in AC, cAMP and PKAC, while other species in the same pathway, such as β_{1}AR 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 normalizedHill and biochemical models
In order to assess the predictive accuracy of the normalizedHill β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 normalizedHill model.
To test the appropriateness of our default selections for Hill coefficients and EC_{50 }values, we examined correlation coefficients between the biochemical model and normalizedHill models while varying default "n" and EC_{50}. 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 EC_{50 }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).
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 normalizedHill model compared to the detailed biochemical model, though n = 1.4 is optimal. B) Strength of normalizedHill 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.
Format: PDF Size: 24KB Download file
This file can be viewed with: Adobe Acrobat Reader
Despite striking similarities between predictions by the normalizedHill and biochemical models, there are also several notable discrepancies that may provide further insight. To highlight these, we reclassified 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 normalizedHill model predicted that reduction of Gsα/βγ increases activation of B1AR and B1ARPG, while the opposite result was obtained from biochemical model predictions. In the normalizedHill model, reduced PKA feedback via reduced GsaBg enhanced fractional activation of B1AR and, thus, B1ARPG. These results are expected given that all normalizedHill 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 Gprotein 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 normalizedHill 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 normalizedHill βadrenergic signaling model that can be addressed with additional reactions (though not done here).
Additional file 3. Figure S2: Qualitative comparison of sensitivity matrices from normalizedHill and biochemical βadrenergic models. A) Individual elements of each sensitivity matrix were reclassified 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 normalizedHill 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).
Format: PDF Size: 38KB Download file
This file can be viewed with: Adobe Acrobat Reader
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 inhibitor1 perturbation in the normalizedHill model since this species inhibits global PP1 activity. In the biochemical model, however, inhibitor1 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 logicbased modeling approaches have also used logical AND/OR interactions and differential equations to describe biochemical networks [9,15,18]. While these logicbased 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 reimplemented using piecewise linear, Hill, or linear activation functions. As shown in Additional File 4, Figure S3 and Additional File 5, Figure S4, the normalizedHill approach exhibits substantially better performance compared with the piecewise linear approach or Hill equations with previouslyused default parameters (n = 3, K = 0.3) [18]. The poor performance of the piecewise linear approach is largely due to the fact that its steadystate 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 steadystate sensitivity analysis, but linear activation functions are not able to predict nonlinear phenomena such as the bistability shown in Figure 1.
Additional file 4. Figure S3: Comparison of biochemical model with normalizedHill model and alternative modeling implementations. Differential equation models of the same βadrenergic network were generated using linear, piecewise 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 normalizedHill 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 Tcell 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.
Format: PDF Size: 41KB Download file
This file can be viewed with: Adobe Acrobat Reader
Additional file 5. Figure S4: Sensitivity matrices and differences from the biochemical model for βadrenergic models implemented with "piecewise linear", Hill, and linear activation functions. A) Traditional Hill activation with n = 2 and EC_{50 }= 0.5. B) Activation was implemented as f_{act}(x < 0.5) = 0, f_{act}(x≥0.5) = 1 as described for "piecewise 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 steadystate values using this activation function. C) Linear activation functions, implemented as f_{act}(x) = x.
Format: PDF Size: 69KB Download file
This file can be viewed with: Adobe Acrobat Reader
To further compare predictions from the normalizedHill 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 normalizedHill 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 normalizedHill model predictions can be refined further, we fit several parameters in the normalizedHill model to timecourse data from the biochemical model using a nonlinear least squares optimization algorithm (lsqnonlin in Matlab). Eleven parameters (4 τ's, 3 weights, 3 EC_{50}'s, and an additional basal receptor activity term) were adjusted to fit model predictions (see Additional file 1, Supplemental Methods).
Figure 7. Comparing predicted dynamics of the normalizedHill 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. NormalizedHill model predictions of the same outputs during a similar NE exposure are shown in the right panel. Activity levels and dynamics of the normalizedHill model traces were fit to those of the biochemical model via leastsquares parameter estimation. In all, 11 parameters were adjusted (4 τ's, 3 W's, 3 EC_{50}'s, and an additional basal receptor activity term) to mimic biochemical model outputs using the normalizedHill model.
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 steadystate adaptation. These analyses demonstrate that the normalizedHill 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.
Additional file 6. Figure S5: Dynamic roles for feedback and feedforward loops in the βadrenergic network in the adjusted normalizedHill Model. Using the adjusted model described in Figure 7, perturbations to feedback and feedforward 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 inhibitor1 coherent positive feedforward loop (thick dashed line; WPKACInhib1 = 0), the GRK negative feedback loop (thin dashed line; WGRKB1ARPG = 0) or the PKAC negative feedback loop (thin solid line; WPKACB1ARPA = 0) were eliminated. Inhibitor1 amplified PLB activation without affecting PLB dynamics. Both GRK and PKAC feedback loops decreased PLB phosphorylation via β1adrenergic receptor desensitization, though PKAC controlled steadystate 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 feedforward and feedback loops were largely maintained.
Format: PDF Size: 27KB Download file
This file can be viewed with: Adobe Acrobat Reader
Model extension to incorporate integrinmediated 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 normalizedHill modeling framework may be more easily extended. We searched the Pathway Commons database http://www.pathwaycommons.org/pc/ webcite for additional PKA substrates, finding 46 pathways and 63 catalysis reactions. The first listed interaction was PKAmediated phosphorylation of βintegrins, which regulates cellular adhesion [32]. Alenghat et al. recently showed that shear stresses applied to β_{1 }integrins using RGDcoated 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 normalizedHill 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 stretchinduced increases in cardiac contractility [34].
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/ webcite, reactions involving mechanical stress, RGDcoupled microbeads, integrins (Itg), and PKAmediated integrin phosphorylation (Itgbp) were added to the model to demonstrate the extensibility of the normalizedHill modeling approach.
Format: PDF Size: 102KB Download file
This file can be viewed with: Adobe Acrobat Reader
Figure 8. Model extension to include integrin activation and phosphorylation in response to mechanical stress. Reactions for stressinduced integrin activation and subsequent phosphorylation by PKA were identified from recent protein interaction data http://www.pathwaycommons.org/pc/ webcite. 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.
Discussion and conclusions
We developed a normalizedHill 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 Inhibitor1mediated coherent positive feedforward loop was predicted to amplify signaling, which is consistent with reports of decreased PLB phosphorylation in mice with Inhibitor1 ablation [35]. In contrast, PKA and GRKmediated negative feedback loops (B1ARPA and B1ARPG, respectively) reduced PLB phosphorylation either by enhancing steadystate 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 normalizedHill 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 integrinmediated mechanotransduction, helping to explain experimentallyobserved cAMP synthesis and increased cardiac contractility during stretch [34].
The normalizedHill differential equation approach is related to logicbased approaches that also attempt to simulate biological networks with limited parameter data. The most wellcharacterized 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. Discretetime 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 continuousvalued species states and continuous time [9,15,17]. Like our approach, piecewise linear differential equations use logical AND/OR operations and represent timedependence with differential equations. However, this approach still uses discrete thresholds for activation of reactions that may not be consistent with experimental data, and the steadystate predictions of these models are binary [39]. As a result of limited possible steadystate 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 wellapproximated 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(EC_{50}) = 0.5. We found that normalizing with just 2 of these 3 constraints, as done previously [18], caused undesirable shifts in EC_{50}'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 (Y_{MAX}). 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 Y_{MAX }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 normalizedHill β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 logicbased 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 normalizedHill 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 optimizationbased parameter estimation. Because only a small number of steps in a pathway are highly sensitive or ratelimiting, 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 multiscale models of electricalmechanical 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 posttranslational 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 highlystimulated condition. Second, not all relationships are wellrepresented 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 normalizedHill model or vice versa. Third, this normalizedHill approach does not explicitly incorporate competitive inhibition, leading to a small number of incorrect predictions when G_{βγ }was perturbed. Finally, multiple posttranslational 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 normalizedHill 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
NormalizedHill 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 (f_{act }or f_{inhib}) 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:
where τ is the time constant for a given species, W is the reaction weight (constrained to 0≤W≤1), and Y_{MAX }is the maximal fractional activation allowing simulations of knockdown (Y_{MAX }< 1) or overexpression (Y_{MAX }> 1). The normalized activating or inhibiting Hill functions, depicted graphically in Figure 1B, have the following form:
where B and K are constrained such that f_{act}(0) = 0, f_{act}(EC_{50}) = 0.5 and f_{act}(1) = 1. From these constraints, we derived that:
We further constrained f_{act}(X) = 1 for X≥1 to ensure that species activities are limited to Y_{MAX}. As default parameters, we used W = 1, EC_{50 }= 0.5, n = 1.4, τ = 1, and Y_{MAX }= 1. Sensitivity of βadrenergic model predictions to the choice of default values for n and EC_{50 }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 EC_{50 }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
Twodimensional sensitivity matrices were generated from normalizedHill model simulations according to the equation S_{ij }= (ΔY_{i}/ΔP_{j})(P_{o,j}/Y_{o,i}), where S_{ij }is the sensitivity of species "i" to perturbation of species "j", ΔY_{i }is the change in steadystate output of species "i", and ΔP_{j }is the change in parameter "j". These sensitivities are normalized to the original output (Y_{o,i}) and parameter value (P_{o,j}) to facilitate comparison between parameters and species. For a given numerical experiment (i.e. a column in a sensitivity matrix), a single species' Y_{MAX }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 n_{i}*n_{j }predictions from n_{j }numerical experiments. We found that different perturbation magnitudes (10 to 75%) produced very similar sensitivity matrices (correlation coefficient > 0.98).
Comparison of normalizedHill model with biochemical model
NormalizedHill 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 normalizedHill model, except that the biochemical model does not incorporate Y_{MAX }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 normalizedHill model species as closely as possible. For quantitative comparison of normalizedHill 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 normalizedHill 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 normalizedHill 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), normalizedHill 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 normalizedHill model were then optimized using nonlinear leastsquares 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 normalizedHill 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 normalizedHill model curves in Figure 7. See Additional File 1, Supplemental Methods for further details.
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.
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

Kiuchi K, Sato N, Shannon RP, Vatner DE, Morgan K, Vatner SF: Depressed betaadrenergic receptor and endotheliummediated vasodilation in conscious dogs with heart failure.
Circ Res 1993, 73(6):10131023. PubMed Abstract  Publisher Full Text

Brodde OE: Betaadrenoceptors in cardiac disease.
Pharmacol Ther 1993, 60(3):405430. PubMed Abstract  Publisher Full Text

Bristow MR, Ginsburg R, Minobe W, Cubicciotti RS, Sageman WS, Lurie K, Billingham ME, Harrison DC, Stinson EB: Decreased catecholamine sensitivity and betaadrenergicreceptor density in failing human hearts.
N Engl J Med 1982, 307(4):205211. PubMed Abstract  Publisher Full Text

Brodde OE, Michel MC: Adrenergic and muscarinic receptors in the human heart.
Pharmacol Rev 1999, 51(4):651690. PubMed Abstract  Publisher Full Text

Engelhardt S, Bohm M, Erdmann E, Lohse MJ: Analysis of betaadrenergic receptor mRNA levels in human ventricular biopsy specimens by quantitative polymerase chain reactions: progressive reduction of beta 1adrenergic receptor mRNA in heart failure.
J Am Coll Cardiol 1996, 27(1):146154. PubMed Abstract  Publisher Full Text

Saucerman JJ, Brunton LL, Michailova AP, McCulloch AD: Modeling betaadrenergic control of cardiac myocyte contractility in silico.
J Biol Chem 2003, 278(48):4799748003. PubMed Abstract  Publisher Full Text

Saucerman JJ, Healy SN, Belik ME, Puglisi JL, McCulloch AD: Proarrhythmic consequences of a KCNQ1 AKAPbinding domain mutation: computational models of whole cells and heterogeneous tissue.
Circ Res 2004, 95(12):12161224. PubMed Abstract  Publisher Full Text

Saucerman JJ, McCulloch AD: Cardiac betaadrenergic signaling: from subcellular microdomains to heart failure.
Ann N Y Acad Sci 2006, 1080:348361. PubMed Abstract  Publisher Full Text

Albert I, Thakar J, Li S, Zhang R, Albert R: Boolean network simulations for life scientists.
Source Code Biol Med 2008, 3:16. PubMed Abstract  BioMed Central Full Text  PubMed Central Full Text

Aldridge BB, SaezRodriguez J, Muhlich JL, Sorger PK, Lauffenburger DA: Fuzzy logic analysis of kinase pathway crosstalk in TNF/EGF/insulininduced signaling.
PLoS Comput Biol 2009, 5(4):e1000340. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Papin JA, Price ND, Palsson BO: Extreme pathway lengths and reaction participation in genomescale metabolic networks.
Genome Res 2002, 12(12):18891900. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Klamt S, Stelling J: Two approaches for metabolic pathway analysis?
Trends Biotechnol 2003, 21(2):6469. PubMed Abstract  Publisher Full Text

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. PubMed Abstract  Publisher Full Text

Ruths D, Muller M, Tseng JT, Nakhleh L, Ram PT: The signaling petri netbased simulator: a nonparametric strategy for characterizing the dynamics of cellspecific signaling networks.
PLoS Comput Biol 2008, 4(2):e1000005. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

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

Morris MK, SaezRodriguez J, Sorger PK, Lauffenburger DA: Logicbased models for the analysis of cell signaling networks.
Biochemistry 2010, 49(15):32163224. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Glass L, Kauffman SA: The logical analysis of continuous, nonlinear biochemical control networks.
J Theor Biol 1973, 39(1):103129. PubMed Abstract  Publisher Full Text

Wittmann DM, Krumsiek J, SaezRodriguez J, Lauffenburger DA, Klamt S, Theis FJ: Transforming Boolean models to continuous models: methodology and application to Tcell receptor signaling.
BMC Syst Biol 2009, 3(1):98. PubMed Abstract  BioMed Central Full Text  PubMed Central Full Text

Saucerman JJ, Zhang J, Martin JC, Peng LX, Stenbit AE, Tsien RY, McCulloch AD: Systems analysis of PKAmediated phosphorylation gradients in live cardiac myocytes.
Proc Natl Acad Sci USA 2006, 103(34):1292312928. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Bers DM: Cardiac excitationcontraction coupling.
Nature 2002, 415(6868):198205. PubMed Abstract  Publisher Full Text

Vila Petroff MG, Egan JM, Wang X, Sollott SJ: Glucagonlike peptide1 increases cAMP but fails to augment contraction in adult rat cardiac myocytes.
Circ Res 2001, 89(5):445452. PubMed Abstract  Publisher Full Text

Hausdorff WP, Caron MG, Lefkowitz RJ: Turning off the signal: desensitization of betaadrenergic receptor function.
FASEB J 1990, 4(11):28812889. PubMed Abstract

Hagemann D, Kuschel M, Kuramochi T, Zhu W, Cheng H, Xiao RP: Frequencyencoding Thr17 phospholamban phosphorylation is independent of Ser16 phosphorylation in cardiac myocytes.
J Biol Chem 2000, 275(29):2253222536. PubMed Abstract  Publisher Full Text

Laflamme MA, Becker PL: Do beta 2adrenergic receptors modulate Ca2+ in adult rat ventricular myocytes?
Am J Physiol 1998, 274(4 Pt 2):H13081314. PubMed Abstract  Publisher Full Text

Yakoby N, Lembong J, Schupbach T, Shvartsman SY: Drosophila eggshell is patterned by sequential action of feedforward and feedback loops.
Development 2008, 135(2):343351. PubMed Abstract  Publisher Full Text

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):349364. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Violin JD, DiPilato LM, Yildirim N, Elston TC, Zhang J, Lefkowitz RJ: beta2adrenergic receptor signaling and desensitization elucidated by quantitative modeling of real time cAMP dynamics.
J Biol Chem 2008, 283(5):29492961. PubMed Abstract  Publisher Full Text

Rapacciuolo A, Suvarna S, BarkiHarrington L, Luttrell LM, Cong M, Lefkowitz RJ, Rockman HA: Protein kinase A and G proteincoupled receptor kinase phosphorylation mediates beta1 adrenergic receptor endocytosis through different pathways.
J Biol Chem 2003, 278(37):3540335411. PubMed Abstract  Publisher Full Text

Packer M: Effect of phosphodiesterase inhibitors on survival of patients with chronic congestive heart failure.
Am J Cardiol 1989, 63(2):41A45A. PubMed Abstract  Publisher Full Text

Endoh M, Hori M: Acute heart failure: inotropic agents and their clinical uses.
Expert Opin Pharmacother 2006, 7(16):21792202. PubMed Abstract  Publisher Full Text

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):14681475. PubMed Abstract  Publisher Full Text

Jones SL: Protein kinase A regulates beta2 integrin avidity in neutrophils.
J Leukoc Biol 2002, 71(6):10421048. PubMed Abstract  Publisher Full Text

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):529538. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

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):H9901000. PubMed Abstract  Publisher Full Text

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):41244135. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Van Tassell BW, Radwanski P, Movsesian M, Munger MA: Combination therapy with betaadrenergic receptor antagonists and phosphodiesterase inhibitors for chronic heart failure.
Pharmacotherapy 2008, 28(12):15231530. PubMed Abstract  Publisher Full Text

Pal R, Datta A, Bittner ML, Dougherty ER: Intervention in contextsensitive probabilistic Boolean networks.
Bioinformatics 2005, 21(7):12111218. PubMed Abstract  Publisher Full Text

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):82098214. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Albert R, Wang RS: Discrete dynamic modeling of cellular signaling networks.
Methods in enzymology 2009, 467:281306. PubMed Abstract  Publisher Full Text

Saroff HA, Yap WT: The action of hemoglobin. Cooperative effects in tetrameric proteins.
Biopolymers 1972, 11(5):957971. PubMed Abstract  Publisher Full Text

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. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Aldridge BB, Burke JM, Lauffenburger DA, Sorger PK: Physicochemical modelling of cell signalling pathways.
Nat Cell Biol 2006, 8(11):11951203. PubMed Abstract  Publisher Full Text

Malmstrom J, Beck M, Schmidt A, Lange V, Deutsch EW, Aebersold R: Proteomewide cellular protein concentrations of the human pathogen Leptospira interrogans.
Nature 2009, 460(7256):762765. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Higham CF: Bifurcation analysis informs Bayesian inference in the Hes1 feedback loop.
BMC Syst Biol 2009, 3:12. PubMed Abstract  BioMed Central Full Text  PubMed Central Full Text

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. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Bassingthwaighte J, Hunter P, Noble D: The Cardiac Physiome: perspectives for the future.
Exp Physiol 2009, 94(5):597605. PubMed Abstract  Publisher Full Text  PubMed Central Full Text