Article Information

• PDF (552 kb)
• Full Text with Large Figures
• Supporting Material
• Cited by in Scopus (3)

PubMed

• Articles by Michael Cooling
• Articles by Edmund J. Crampin

Related Articles

• A Signal Transduction Pathway Model Prototype I: From Agonis...

  A Signal Transduction Pathway Model Prototype I: From Agonist to Cellular Endpoint Biophysical Journal, Volume 87, Issue 3, 1 September 2004, Pages 1406-1416 Thomas J. Lukas Abstract The postgenomic era is providing a wealth of information about the genes involved in many cellular processes. However, the ability to apply this information to understanding cellular signal transduction is limited by the lack of tools that quantitatively describe cellular signaling processes. The objective of the current studies is to provide a framework for modeling cellular signaling processes beginning at a plasma membrane receptor and ending with a measurable endpoint in the signaling process. Agonist-induced Ca mobilization coupled to down stream phosphorylation events was modeled using knowledge of in vitro and in vivo process parameters. The simulation process includes several modules that describe cellular processes involving receptor activation phosphoinositide metabolism, Ca-release, and activation of a calmodulin-dependent protein kinase. A Virtual Cell-based simulation was formulated using available literature data and compared to new and existing experimental results. The model provides a new approach to facilitate hypothesis-driven investigation and experimental design based upon simulation results. These investigations may be directed at the timing of multiple phosphorylation/dephosphorylation events affecting key enzymatic activities in the signaling pathway. Abstract | Full Text | PDF (189 kb)

• Calcium Signaling Pathways

  Calcium Signaling Pathways Biophysical Journal, Volume 94, Issue , 1 February 2008, Pages 150-157 Full Text | PDF (157 kb)

• PIP2 Hydrolysis and Calcium Release Are Required for Cytokin...

  PIP2 Hydrolysis and Calcium Release Are Required for Cytokinesis in Drosophila Spermatocytes Current Biology, Volume 15, Issue 15, 9 August 2005, Pages 1401-1406 Raymond Wong, Irene Hadjiyanni, Ho-Chun Wei, Gordon Polevoy, Rachel McBride, Kai-Ping Sem and Julie A. Brill Summary Summary | Full Text | PDF (359 kb)

• …more

Copyright © 2007 The Biophysical Society. All rights reserved.
Biophysical Journal, Volume 93, Issue 10, 3421-3433, 15 November 2007

doi:10.1529/biophysj.107.110031

Biophysical Theory and Modeling

Modeling Hypertrophic IP3 Transients in the Cardiac Myocyte

Michael Cooling^*, ,  , Peter Hunter^* and Edmund J. Crampin^*, †

^* Auckland Bioengineering Institute and the University of Auckland, New Zealand
^† Department of Engineering Science, University of Auckland, New Zealand

Address reprint requests to M. Cooling, Tel.: 64-9-373-7599, ext. 82165.

Pathological cardiac hypertrophy has been identified as a major risk factor leading to heart failure 1. An understanding of the biological pathways involved in cardiac hypertrophy may yield benefits in understanding the progression of, and the development of therapies for, cardiac disease.

A range of signal transduction pathways are implicated in hypertrophy, including the inositol 1,4,5-trisphosphate (IP3)-calcineurin pathway, the mitogen-activated protein kinase p42/p44, mitogen-activated protein kinase p38, jun N-terminal kinase, and the phosphoinositide 3-kinase-wingless-related mouse mammary tumor virus integration (PI3K-Wnt) pathways 2,3,4,5,6,7,8,9,10. In vivo, none has a monopoly on hypertrophic effects and many pathways interact 3, the details of which have yet to be fully elucidated experimentally.

The role of the IP3-calcineurin pathway in cardiac hypertrophy was discovered during the 1990s and, despite some initial controversy, calcineurin's role has now been accepted. Extracellular agonists stimulate G_q-family G-protein-coupled receptors (GPCRs). On stimulation, the G-protein α-subunit activates phospholipase C (PLC) 11,12, and causes hydrolysis of membrane-bound phosphatidyl inositol-4,5-bisphosphate (PIP2). This forms the secondary messenger IP3. IP3 production causes an increase in the intracellular calcium concentration ([Ca²⁺]_i), which in turn activates the phosphatase calcineurin. Calcineurin activation leads to changes in gene expression by facilitating translocation to the nucleus of cytosolic transcription factors of the nuclear factor of activated T cells (NFAT) family. There such factors bind to nuclear DNA cooperatively with transcription factors from the other hypertrophic pathways to facilitate transcription 13. NFAT's ability to act as a signal integrator for hypertrophy 14 in this manner is a key reason that the IP3-calcineurin pathway is significant.

These findings have motivated the use of existing immunosuppressive calcineurin inhibitors (such as Cyclosporin A and FK506) as experimental attenuators of hypertrophy. The results of that approach have been conflicting, with similar rodent models exhibiting no effect, attenuation, or even complete prevention of hypertrophy, and occasionally lethal side-effects from the immunosuppressors 15. Thus, a more detailed understanding of the pathways for hypertrophy is necessary for the development of practical therapeutic treatments.

The biological processes of the IP3 production system represent only a tiny proportion of the activities within a cell, yet even these are complex with many players and diverse interactions. For such systems a computational model is a useful tool 16, bringing together quantitative and qualitative data to allow interactive exploration of our understanding. Currently no mathematical model for the IP3-calcineurin pathway exists in the cardiac myocyte.

This work focuses on modeling the signal transduction pathway from extracellular agonist to the production of IP3 by the extracellular agonists endothelin-1 17 (ET-1) and angiotensin-II 18 (Ang-2). The responses of IP3 to the two agonists are different, despite the signal transduction pathway being almost identical. We construct the first mathematical model of the IP3 production and degradation signal transduction system in the mouse atrial cardiac myocyte, as stimulated by the hypertrophic agonist ET-1. Sensitivity analysis was undertaken to determine which parameters are significant in determining features of the IP3 transient. We investigate possible causes of the different IP3 response from stimulation with Ang-2, and show that the differences in system behavior are explainable in terms of receptor kinetics.

We constructed a mathematical model of the IP3 production system as stimulated by the α-adrenergic agonist ET-1 in the atrial cardiac myocyte. Since experimental processes are often time consuming and costly, it is important to assess from current knowledge what the most influential parameters are that determine properties of the system. We applied a global sensitivity analysis (the Morris Method) to determine the important parameters of the system by their influence on five characteristics of stimulated IP3 transients.

Using the ET-1 model as a baseline, we then refined the model to stimulation with the alternative agonist Ang-2. Ang-2 is a hormone produced via the angiotensin-renin system, and also produced locally on stretch of the cardiomyocyte 18,19. Ang-2 signals are received by two main isoforms of the Ang-2 receptor-Type 1 (AT-1) and Type 2 (AT-2). Cells from different animals contain different proportions of each although the proportion is approximately equal in rodents such as the rat 20. There is evidence linking AT-1 receptors to PLCβ, IP3 production, and calcium release via G_q family GPCRs 21,22,23. Transients resulting from stimulation with Ang-2 are much shorter-lived than those from ET-1, despite the fact that the same pathway is activated in both cases. We sought to determine whether the different transient behavior between the two agonists can be explained in our model by altering key parameters of the receptor alone. We found that the necessary changes to the model could be both motivated and understood by considering the results of the previously completed sensitivity analysis.

We developed a mathematical model of the significant biological processes relating to the transduction of extracellular ligand signals through G_q family receptors and subsequent IP3 production in the mouse atrial cardiac myocyte. A reaction scheme of the model is shown in Fig. 1. The model encompasses the extracellular ligand (L), the GPCRs (R), activation of phospholipase Cβ (P) by G_αGTP (G_t) subunits, and both the basal and stimulated hydrolysis of PIP2 to form IP3. Conceptually, the model can be divided into three modules:

Display large version of this figure

Figure 1

Reaction scheme of the IP3 production system. The extracellular ligand (L) binds to receptors (R), whether precoupled with G_αGDP (G_d, yielding R_lg) or not (R_l). Fully activated receptors (R_lg) release G_αGTP (G_t), which, along with calcium (Ca), stimulates PLCβ (P). In the unstimulated state, PLCβ-Ca²⁺ (P_c) hydrolyzes PIP2 to produce IP3 via reaction R14. When stimulated, PLCβ-Ca²⁺-G_αGTP (P_cg) hydrolyses PIP2 at a faster rate via reaction R15. Free IP3 is degraded via reaction R16.

The mathematical model was formulated using mass action kinetics, and implemented as a system of ordinary differential equations. For convenience, the model was encoded into the computer-readable protocol Cellular modeling Markup Language (CellML) 24. A complete list of the model equations, and a list of model parameters and initial conditions are provided in Table 3,Table 4, respectively. This signal transduction pathway involves interactions between membrane-bound and freely diffusing cytosolic species. Cytosolic species are represented by concentration (in μM) while membrane-bound species are given as area densities (number per unit of membrane area, in μm⁻²). The model is also available as CellML code in the online CellML Repository (http://www.cellml.org/models).

Extracellular ligand (L) binds to cell-surface receptors (R), shown by reversible reactions R1 and R4 in Fig. 1. For G_q family receptors, this binding causes a conformational change which replaces GDP with a GTP on an at-tached G-protein's α-subunit. This causes the G_α subunit to dissociate (reaction R6) and stimulate other proteins inside the cell. For a receptor to transduce the ligand signal across the membrane, it requires both an attached G-protein and a ligand-binding event. The model utilizes the receptor components from a nonexcitable cell model of Lukas 25 that uses the precoupled-receptor concept 26: should the G_α subunit bind first (as shown in reaction R2), then the receptor is known as precoupled (R_g), and exhibits a higher affinity for the ligand than if G_α had bound after the ligand (via reaction R3). Only the state where both ligand and G_α subunit are bound is considered an activated receptor (R_lg). Active receptors are targeted for phosphorylation and eventual invagination and recycling by the cellular machinery. Once phosphorylated, the receptor is no longer available for further signal transduction. For the purposes of the model, this is abstracted into a one-time phosphorylation step (reaction R5).

The model is designed to simulate experiments where the cell is initially unstimulated. At a predetermined stimulation timepoint (t_s), ligand is added, which is represented by a (smoothed) Heaviside step function to the desired stimulation concentration (L_s). The amount of ligand bound to receptors is assumed negligible compared to the total amount of ligand available, hence the concentration of free ligand in the model does not decrease as ligand binds to the receptors.

All reactions were modeled as kinetic fluxes following the Law of Mass Action. For example, the flux (J₁) for reaction R1

Endothelins bind to their receptors with high enough affinity for the reaction to be considered irreversible under physiological conditions 27, hence it seems likely that when the active G_α subunit dissociates during reaction R6, the ligand and receptor complex (R_l) remains together.

The active G_αGTP messenger (G_t) released from stimulated receptors binds to the enzyme PLCβ (P), shown in Fig. 1 as the reversible reactions R9 and R10. This binding can be considered independently from the binding of the costimulatory calcium ions (Ca) shown in reversible reactions R8 and R11. Free G_αGTP degrades to G_αGDP via a self-GTPase activity of the α-subunit, removing its ability to stimulate PLCβ, represented as reaction R7. This activity in PLCβ-bound G_αGTP and dissociation from PLCβ is assumed for modeling purposes to occur in the same step and is shown by reactions R12 and R13 depending on whether Ca²⁺ is also bound or not. When bound to PLCβ, this step is assumed to occur at the same rate, irrespective of the additional presence or absence of bound calcium (from P_cg or P_g, respectively).

The formulation of the reactions R7-R13 follows the same form described in GPCR cycling above, although for some reactions an additional conversion factor is required: Ca²⁺ is a cytosolic species, measured in μM whereas the membrane-associated fluxes such as J₁ above are membrane density fluxes (μm⁻² s⁻¹). To convert to concentration fluxes, it is convenient to define a conversion factor

(1)

A number of parameters were taken from the existing model by Lukas 25, where appropriate. For the binding of Ca²⁺ to PLCβ-G_αGTP (Ca binding to P_g in reaction R11) only the dissociation constant (K_d,11) was known. We have made the assumption that the forward rate constant for reaction R11 is twice that of the similar reaction R8, due to the binding of G_αGTP. From this assumption the reverse rate constant for reaction R11 was calculated.

The resting level of calcium in the cardiac myocyte was assumed to be 100 nM 28.

The species PLCβ-Ca²⁺ (P_c) and PLCβ-Ca²⁺-G_αGTP (P_cg) act as hydrolyzers, converting the substrate PIP2 to IP3, shown by reactions R14 and R15 in Fig. 1, respectively. These catalytic reactions were modeled using Michaelis-Menten kinetics according to the quasi-steady-state approximation, their equations therefore being of the form (J₁₄ is shown as an example)

IP3 binds to IP3 receptors (IP3Rs) in the cell, which are also regulated by calcium. It was assumed that the amounts of both IP3 and calcium bound are negligible compared to their cytosolic concentrations, and therefore IP3Rs do not effect the concentrations of those small species. IP3R activation causes a minor cytosolic calcium spike from intracellular stores, followed by a much larger sustained calcium influx from outside the cell—this latter phenomenon being known as capacitative calcium entry 29. In the protocol of interest to this work, the extracellular entry does not occur, and the minor initial spike is assumed to have little effect on the model calcium. The IP3Rs, having no effect on the rest of the model, are therefore not included in this formulation.

The hydrolysis of PIP2 also produces the membrane-bound diacylglycerol species, but that is not considered in this model.

PIP2 is manufactured by the cell, and resupplied after use to the plasma membrane via the breakdown of inositol phosphates 30. Since there is no evidence that the level of PIP2 is rate-limiting 31, in the model PIP2 is fixed at its initial value.

A degradation reaction (R16) encompasses the dephosphorylation or phosphorylation of IP3 to form inositol bisphosphate and inositol tetraphosphate, respectively. For the purposes of this model, the metabolic recycling of IP3 from those products is not considered. The basal level of IP3 is therefore a balance between the basal production and the degradation rates (reactions R14 and R16, respectively).

A published report 32 of a ∼30 nM peak in IP3 on stimulation by 100 nM of ET-1, together with estimates of the peak-to-basal ratio of ∼2:1 33,34, suggest that the basal concentration of IP3 in the mammalian myocyte is ∼15 nM. This is the estimate used in this model. The kinetics of the basal IP3 production reaction (via k_f,14) were adjusted to provide this value. The balancing degradation rate was known with more confidence and kept as in the literature 25.

The model was refined to simulate the response to ET-1 stimulation (see Results for a comparison of model behavior to experimental data). We then sought to determine which parameters controlled the IP3 transient produced by the system by undertaking a global sensitivity analysis following the Morris method 35. The Morris method is a screening algorithm: it ranks the parameters of the model, by their average effect on a particular model simulation output, over a given range of parameter values. The global nature of the analysis is applicable given the system's high nonlinearity 36 and the method has been shown to be as efficient as commonly-used variance based techniques in detecting factors of low sensitivity. It is also computationally cheaper than such techniques 37 as the required number of model evaluations varies linearly with the number of model factors.

The method yields a score for each parameter of the model estimating the partial differential of an objective function with respect to that parameter. We chose five objective functions to analyze the model by, to capture the behavior of the resulting IP3 transient (these measures are illustrated graphically on a representative IP3 transient in Fig. 4):

To determine the sensitivity of these measures to the model parameters, we define for each parameter a range of values over which it can vary. These may vary widely, whether by virtue of the very different quantities which they represent, or simply by different physiological constraints. The amount of variation of a parameter inside a particular range is normalized to a fraction of that range, designated δ. The parameter ranges defined for our analysis are listed in Table S1 in the Supplementary Material . Parameter ranges were normalized and parameters allowed to take one of 16 equally spaced values (Morris parameter p=16). For each parameter, changes in each of the objective functions were measured for a δ change over the range defined, where δ is 38

This yields a sliding window of 53.3% of each parameter's range. Parameters therefore begin each of a series of analysis trajectories at one of 16 values within their defined range, and at each step of each trajectory, one parameter is varied by 53.3% of that range, and the subsequent effect on each the objective functions are measured.

Parameter change effects are divided by δ, to give an estimate for the parameter's effect on each objective measure, over the parameter's entire defined range. We ran two sets of 8000 trajectories as described by Morris 35 for each parameter over the five objective functions (Morris parameter r=8000), as we found this gave repeatable clustering and ordering (unless parameters were very similar in significance) results between the two sets.

For each test, the model was run to achieve steady state before agonist stimulus was applied. Up to a further 130,000s was simulated (∼36h) to allow IP3 to peak and resume the baseline level again for all parameter sets tested.

We then combined both sets of tests to give an overall population for each parameter of 16,000 sensitivity measures, and used the Campolongo extension 38 of the method to calculate the mean absolute effect (μ*) as a direction-independent measure of sensitivity of the defined normalized ranges.

Tables S2–S6 in the Supplementary Material show the parameters ranked by their effects thus calculated, for each of the five objective functions.

We sought to determine whether the different IP3 transient produced by Ang-2 (as opposed to ET-1) stimulation could be explained by receptor-related parameters alone. A number of model parameters were therefore adjusted to match known published kinetics for the Ang-2 receptors. It is known that the rodent Ang-2 receptor has a larger dissociation constant than the ET-1 receptor: ∼1.5nM in rat and ferret, and that the receptor density is much lower than that for ET-1, being approximately one-fifth (see Table 25 in 28). These literature-motivated parameter changes alone were not sufficient to explain the differences in the resulting IP3 transient. We therefore investigated the most significant remaining receptor-related parameters.

We used the results of the sensitivity analysis to identify the 10 most significant parameters with respect to the IP3 transient (see Table 1). Two of these (k_f,4 and k_f,5) relate to kinetics associated with the receptors themselves, and hence are likely to be different between ET-1 and Ang-2 receptors.

Experimental research by Abdellatif et al. 39 measured accumulation of IP3 in cultured rat neonatal ventricular cardiac myocytes on stimulation of 100μM Ang-2, using ion-exchange chromatography. During fitting to this data, we determined predicted new values for those significant receptor parameters using the Levenburg-Marquardt algorithm (implemented as the “levmar” nonlinear least-squares fitting package obtained from http://www.ics.forth.gr/(lourakis/levmar). Multiple runs were conducted in an attempt to evade local minima, and a close fit with the experimental data was obtained. Varying additional parameters did not significantly increase the goodness-of-fit.

Model parameters were refined to the mouse atrial myocyte by fitting the IP3 transient resulting from stimulation at various ET-1 strengths to experimental data from Jiang et al. 40. Based on the available data, we performed the fit based on IP3 transient observations. To our knowledge data on the ET-1 stimulated responses of proteins upstream of IP3 do not currently exist. Model parameters were adjusted until the output from the model closely matched the experimental observations, to match the mouse atrial myocyte on stimulation by the specific extracellular ligand ET-1.

To achieve the fit, several parameters’ values had to be changed from those derived from other sources. These changes are presented in Table 4.

A comparison of an IP3 transient from the model with experimentally determined IP3 time-course data from Jiang et al. on 100μM of ET-1 stimulation is shown in Fig. 2. Simulation results were generated from the model defined in Table 3, with the parameters defined in Table 4, and L_s=0.100μM, t_s=30.0s, to match the experimental protocol. The numerical solver used was an implicit Runge-Kutta method based on Radau quadrature, as described by Hairer and Wanner 41. Absolute and relative error tolerances were set at 10⁻⁹. As can be seen, the model behavior closely matches experimental results.

Display large version of this figure

Figure 2

IP3 transient-based objective measures. Graphical depiction of the measures used to define the objective functions for the sensitivity analysis. Tau is the time point at which the transient dips beneath 1/e×amplitude.

A further comparison with IP3 transient data is shown in Fig. 3. Jiang et al. measured the IP3 concentration in cells 30min after stimulation with various concentrations of ligand. To simulate this, the model's ligand stimulation strength (L_s) was varied at regular intervals from 1×10⁻⁴μM to 1.0μM, and the concentration of IP3 measured 1800s (30min) after ligand addition. Figure 2 and Figure 3 demonstrate that appropriate IP3 transient behavior is produced by the model.

Display large version of this figure

Figure 3

IP3 transient curve on ET-1 stimulation. The simulated IP3 transient closely matches experimental observations 40 on application of 100 nM ET-1. Simulations were performed with the equations as in Table 3 and the parameters and initial conditions as in Table 4. Additionally, t_s=30s, and L_s=0.100μM.

Fig. 4

Display large version of this figure

Figure 4

IP3 dose-response curve on ET-1 stimulation. The model also closely matches experimental observations 40 for the concentration of IP3, after 30min, on stimulation by various concentrations of ET-1. Simulations were performed with the equations as in Table 3 and the parameters and initial conditions as in Table 4. Additionally, t_s=100s, and the model was run until (t−t_s)=1800s (30min of stimulation) for each of L_s=1×10⁻⁴, 5×10⁻⁴, 1×10⁻³, 5×10⁻³, 1×10⁻², 5×10², 1×10⁻¹, and 5×10⁻¹μM.

The most significant parameters in terms of impact on the IP3 transient are listed in Table 1. Parameters were sorted by summing their proportional significance for each objective measure. A complete list of all model parameters and their significance scores can be found in Table S7 in the Supplementary Material .

The single most important parameter is the receptor phosphorylation rate constant (k_f,5). The quantity k_f,5 governs the reaction that terminates the signal by reducing the total available receptor density, and thus reducing the possibility of sarcolemmal signal transduction in response to extracellular hormonal stimulation. It is of critical importance for setting the Time-to-Baseline, as without a phosphorylation step the system would simply equilibrate based on the stimulation concentration of ligand and the transient would be constant. Decreasing the phosphorylation rate via the k_f,5 parameter increases the tail as shown in Fig. 5. As can also be seen in that figure, a higher k_f,5 reduces the Time-to-Peak. While lowering k_f,5 increases the length of the transient's tail it also increases Tau. This is to be expected, however the increase in Tau is proportionally less than the increase in the tail length. Hence, over the range tested, the Tau-to-Tail Ratio reduces as k_f,5 increases. These three effects make the phosphorylation rate of the receptors crucial to determining the transient's behavior.

Display large version of this figure

Figure 5

Decreased k_f,5 results in a higher peak and a longer time-to-baseline from peak. Simulations were performed with the equations as in Table 3 and the parameters and initial conditions as in Table 4. Additionally, t_s=100s, L_s=0.100μM, and k_f,5 was varied as depicted in the legend.

The following two parameters—the ligand strength (L_s) and the rate constant for the binding of ligand to precoupled receptor (k_f,4)—have virtually identical effects across all objective function measures (see Table S7 in the Supplementary Material for quantitative details). We consider each in turn.

It might be expected that the applied ligand concentration (L_s) is significant in determining the response. L_s is significant largely due to its effect on the Time-to-Peak, which is demonstrated in Fig. 6. For stimulation with lower ligand concentration, the transient flattens and takes longer to achieve its maximum. This is due to the effect that ligand concentration has on the phosphorylation of available receptor. With higher ligand concentration, more active receptor is immediately available for phosphorylation, which turns off the signal faster. For low ligand strength, less IP3 signal amplitude is achieved, but the receptor pool is not consumed as quickly. This can be seen from a plot of the percentage of phosphorylated receptor over time, as shown in Fig. 7 for the highest and lowest ligand strengths in Fig. 6.

Display large version of this figure

Figure 6

Decreasing ligand strength strongly effects Time-to-Peak. Lower ligand strength produces a slower rate of active receptor production, thus the transient takes more time to achieve its maximum. Simulations were performed with the equations as in Table 3 and the parameters and initial conditions as in Table 4. Additionally, t_s=100s and L_s was varied as depicted in the legend.

Display large version of this figure

Figure 7

Higher ligand strength more readily increases phosphorylated receptor density. Although the long-term gradients for phosphorylated receptor density are similar for both high and low concentrations of ligand strength, a higher strength yields a much larger initial gain. Simulations were performed with the equations as in Table 3 and the parameters and initial conditions as in Table 4. Additionally, t_s=100s and L_s was varied as depicted in the legend.

The forward rate constant for the binding of ligand to precoupled receptor (k_f,4) also heavily influences the formation of activated receptors. This parameter has a large effect on the Time-to-Peak. Examination of the changes in the membrane density of various species reveals that the IP3 transient closely follows the shape of the density curve on the active receptors (the product of reaction R4) (see Fig. 8). The parameter k_f,4's effect on Time-to-Peak can thus be explained by its influence on the production of active receptor, which in turn leads to the peak for IP3.

Display large version of this figure

Figure 8

The IP3 transient follows the active receptor transient. Simulations were performed with the equations as in Table 3 and the parameters and initial conditions as in Table 4. Additionally, t_s=100s and L_s=0.100μM.

Both parameters also significantly influence the Time-to-Peak and the Tau-to-Tail Ratio. A higher peak due to high active receptor formulation (as increased by k_f,4, L_s and also the parameter G_d) leads sooner to a more rapidly decreasing transient, as receptors are more readily available to be phosphorylated faster. This gives a sharper peak earlier, with a longer post-peak tail and therefore a lower Tau-to-Tail Ratio, compared to a flatter transient with a longer time to peak, a shorter tail, and higher Tau-to-Tail Ratio when those parameters are lower. This is illustrated for k_f,4 specifically in Fig. 9. The same effect occurs when L_s or G_d are decreased (not shown).

Display large version of this figure

Figure 9

Higher k_f,4 reduces the Tau-to-Tail Ratio. Increased k_f,4 gives sharper peaks, but longer tails from the peak, and a lower Tau-to-Tail Ratio compared to flatter transients with a higher ratio when k_f,4 is decreased. Simulations were performed with the equations as in Table 3 and the parameters and initial conditions as in Table 4. Additionally, t_s=100s, L_s=0.100μM, and k_f,4 was varied as depicted in the legend.

The next most important parameter is the forward rate constant (k_f,16) for the degradation of IP3, which is the main consumer of IP3 in this system. The ratio of sarcolemmal surface area to cytosolic volume (R_pc) is also important as it determines the cytosolic flux in IP3 concentration resulting from a given surface area on the plasma membrane (a smaller cytosolic volume would yield a higher concentration change for the same IP3-producing sarcolemmal surface area). Hence this parameter has a direct effect on the concentration of IP3 produced for a given stimulus. Both these parameters mainly effect the IP3 Baseline Concentration and the Peak Concentration—one parameter affecting via IP3 degradation and the other affecting production, respectively.

The top five parameters together represent the most important determinants for each of the five objective functions—k_f,5 for Time-to-Baseline, L_s and k_f,4 both equally for each of Time-to-Peak and Tau-to-Tail Ratio, and both k_f,16 and R_pc for each of Baseline and Peak Concentration.

The value G_d represents the densities of a species required to form that mobile messenger—the amount of free G_αGDP. G_d's effects on Tau-to-Tail Ratio and Time-to-Peak is similar to the effect of k_f,4 or L_s, as has already been discussed, being a multiplicative factor in the forward rates of both reactions R2 and R3, which form active receptors.

K_d,2 has similar but less pronounced effects to G_d, which is due to it being a determinant of the rate of precoupled receptor formation, and therefore influences the rate of production of active receptors.

The basal IP3 production forward rate constant (k_f,14) is also important, largely for its effect on the Baseline Concentration. This is expected since the basal IP3 concentration is a balance of reaction R14 and the degradation rate as previously discussed.

The next two parameters of importance are the forward and reverse rate constants for the binding of calcium to PLCβ. The PLCβ-Ca²⁺ so produced is the only enzyme capable of producing IP3 in this system without ligand stimulation. Hence their effect on Baseline and Peak Concentrations, while less than that of other parameters such as k_f,16 and R_pc, is notable.

Together these are the most significant control parameters of the IP3 model. They are also, therefore, the key parameters to be determined experimentally for a species of interest, as they have the strongest effect on the time course of the IP3 transient after stimulation of the G_q protein-coupled receptors.

The literature-derived and refinement-fitted parameter changes required to adjust the model (as discussed in Methods) for Ang-2 cell stimulation are shown in Table 2. A graph of the experimental observations together with the fitted IP3 transient is shown in Fig. 10.

Display large version of this figure

Figure 10

Simulated fit to the observations of Abdellatif et al. 39. The simulated transient matches the observations 39 for stimulation with 100 nM of Ang-2. Simulations were performed with the equations as in Table 3 and the parameters and initial conditions as in Table 4, and the Ang-2 values listed in Table 2. Additionally, L_s=100μM and t_s=100s.

Table 3 IP3 model equations

	if Ls; if otherwise.
	kr,1=kf,1×Kd,1
	J1=kf,1×R×L−kr,1×Rl
	kr,2=kf,2×Kd,2
	J3=kf,3×Rl×Gd−kr,3×Rlg
	kr,4=kf,4×Kd,4
	J4=kf,4×L×Rg−kr,4×Rlg
	J5=kf,5×Rlg
	J6=kf,6×Rlg
	J7=kf,7×Gt
	J9=kf, 9×P×Gt−kr,9×Pg
	J8=kf, 8×P×Ca−kr,8×Pc
	J10=kf,10×Pc×Gt−kr,10×Pcg
	kr,11=kf,11×Kd,11
	J11=kf,11×Pg×Ca−kr,11×Pcg
	J12=kf,12×Pcg
	J13=kf,13×Pg
	J16=kf,16×IP3

Three parameters were adjusted during the refinement. An increase in k_f,5 is needed to reduce the length of the transient to the shorter Ang-2 response. A minor adjustment to k_f,4 (less than an order of magnitude) was also made. A third parameter was also modified—for forward rate constant for the dissociation of the messenger G_αGTP from active receptors (k_f,6). This modification was necessary due to the increased k_f,5, which results not only in a shorter tail, but a lower peak. By increasing k_f,6, more signal can be transferred to the PLCβ module via the messenger G_αGTP, raising the IP3 transient peak without greatly increasing the tail. This combination of parameter changes, in concert with those derived from the literature, was sufficient to explain the experimental data.

Thus, once we allow for different receptor densities and known kinetic differences, the difference between the ET-1- and Ang-2-induced IP3 transients can be explained solely in terms of additional receptor kinetic changes in the phosphorylation rate (via rate constant k_f,5), the binding of extracellular ligand to the precoupled receptor (k_f,4) and the release of active G-protein (k_f,6).

We have presented the first model of the IP3 production system in the atrial cardiac myocyte. We used this model to 1), assess the key drivers of the IP3 transient system; and 2), hypothesize the kinetic mechanism that allows two different receptors to provoke different behaviors in the same downstream pathway.

The sensitivity analysis showed that several parameters that we expected to be important were in fact significant, increasing our confidence in the model. It also rated some parameters unexpectedly high—indications which have led to further insights on the functioning of the IP3 production system.

That the phosphorylation rate of the receptor (via k_f,5) is of high importance is expected. There are two other switch-off points for the stimulated IP3 transient—the G_αGTP self-hydrolysis to G_αGDP when mobile (via k_f,7), or when bound to PLCβ-Ca²⁺ (via k_f,12 and k_f,13). The fact that the phosphorylation rate of the receptor is more important than these highlights the importance of the interface of the cell to the surroundings above the internal mechanisms, at least in this system.

It is also intuitive that the ligand strength (L_s) is important. What is perhaps less obvious but shown by the sensitivity analysis is that this effect seems to be due largely to the ligand's effect on precoupled receptor density as opposed to its effect on noncoupled receptor density. The impact of L_s is identical to that of the forward rate constant for the binding of that ligand to precoupled receptor. This is explainable by both being multiplicative factors in that reaction (R4). It is such reaction that is important for the formation of active receptors (R_lg) on stimulation by ligand. The noncoupled reaction (R1) is less important on the stimulation of ligand as a further reaction must then be carried out (reaction R3) for the G_αGDP to bind to the ligand-receptor complex before the receptor is active. In the precoupled state, the receptor becomes active much more readily and determines the responsiveness of the system far more significantly.

It seems likely that the effect on precoupled receptor density is also why the amount of G_αGDP and the dissociation constant for receptor precoupling (K_d,2) are significant at numbers six and seven in the top 10—which may have at first seemed puzzling. G_αGDP is a reactant for the formation of precoupled receptors (R_g), which in turn are reactants for that important reaction R4.

That the degradation rate of IP3 (via k_f,16) and geometric ratio of membrane area-to-cytosolic volume (R_pc) are important is also expected when their role in the degradation and production of IP3, respectively, is considered, as was noted above. That the parameters (k_f,8 and k_r,8) related to the formation of the basal PIP2-hydrolyser (P_c) and the subsequent hydrolysis reaction (k_f,14) are important is not surprising, given that the Baseline Concentration was one of the objective functions. The availability of P_c influences the resulting stimulated transient also, as mass of PLCβ is conserved when forming the stimulated hydrolyser P_cg, and k_f,14 influences the difference between P_c activity and P_cg activity, having implications for the tail of the IP3 transient. Hence overall these are among the more significant parameters of the system.

Of the top seven parameters, five relate to the receptor kinetics. This highlights again the importance of the cell/environment interface in this system and indicates how a change in receptor density and kinetics from ET-1 to Ang-2 can alter the resulting IP3 transient to the extent observed experimentally.

Taken together, only parameters relating to the receptor kinetics and density have been altered to yield the fits to experimental data. Hence the model is consistent with the idea that a shorter IP3 transient from Ang-2 compared to ET-1 can be produced simply by a change in the type and number of receptors without the need for any other mechanisms.

The change in phosphorylation rate of active receptor (via constant k_f,5) is consistent with the observation that the receptor phosphorylation rates differ even between subtypes of the same receptor 42. As discussed above, this phosphorylation rate is the most important factor in determining the time for the transient to reach the baseline level again after stimulation. Therefore in determining the possible differences in receptor kinetics between ET-1 and Ang-2 binding receptors, an increase in the phosphorylation rate (to reduce the transient's time length) for Ang-2 receptors seems likely.

The increase required in the dissociation rate of active G-protein subunits from the receptor (k_f,6) suggests that, in cardiomyocytes, along with an increased phosphorylation rate (compared with endothelin receptors), the AT-1 receptor may also have an increased rate of G_αGTP dissociation from the activated receptor compared to that of the ET-1A receptor. This prediction, and the increase in the rate of binding of ligand to precoupled receptor (k_f,4), could be investigated experimentally.

We also refined the model to a second Ang-2 dataset from rat neonatal ventricular myocytes, that of Sadoshima and Izumo 43, which measured IP3 transients at two different concentrations of Ang-2. The experimental observations were too coarsely grained for our purposes to be confident about the peak of the transient, but best fits were achieved with similar values for k_f,5 and k_f,6 (0.0421s⁻¹ compared with 0.0622s⁻¹ and 14.9s⁻¹ compared with 22.2s⁻¹, respectively), increasing our confidence in the roles of the receptor phosphorylation and active G-protein dissociation rates.

The importance of the receptor phosphorylation rate, as evidenced both by the sensitivity analysis and its role in the fitting of the Ang-2 dataset, highlights the value of experimental measurement of this parameter for receptors of interest. Often experiments focus on determining the binding constants to certain ligands. As we have shown, for signal transduction the resulting intracellular signal is heavily influenced by the temporal dynamics of receptor activation, which is a balance between the binding of the ligand and the molecul