<title>Introduction</title> <p>The G-protein-coupled signal transduction system integrates a wide range of intercellular signals and actuates downstream pathways. G-protein-coupled receptors (GPCRs) are composed of seven α-helices that span the plasma membrane, an extracellular domain that is activated by an agonist and an intracellular domain that binds a guanine nucleotide heterotrimer made up of different α, β, and γ subunit isoforms. This receptor system accounts for 40–50% of modern medicinal drug targets but only 10% of the known receptors are targeted by drugs ##REF##14625380##[1]##. Though the system is physiologically and pharmacologically important, the mechanism by which the system integrates multiple signals is not well understood ##REF##12790797##[2]##.</p> <p>We address the G-protein-mediated route to calcium release in RAW264.7 cells. When activated by a specific ligand, the G protein heterotrimer dissociates to free Gα-GTP and Gβγ. Specific Gα and Gβγ isoforms are able to bind specific isoforms of phospholipase C β (PLCβ) and catalyze the synthesis of inositol (1,4,5)-triphosphate (IP3) and diacylglycerol (DAG) from phosphatidylinositol (4,5)-bisphosphate (PIP2) ##REF##2830256##[3]##,##REF##8389480##[4]##. In addition to its catalytic activity, PLCβ acts as a GTPase for Gα-GTP ##REF##10449728##[5]##. IP3 binds to specific receptor-channels on the membrane of the ER to release Ca<sup>2+</sup> into the cytosol ##REF##15189149##[6]##. DAG and Ca<sup>2+</sup> bind to and activate protein kinase C (PKC) which may phosphorylate and inactivate specific PLCβ isoforms ##REF##12954613##[7]##. G protein receptor kinase (GRK) is activated once it is phosphorylated by PKC ##REF##14499340##[8]## and is localized to the plasma membrane by Gβγ ##REF##7744811##[9]##. Though phosphorylation has not been shown to be necessary for GRK activation, we have assumed so in our model because phosphorylation by PKC may release the inhibition of GRK2 by being bound to calmodulin ##REF##14499340##[8]##. Activated GRK can then phosphorylate specific GPCRs which leads to receptor inactivation—perhaps directly or by arrestin activity ##REF##14499340##[8]##. In this complex signal transduction network, Gα and Gβγ subunits have different patterns of specificity for PLCβ isoforms and calcium is an important cofactor in several important feedback loops ##UREF##0##[10]##.</p> <p>The two extracellular signaling ligands we consider here are C5a and UDP. The small peptide C5a is a potent anaphylotoxin and a strong chemoattractant for many immune system components ##REF##15638737##[11]##. The calcium response due to stimulation by C5a is predominantly coupled through Gαi-linked heterotrimers. Macrophage cells and their precursors, monocytes, express several receptors that are specific to extracellular nucleotides and it has been shown that the P2Y6 receptor, which is sensitive to UDP, regulates the production and secretion of the chemokine interleukin 8 (IL-8) in monocytes ##REF##11349132##[12]##. The UDP response is mediated by Gαq-linked heterotrimers, but other receptors in the P2Y family may respond to UDP and couple the signal through other G protein isoforms ##REF##11390975##[13]##.</p> <p>Four recent models have sought to explore various aspects of the G protein coupled signal transduction system in detail. Lukas et al. compare measured calcium response over a range of bradykinin doses to their model predictions ##REF##15345523##[14]##. Mishra and Bhalla built a model to investigate the role of IP4 as a signal coincidence detector in the GPCR pathway ##REF##12202356##[15]##. The model by Lemon et al. predicts the calcium response to UTP stimulation and is the closest in focus to our model ##REF##12782119##[16]##. A recent model of calcium dynamics in RAW cells has been proposed that is quite similar to this model, but does not deal with crosstalk between receptors or formal statistical uncertainty in model predictions ##REF##17483189##[17]##,##REF##17483174##[18]##.</p> <p>Several hypotheses for the mechanism of crosstalk and synergy among GPCR-mediated pathways have been proposed. Crosstalk among GPCR-mediated pathways is important both physiologically and pharmaceutically. Quitterer et al. propose that crosstalk is mediated by Gβγ exchange between Gαi-coupled and Gαq-coupled receptors ##REF##10485876##[19]##. Zhu et al. speculated that PLC is under either conditional or dual regulation of Gβγ and Gα ##REF##8610126##[20]##. Though these hypothetical mechanisms for crosstalk among G protein coupled receptor systems are conceptually plausible we have not found these or any other of the many competing hypothetical mechanisms tested in the context of a quantitative mathematical model ##REF##12790797##[2]##.</p> <p>In this paper Bayesian statistical inference is used to provide a rigorous connection between the mathematical model derived from mass-action kinetics, prior information from in-vitro biochemical studies and heterogeneous experimental data. The prior distribution over the parameters represents our uncertainty before observing a set of experimental data. A broad, high variance, prior distribution means we are quite uncertain and a concentrated, low variance, prior means we are more certain about the parameter a priori. The objective of our inference is the posterior distribution over the parameters because it is an informed estimate of both the value of the parameter and the uncertainty in the parameter value. The posterior distribution over the parameters is then used as a tool for experiment design to estimate the model-based posterior distribution over observable quantities such as the cytosolic calcium concentration and to drive the design of new experiments. This statistical approach is possible in a model of this size because of the abundance and quality of the data collected for this study.</p>
<title>Materials and Methods</title> <p>The model equations are given in ##SUPPL##7##Figure S7##. The initial conditions and parameter values are in ##SUPPL##11##Table S1## and ##SUPPL##12##Table S2##, respectively. All the data used in this work and a stand-alone implementation of the model is provided at <ext-link ext-link-type="uri" xlink:href="http://genomics.lbl.gov/supplemental/flaherty-gpcr/">http://genomics.lbl.gov/supplemental/flaherty-gpcr/</ext-link>. The model was simulated using CVODE ##UREF##1##[42]## and the GNU Scientific Library. Further details on materials and methods are available in ##SUPPL##0##Dataset S1##.</p> <title>Experimental Methods</title> <p>Intracellular free calcium in cultured adherent RAW264.7 cells was measured in a 96-well plate format using the Ca<sup>2+</sup>-sensitive fluorescent dye Fura-2 ##REF##2538708##[43]##,##REF##3838314##[44]##. A Molecular Devices FLEXstation scanning fluorometer was used to measure fluorescence using a bottom read of a 96-well plate. Each well was sampled approximately every 4 seconds. The measurement protocol is described in AfCS experimental protocol ID #PP00000211 (available from <ext-link ext-link-type="uri" xlink:href="http://www.signaling-gateway.org">http://www.signaling-gateway.org</ext-link>). The parameters in ligand concentration model were estimated using FITC solution in the FLEXstation scanning fluorometer as described in Molecular Devices Maxline Application Note #45 and in ##SUPPL##10##Protocol S1## (see also ##SUPPL##5##Figure S5##).</p> <title>Statistical Inference</title> <p>Twenty of the 84 parameters were chosen to be estimated from data based on relevance to the experimental hypothesis. Only those parameters that related to the knockdown experiments in the dataset were estimated and are denoted with a star in ##SUPPL##12##Table S2##. We used data to estimate only the two forward rate constants in the enzymatic mass-action equations because the forward and reverse rate constants for a given reaction will be highly correlated in the posterior distribution making estimation by Markov chain methods computationally expensive. An analysis of the sensitivity of the model to each parameter is shown in ##SUPPL##9##Figure S9##.</p> <p>For each estimated parameter we constructed an independent Gaussian prior on a log scale with a mean chosen based on relevant literature and a standard deviation of 0.25. We found that this prior variance was sufficiently permissive to allow exploration of the space while still constraining the rates to be physically reasonable. The prior distribution over the parameters allows the incorporation of both soft and hard constraints in the parameter estimates. Parameter sets with zero measure are not permitted in the posterior distribution and parameter sets with small measure must be assigned a large likelihood in order to have a large posterior probability.</p> <p>The likelihood is a function of the parameters (<italic>θ</italic>) and links the prior distribution with the posterior distribution under Bayes rulewhere <italic>y</italic> denotes the observed data.</p> <p>In our model, the likelihood function is a Gaussian distribution according to the non-linear regression equation <italic>y</italic> = <italic>f</italic>(<italic>θ</italic>)+<italic>ε</italic>, <italic>ε</italic>∼<italic>N</italic>(0,<italic>σ</italic>
<sup>2</sup>), where <italic>f</italic>(<italic>θ</italic>) is the deterministic model prediction. The posterior distribution is of interest because it informs us as to the most probable setting of the parameters as well as the uncertainty in the values.</p> <p>The Metropolis-Hastings algorithm ##UREF##2##[45]## was used to estimate the posterior density of the parameters Pr(<italic>θ</italic>|<italic>y</italic>). Three independent chains were simulated from different initial parameter values (see ##SUPPL##1##Figure S1##). To assess convergence of the posterior distribution estimate, we used the Gelman-Rubin potential scale reduction factor (PSRF) ##UREF##3##[46]##. The multivariate PSRF is 2.44 and 95% of the individual PSRFs were less than 1.5. A PSRF value of one indicates that the distribution has converged and values near one are close to converged.</p> <p>Posterior prediction confidence intervals were constructed using the percentiles from the predictive distribution approximated with 2000 Monte Carlo samples from Pr(<italic>y</italic>
<sub>new</sub>|<italic>θ<sub>i</sub></italic>) at each of 100 simple random samples from Pr(<italic>θ</italic>|<italic>y</italic>) obtained fromwhere Pr(<italic>y</italic>
<sub>new</sub>|<italic>θ<sub>i</sub></italic>)∼<italic>N</italic>(<italic>f</italic>(<italic>θ</italic>),<italic>s</italic>
<sup>2</sup>) and <italic>s</italic>
<sup>2</sup> is the pooled variance estimate, which is computed as an average of the variances of all the time points in each of the 29 wild-type experiments. These average variances were weighted by the number of technical replicates in each experiment and then averaged to yield the estimate <italic>s</italic>
<sup>2</sup>. A small factor of 1 nM<sup>2</sup> was added to each variance estimate to bound variance estimates away from zero.</p>
<title>Results</title> <p>There are two main features of the structure of our model, shown in ##FIG##0##Figure 1##, which contribute to crosstalk in the system and produce the key dynamical features in the calcium response: isoform specificity and calcium-dependent feedback. As we will show, by including multiple isoforms of PLCβ and Gα as well as the negative feedback mediated by PKC, GRK and the IP3 receptor itself, we are able to predict the synergistic interaction between C5a and UDP observed in the experimental data.</p> <p>Our representation of the G-protein-coupled signal transduction system includes C5a and P2Y6 receptors, Gαi2, Gαq, Gβγ, PLCβ3, PLCβ4, PIP2, DAG, IP3, PKC, GRK2, calcium buffer, a Na<sup>2+</sup>/Ca<sup>2+</sup> exchanger, a sarco(endo)plasmic reticulum Ca<sup>2+</sup>-ATPase (SERCA) pump, IP3 receptors and RGS. The model is composed of 53 coupled ordinary differential equations with 84 parameters and 24 non-zero initial conditions. The complete model equations are shown in ##SUPPL##7##Figure S7## and a more detailed model diagram is shown in ##SUPPL##6##Figure S6##. The parameters and initial conditions are in ##SUPPL##12##Table S2## and ##SUPPL##11##Table S1##, respectively. Where available, we have relied on <italic>in-vitro</italic> or <italic>in-vivo</italic> biochemical experiments for the reactions and parameter values (see Supporting Information). In cases where the biochemical parameter values were not known, we chose physically reasonable values. Twenty of the 84 parameters most relevant to the knock-down and wild-type data were estimated from cytosolic calcium measurements as described in the <xref ref-type="sec" rid="s4">Methods</xref> section. Most reactions were assumed to be governed by mass-action kinetics, but for a few proteins—such as RGS—the mechanism of regulation is not known in enough detail and we have approximated with Michaelis-Menten kinetics or a phenomenological function.</p> <p>We briefly discuss the reactions involving the Na<sup>2+</sup>/Ca<sup>2+</sup> exchanger, SERCA pump, IP3 receptors, RGS and calcium buffer because they are important for the faithful representation of the system in our model. Regulators of G protein signaling (RGS) are GTPase proteins that down-regulate the extent of signaling ##REF##9462506##[21]##; RGS2 at least is expressed in RAW264.7 macrophage cells and therefore an RGS activity is included in our model. The mechanism of activation of RGS2 as it relates to Gαi and Gαq signaling is not entirely known and is difficult to assess because antibodies that specifically recognize RGS2 are not widely available ##REF##10966476##[22]##–##REF##11906816##[24]##. We have assumed constitutive activity and expect as more information becomes available a more accurate model of the regulation of RGS2 and other RGS isoforms will be possible. The SERCA pump helps to bring the cytosolic Ca<sup>2+</sup> concentration back to the resting level after stimulation. We have modeled the SERCA pump as in the Keizer and DeYoung model ##REF##1324019##[25]##. The IP3 dependent opening of ER calcium channels was found to be cooperative ##REF##2452482##[26]## and we have used the Meyer and Stryer model for the IP3-gated channel with a Hill coefficient of four ##REF##1324019##[25]##,##REF##2455890##[27]##. Finally, many other proteins such as calmodulin and the fluorescent indicator Fura-2 bind Ca<sup>2+</sup>. Because our measurements reflect these effects, we have included a general buffer for cytosolic calcium.</p> <title>Isoform Specificity</title> <p>Complement factor 5a activates the C5a receptor which is a Gαi-coupled receptor ##REF##8662841##[28]##. The released Gβγ dimer activates PLCβ2 and PLCβ3 which are lumped and called PLCβ3 in our model because: (i) the activity of Gβγ-activated PLCβ3 has been shown to be greater than Gβγ-activated PLCβ2 in <italic>in-vitro</italic> studies and (ii) Gαq activates both PLCβ2 and PLCβ3 so the structural connections from Gβγ and Gαq to PLCβ2 and PLCβ3 in the model are identical ##REF##8389480##[4]##,##REF##8383116##[29]##. PLCβ1 is activated by Gβγ and Gαq, but RAW264.7 macrophage cells do not express this isoform, so we have not included it in the model. PLCβ3 then catalyzes the hydrolysis of phosphatidylinositol (4,5)-bisphosphate (PIP2) into inositol 1,4,5-trisphosphate (IP3) and diacylglycerol (DAG).</p> <p>UDP stimulates the P2Y6 receptor and the associated Gαq-GTP activates both PLCβ3 ##REF##8387502##[30]## and PLCβ4 ##REF##7929227##[31]##. The GTPase rate of Gαq is increased 1000-fold when bound to PLCβ ##REF##10449728##[5]##. Due to this rapid hydrolysis rate, we have assumed, in our model, that PLCβ3 or PLCβ4 bound Gαq-GTP may only hydrolyze one molecule of PIP2 before releasing Gαq-GDP. Additionally, the Gβγ released by the P2Y6 receptor also activates PLCβ3 ##REF##8387502##[30]##, but does not activate PLCβ4 ##REF##8125982##[32]##.</p> <p>Our model assumes that PLCβ3 does not simultaneously bind Gβγ and Gαq. Indeed, a biochemical study of PLCβ2 activity in reconstituted membrane fractions strongly argues that Gαq and Gβγ do not simultaneously bind this effector ##REF##9931014##[33]##. While this was specifically demonstrated for PLCβ2, we implicitly assume the same holds for PLCβ3 because we lump the two in our model. This is a mechanistic assumption of our model and an interesting issue for future testing with directed experiments.</p> <title>Calcium-Dependent Feedback</title> <p>Though important for response specificity, the dynamical control of calcium release is not limited to the forward pathway in this system. Calcium participates in feedback processes that both enhance and inhibit its own release at multiple points in the pathway. There are four main nodes of calcium-dependent feedback control in our model: PLCβ, IP3 receptor, protein kinase C (PKC) and G protein receptor kinase (GRK).</p> <p>Calcium enhances its own release by binding to the EF-hand domain on PLCβ and is required for PLCβ to hydrolyze PIP2 into IP3 and DAG ##REF##11395409##[34]##. Because the dissociation constant for PLCβ-Ca<sup>2+</sup> in our model is larger than the basal concentration of cytosolic calcium, as more Ca<sup>2+</sup> is released from the ER, more PLCβ-Ca<sup>2+</sup> becomes available to bind Gαq or Gβγ. This positive feedback mechanism accelerates the release of Ca<sup>2+</sup>.</p> <p>In our model, Ca<sup>2+</sup> and IP3 cooperatively open the channel between the ER and the cytosol. It is believed that Ca<sup>2+</sup> initially stimulates the IP3 receptor with maximal stimulatory effect at 100–300 nM ##REF##15189149##[6]##. At higher concentrations, Ca<sup>2+</sup> has an inhibitory effect. We use the IP3 receptor model structure in the Keizer and DeYoung model for this component ##REF##1324019##[25]##.</p> <p>Protein kinase C (PKC) has been shown to phosphorylate PLCβ3 which inhibits PLCβ3 activation due to Gαq and Gβγ ##REF##10893237##[35]##,##REF##12587975##[36]##. PKC is activated when bound to DAG and Ca<sup>2+</sup>
##REF##12954613##[7]##,##REF##15282562##[37]##. Because the preferred order of binding is not entirely known, PKC, DAG and Ca<sup>2+</sup> form a thermodynamic cycle of reversible reaction with only the PKC-DAG-Ca<sup>2+</sup> form active. In our model, the dissociation constant of PKC and Ca<sup>2+</sup> is much greater than the basal Ca<sup>2+</sup> concentration, and upon binding DAG, the PKC-DAG complex has a higher affinity for Ca<sup>2+</sup> making the order of binding preferentially PKC to DAG then PKC-DAG to Ca<sup>2+</sup>. It is not known whether PLCβ4 is also regulated by PKC. We have assumed, in our model, the same mechanism of PKC regulation of PLCβ3 and PLCβ4.</p> <p>The final key calcium-dependent feedback loop in our model is mediated by G protein receptor kinase (GRK). GRK2 phosphorylates and inactivates ligand-bound C5a receptors when activated by PKC and Gβγ. In sequence, PKC phosphorylates GRK2 which causes translocation to the plasma membrane ##REF##14499340##[8]##. When properly localized, GRK2 may bind Gβγ and then phosphorylate the C5a-C5a receptor complex to inactivate it ##REF##10508278##[38]##. This simplified representation of the receptor desensitization mechanism does not include arrestin activity, multiple receptor phosphorylation sites and other fine grain or slower biochemical interactions that may be present <italic>in-vivo</italic>.</p> <title>
<xref ref-type="sec" rid="s4">Single Ligand Experiments</xref>
</title> <p>Having specified the structure of our model, we direct our attention to the parameters. We estimate 20 of the 84 parameters in our model using a dataset composed of 96 Fura-2 time series measurements as described in the Materials and <xref ref-type="sec" rid="s4">Methods</xref> section. Each experiment consists of 3–4 samples from different wells in a 96 well plate. There are 15 experiments spanning 9 doses of C5a and 14 experiments spanning 11 doses of UDP on wild-type cells in the dataset (see ##SUPPL##3##Figure S3##). The dataset also contains calcium measurements on 5 different shRNAi knockdown cell lines constructed by lentiviral infection (see ##SUPPL##4##Figure S4##). The time interval between samples is approximately 3–4 seconds and each time series is approximately 100–300 seconds of post-stimulation data. ##TAB##0##Table 1## shows a summary of the knockdown data used for statistical parameter estimation for this model in addition to the wild-type experiments.</p> <p>We find that our model is generally quantitatively consistent with the experimental data within measurement uncertainty. Where the model is less consistent with the data – specifically for the GRK knockdown experiment – we find the deviation has a reasonable biological explanation. The summary of the dataset and the fit of the model to each single ligand experiment are available in the Supporting Information. We briefly discuss some issues relating to goodness of fit and the Bayesian parameter estimation here.</p> <p>While most optimization procedures produce a point estimate of the parameters that maximize the goodness of fit of the model to the observed data, the Bayesian procedure we have employed here estimates the entire posterior distribution of the parameters given the data. This information is valuable for qualitatively and quantitatively evaluating the precision of the parameters estimates. ##FIG##1##Figure 2## shows, as a qualitative evaluation, that while the a-priori forward and reverse binding rates for the receptors (C5aR and P2YR) are uncorrelated they are correlated in the posterior distribution. The calcium measurements have informed and constrained the posterior estimates of the dissociation constants to be approximately 5 nM and 250 nM for the C5aR and P2YR respectively. We have quantitatively computed marginal highest posterior density (HPD) confidence intervals for each of the twenty parameters we have estimated from the data. Those estimates are shown in ##SUPPL##13##Table S3##. Those parameters with large HPD intervals are not well informed by the measurements and are candidates for directed biochemical experiments.</p> <title>Wild-Type Experiments</title> <p>The calcium response to C5a adapts and returns to the basal level, but the UDP response has a sustained elevated calcium level that slowly decays. ##FIG##2##Figure 3## shows two representative experiments of the response of the wild-type cell to stimulation with C5a and UDP. We expect that the fit to this data will be good because 20 key model parameters were fit using an experimental dataset that included these experiments – the fit is indeed accurate. The point estimate curve is constructed from the maximum a-posteriori parameters from an MCMC chain. The prediction intervals are estimated by Monte Carlo sampling from the posterior parameter distribution and the measurement error distribution conditional on the parameters. The prediction confidence intervals generally cover the observed data.</p> <title>Knockdown Experiments</title> <p>Lentiviral infection is used to introduce small hairpin RNAs to interfere with the translation of the key signaling proteins GRK2, Gαi2, Gαq, PLCβ3 and PLCβ4 ##REF##16945906##[39]##. There are three main sources of uncertainty in the knockdown experiment model predictions: parametric uncertainty, measurement uncertainty and knockdown efficiency uncertainty. We have dealt with the first two sources in the previous section on wild-type experiments. Here we address prediction variability due to knockdown efficiency uncertainty by using nominal parameter values.</p> <p>
##FIG##3##Figure 4## shows simulations and experimental data for three representative knockdown experiments. The upper-left panel of ##FIG##3##Figure 4## shows a GRK knockdown line stimulated with 250 nM C5a. Because GRK2 desensitizes the C5a receptor, we expect that by eliminating the feedback mechanism, the calcium peak will be higher and more sustained. The experimental data as well as the model indeed show that effect. Quantitatively, the model prediction shows a greater effect than the experimental data. A likely reason is that the model only considers one isoform of GRK while there are four isoforms expressed in the RAW264.7 cell line (GRK1,2,4,6). If more than one isoform can desensitize the C5a receptor, the effective knockdown in desensitization function will be less than as measured by western blot analysis on GRK2.</p> <p>While GRK does not desensitize the P2Y receptor in our model, it is a buffer for Gβγ released from Gαq. Reducing the amount of GRK will shift the equilibrium towards more Gβγ bound to PLCβ3 and thus more calcium release even though GRK does not directly feed back on the P2Y6 receptor. The top-right panel in ##FIG##3##Figure 4## shows that, based on the model, the peak intracellular calcium concentration is expected to be very slightly higher in the GRK2 knockdown line when stimulated by 25 µM UDP. A comparison of the experimental peak heights of the wild-type and GRK knockdown cell line data by t-test cannot reject the null hypothesis that the peak heights are equal (<italic>p</italic> = 0.9963). The effect of the GRK knockdown is expected to be so slight that the effect size is overwhelmed by the measurement error in the data. The effect of the uncertainty in the GRK2 knockdown fraction impacts the range of the confidence intervals of the predicted C5a response much more than the confidence intervals of the predicted UDP response which is consistent with GRK2 being a more significant component of the C5a response.</p> <p>Our model structure has PLCβ3 stimulated by either Gβγ or Gαq. Because the C5a response signals only through PLCβ3 the effect of the knockdown is expected to be more pronounced for the C5a response than for the UDP response. The bottom-left panel of ##FIG##3##Figure 4## confirms that the model prediction is consistent with the representative experiment. The UDP response activates PLCβ3 through Gβγ, but also activates PLCβ3 and PLCβ4 with Gαq. Therefore, we expect that the calcium response should be more robust to perturbations in just one of the PLCβ isoforms. The UDP response in the PLCβ3 knockdown line (bottom right panel of ##FIG##3##Figure 4##) shows that our model predicts the knockdown effect to be small relative to the total magnitude of the response in part due to the redundancy in the use of PLCβ isoforms in the UDP response.</p> <p>Because this dataset was used for parameter estimation, the fit of model to the data may overstate the accuracy of the model. Nonetheless, the good fit does suggest that the model warrants being tested in truly predictive experiments; we describe such experiments in the following section.</p> <title>Double Ligand Experiments</title> <p>We examine our model response to a simultaneous stimulation by C5a and UDP because it has been shown experimentally that macrophage cells respond synergistically to such conditions ##REF##16699502##[40]##. To quantify the amount of synergy or non-additivity that is present in the calcium response, a <italic>synergy ratio</italic> is computed for each ligand dose pair. The numerator of the ratio is the peak offset from baseline of the intracellular calcium concentration. The denominator of the ratio is the sum of the peak offsets when the cell or model is stimulated with only one ligand. A synergy is present when the ratio is greater than one implying the peak height is greater than expected from an additive combination of ligand effects. While this is certainly not the only possible measure of synergy it is widely adopted and has been used in previous studies on calcium synergy ##REF##16699502##[40]##.</p> <p>The left panel of ##FIG##4##Figure 5## shows the results of model simulations at nominal parameters for a grid of doses of C5a and UDP. In the dose response surface, there is a ridge of synergistic calcium release for a moderate dose of UDP. We tested the model prediction with the experiment design measuring the synergy ratio at the points denoted as black open circles in the left panel of ##FIG##4##Figure 5##. A χ<sup>2</sup> goodness-of fit test comparing the model expected synergy ratio to the observed synergy ratio fails to reject the null hypothesis that the data were generated by the model mechanism (<italic>p</italic>-value≈1.0). The root-mean-squared error (RMSE) deviation between the predicted and actual experimental data is 0.492. By way of comparison, the RMSE between the data and the null model of no synergy is 1.044. We therefore conclude that the model predictions are consistent with the experimental observations. It should be noted that measurements of synergy in RAW cells are noisy and the ridge occurs at low doses of UDP. Notwithstanding, the phenomenon has been reported ##REF##16699502##[40]## and has been observed by us in this cell line.</p> <p>The right panel of ##FIG##4##Figure 5## shows the same synergy dose response surface but for a GRK knockdown cell line. The synergy ridge observed in the wild-type cell simulation is changed in the GRK knockdown simulation indicating the C5a receptor desensitization mechanism mediated by GRK is important for the synergistic release of calcium. In the next section we pursue this conclusion in more detail, developing a conceptual explanation of the mechanism of crosstalk and synergy within our model.</p>
<title>Discussion</title> <p>G-protein-coupled receptors form a complex network of interacting proteins that generally exhibits the properties of a system in which each receptor signal is buffered from the others. For a minority of ligand combinations, however, crosstalk between pairs of receptors is apparent. Due to the complexity and importance of the system many hypothetical mechanisms have been proposed to explain the crosstalk ##REF##12790797##[2]##. In particular, simultaneous Gβγ and Gαq binding to PLCβ ##REF##8610126##[20]## and Gβγ exchange between Gαi and Gαq-coupled receptors have been proposed as potential mechanisms ##REF##10485876##[19]##. While our model does not eliminate these potential mechanisms, we do show that the mechanism represented in our model is consistent with a full range of experimental data including a variety of doses of C5a and UDP, C5a and UDP stimulation of five different knockdown cell-lines and double-ligand dose response experiments.</p> <p>To our knowledge, this is the first multireceptor GPCR model and the first to address the complex phenomenon of crosstalk between GPCR receptor pathways that has been statistically estimated and validated with experimental data. This important phenomenon plays a role in processes as diverse as chemotaxis and perhaps drug interactions. In our model, the primary mechanism of synergy is due to the cooperative opening of the IP3 receptor. The robustness of the synergy is due to the feedback of GRK on the C5a receptor and the specificity of the synergy is due to the interaction patterns between specific Gα isoforms and PLCβ isoforms. The simultaneous binding model ##REF##8610126##[20]## accounts for the specificity of synergy, but not the robustness pattern of the synergy.</p> <p>We observe in the model that if the Gαq-PLCβ3-Ca<sup>2+</sup> and Gαq-PLCβ4-Ca<sup>2+</sup> binding reactions are inhibited, the system still exhibits synergy. We conclude from this observation that the crosstalk mechanism is mediated by Gβγ. If the binding reaction of Gβγ to phosphorylated GRK2 is removed, the synergy is eliminated. Furthermore, if the GRK2-mediated phosphorylation of complexed C5a receptors is removed, the double ligand response is additive. We deduce then that the synergy mechanism involves GRK2 phosphorylation of complexed C5a receptors. However, GRK2 phosphorylation does not entirely explain the synergy mechanism.</p> <p>In our model, the calcium released from the IP3 receptor is a function of the number of receptor molecules complexed to IP3 raised to the fourth power ##REF##1329108##[41]##. Therefore, for a small range of IP3 concentration, the amount of Ca<sup>2+</sup> released is more than additive (see ##SUPPL##8##Figure S8##). We conclude from our analysis of the model that the synergy ridge in ##FIG##4##Figure 5## arises because the GRK2 mediated mechanism holds the IP3 concentration in this non-additive region for most concentrations of C5a. The UDP response does not have the GRK2 mediated feedback and thus only shows a synergistic response for a small range of UDP concentration. If the GRK2 desensitization is removed from the model, the synergy ridge is removed and synergy is only present at low doses of C5a and UDP (see ##FIG##4##Figure 5##).</p> <p>The Bayesian method we have used for this model has several advantages for the estimation of model parameters in complex mechanistic system models. We have used an informative prior to exclude negative rate constants from the permitted parameter space. We have also used the prior distribution to center our a priori expectations of the rate constant at values obtained from <italic>in-vitro</italic> and other biochemical experiments. The Bayesian update rule allowed us to estimate parameters with our best current dataset and then update those estimates as new data became available from the calcium assay. In this way, we were able to iteratively refine and recalibrate our model with the most recent data available during data collection period for this project. The posterior distribution provides not only an estimate of the rate constants, but the entire distribution, from which we can calculate highest posterior confidence intervals and posterior correlations between parameters. For example, the posterior correlation between the binding and unbinding rates for the UDP-P2Y receptor complex were highly correlated, but those two constants were uncorrelated with the corresponding rates for the C5a-C5a receptor complex reaction even though we imposed no correlations a priori. Finally, the algorithmic methods for collecting ensembles of samples from the posterior distribution have improved considerably in recent years in terms of speed and robustness</p> <p>We have shown that the signal transduction system as it is represented by our model does not require simultaneous binding of Gαq and Gβγ to PLCβ3 to cause a synergistic Ca<sup>2+</sup> response due to simultaneous stimulation by C5a and UDP. We have shown that our representative model is consistent with this experimental dataset in RAW264.7 macrophage cells, but we have not excluded all other potential mechanisms that may be absent or regulated differently in this cell line compared to other macrophage cell lines. Indeed there are a few examples of statistical discrepancies between the model and experiments in our dataset (##SUPPL##14##Table S4##). These differences are substrate for further experimentation and modeling. The purpose of our model is to provide a quantitative tool to aid in reasoning about such complex interacting systems so that meaningful experiments can be designed to explore and understand the biological mechanism.</p>
<p>Conceived and designed the experiments: PF MIJ APA. Performed the experiments: PF MLR. Analyzed the data: PF TD. Contributed reagents/materials/analysis tools: RAR TIR. Wrote the paper: PF MIJ APA.</p> <p>Macrophage cells that are stimulated by two different ligands that bind to G-protein-coupled receptors (GPCRs) usually respond as if the stimulus effects are additive, but for a minority of ligand combinations the response is synergistic. The G-protein-coupled receptor system integrates signaling cues from the environment to actuate cell morphology, gene expression, ion homeostasis, and other physiological states. We analyze the effects of the two signaling molecules complement factors 5a (C5a) and uridine diphosphate (UDP) on the intracellular second messenger calcium to elucidate the principles that govern the processing of multiple signals by GPCRs. We have developed a formal hypothesis, in the form of a kinetic model, for the mechanism of action of this GPCR signal transduction system using data obtained from RAW264.7 macrophage cells. Bayesian statistical methods are employed to represent uncertainty in both data and model parameters and formally tie the model to experimental data. When the model is also used as a tool in the design of experiments, it predicts a synergistic region in the calcium peak height dose response that results when cells are simultaneously stimulated by C5a and UDP. An analysis of the model reveals a potential mechanism for crosstalk between the Gαi-coupled C5a receptor and the Gαq-coupled UDP receptor signaling systems that results in synergistic calcium release.</p> <title>Author Summary</title> <p>The G protein signal transduction system transmits a wide variety of extracellular signals including light, odors, and hormones, to intracellular effectors in diverse cell types in eukaryotes. G-protein-coupled receptors are involved in many diseases including inflammation, cardiac dysfunction, and diabetes, and are the targets of 40–50% of modern drugs. Despite the physiological and pharmacological importance of this signal transduction system it is not known how the system buffers and integrates information at a biochemical level. The multiple receptors expressed by every cell pass their signals through a common set of downstream effectors distinguished by multiple isoforms with slightly different specificities and activities. The coupling among these pathways causes interactions among the signals sent by the different classes of receptors. We have developed a mechanistic model of the G protein signal transduction system from the receptor to the central intracellular second-messenger calcium. We have used statistical methods to integrate a diverse set of experimental data into our model and quantify confidence in our model predictions. We used this model, trained on single receptor data, to predict the signal processing of two G-protein-coupled-receptor signals. Validation experiments support our hypothesized mechanism for dual receptor signal processing and the predictions of the model.</p>
<title>Supporting Information</title>
<p>The authors would like to thank Mel Simon and the Alliance for Cellular Signaling for the experimental data used in this study and many informative and interesting conversations.</p>
<fig id="pcbi-1000185-g001" position="float"><object-id pub-id-type="doi">10.1371/journal.pcbi.1000185.g001</object-id><label>Figure 1</label><caption><title>The model for crosstalk between the Gαi and Gαq pathways depends on both differential specificity and activity for Gαi, Gαq, and Gβγ interactions with PLCβ3 and PLCβ4 to catalyze PIP2 hydrolysis and calcium-dependent feedback control mediated by GRK and PKC.</title><p>Selected model parameters are informed by calcium measurements taken for various ligand doses on wild-type and cell lines with shRNAi knockdowns on the proteins shown in red.</p></caption></fig> <fig id="pcbi-1000185-g002" position="float"><object-id pub-id-type="doi">10.1371/journal.pcbi.1000185.g002</object-id><label>Figure 2</label><caption><title>This figure shows that the single and pairwise marginal posterior distributions for the ligand binding reactions for the P2YR and C5aR receptors.</title><p>The vertical line in the single marginal posterior distributions shows the point estimate that were selected. The posterior distributions show the dissociation constants for the reactions are tightly constrained by the data, while the values of the forward and reverse rates that make up the ratio are not as well constrained by the data. Additionally, as expected the UDP binding rates are not correlated with the C5a binding rates. Marginal posterior distributions for all parameters and a discussion of the point estimate selection can be found in ##SUPPL##2##Figure S2##.</p></caption></fig> <fig id="pcbi-1000185-g003" position="float"><object-id pub-id-type="doi">10.1371/journal.pcbi.1000185.g003</object-id><label>Figure 3</label><caption><title>Model simulations are compared to experimental data.</title><p>The point estimate is computed using the posterior distribution of the parameter as estimated by Markov chain Monte Carlo given the data from 96 experiments on C5a and UDP at various doses in combination with 5 different shRNAi knockdown cell lines. The 95% posterior predictive intervals are estimated by Monte Carlo simulations including both parameter and measurement uncertainty. The measured mean and approximate 95% confidence intervals of four replicates is shown by a black dot and error bar. (A) C5a at 250 nM was introduced at 20s and the experimentally observed pulse in cytosolic calcium concentration is shown. (B) The qualitative shape of the calcium pulse for 25 µM UDP is different than for 250 nM C5a. The pulse does not completely adapt and return to the prestimulated level. For both ligands, the model prediction confidence intervals overlap the data error bars that indicate the model fit is consistent with the data within the measurement uncertainty.</p></caption></fig> <fig id="pcbi-1000185-g004" position="float"><object-id pub-id-type="doi">10.1371/journal.pcbi.1000185.g004</object-id><label>Figure 4</label><caption><title>The model simulation results for GRK and PLCβ3 knockdown cell lines stimulated with C5a and UDP are shown.</title><p>The experimental mean±1 s.d. of 3–4 replicates within one experimental run is shown in black. The knockdown simulation result with nominal knockdown fraction and parameters is shown in red and the wild-type simulation result is shown in green for comparison. Upper and lower model 99% confidence intervals (shown as blue dashed lines) are simulated using the upper and lower knockdown fraction values from ##TAB##0##Table 1##. As expected the Ca<sup>2+</sup> response to C5a in the GRK knockdown line (A) was increased compared to wild-type. The quantitative deviation between the model and data is possibly due to the availability of multiple redundant GRK isoforms. (B) The expected effect of the GRK knockdown on the UDP response is an increase in the cytosolic calcium levels. Because GRK2 does not directly desensitize the P2Y receptor in this model, the effect is likely due to a reduction of sequestration of Gβγ by GRK. (C) The signal transduction of the C5a response is predominantly through the PLCβ3 isoform. The effect of the PLCβ3 knockdown is much greater for C5a than for UDP (D).</p></caption></fig> <fig id="pcbi-1000185-g005" position="float"><object-id pub-id-type="doi">10.1371/journal.pcbi.1000185.g005</object-id><label>Figure 5</label><caption><title>The model is used as a predictive tool to infer the effect of stimulating the cell simultaneously with UDP and C5a that signal through the Gαq and Gαi pathways, respectively.</title><p>Synergy was measured as the ratio of peak height offset from baseline attained from simultaneous stimulation to the peak height offset calculated by the sum of the responses to each ligand individually. (A) Expected synergy ratio as a function of UDP and C5a dose (truncated at 1.5). The simulations show a ridge of synergy at a moderate UDP dose for most C5a doses. The black circles indicate dose combinations points of experiments that were conducted to test the model. (B) Expected synergy ratio as a function of UDP and C5a dose for a simulated GRK2 knockdown cell line. Without the GRK-mediated negative feedback to keep the IP3 generation from the C5a receptor within the non-linear range of calcium release the ridge in the synergy dose response is diminished. The synergy in the GRK knockdown simulation is not entirely eliminated because the shRNAi knockdown of GRK does not constitute a complete loss-of-function and low concentrations of ligand are still able to synergize. Furthermore, the asymmetric synergy dose response surface is more symmetric in the GRK knockdown simulation because the asymmetric calcium-dependent feedback mechanism is reduced.</p></caption></fig>
<table-wrap id="pcbi-1000185-t001" position="float"><object-id pub-id-type="doi">10.1371/journal.pcbi.1000185.t001</object-id><label>Table 1</label><caption><title>Dataset used for parameter estimation.</title></caption><alternatives><table frame="hsides" rules="groups"><colgroup span="1"><col align="left" span="1"/><col align="center" span="1"/><col align="center" span="1"/><col align="center" span="1"/><col align="center" span="1"/><col align="center" span="1"/><col align="center" span="1"/><col align="center" span="1"/><col align="center" span="1"/><col align="center" span="1"/><col align="center" span="1"/><col align="center" span="1"/></colgroup><thead><tr><td align="left" rowspan="1" colspan="1">Cell Line</td><td colspan="2" align="left" rowspan="1">Measured Fraction Knockdown</td><td colspan="3" align="left" rowspan="1">Model Value</td><td colspan="6" align="left" rowspan="1">Sample Size</td></tr><tr><td align="left" rowspan="1" colspan="1"/><td colspan="2" align="left" rowspan="1"/><td colspan="3" align="left" rowspan="1"/><td colspan="3" align="left" rowspan="1">C5a</td><td colspan="3" align="left" rowspan="1">UDP</td></tr><tr><td align="left" rowspan="1" colspan="1"/><td align="left" rowspan="1" colspan="1">qRT-PCR</td><td align="left" rowspan="1" colspan="1">Western</td><td align="left" rowspan="1" colspan="1">Nominal</td><td align="left" rowspan="1" colspan="1">Lower</td><td align="left" rowspan="1" colspan="1">Upper</td><td align="left" rowspan="1" colspan="1">&lt;10 nM</td><td align="left" rowspan="1" colspan="1">10–100 nM</td><td align="left" rowspan="1" colspan="1">&gt;100 nM</td><td align="left" rowspan="1" colspan="1">&lt;1 µM</td><td align="left" rowspan="1" colspan="1">1–10 µM</td><td align="left" rowspan="1" colspan="1">&gt;10 µM</td></tr></thead><tbody><tr><td align="left" rowspan="1" colspan="1">Wild-type</td><td align="left" rowspan="1" colspan="1">–</td><td align="left" rowspan="1" colspan="1">–</td><td align="left" rowspan="1" colspan="1">–</td><td align="left" rowspan="1" colspan="1">–</td><td align="left" rowspan="1" colspan="1">–</td><td align="left" rowspan="1" colspan="1">4</td><td align="left" rowspan="1" colspan="1">8</td><td align="left" rowspan="1" colspan="1">3</td><td align="left" rowspan="1" colspan="1">5</td><td align="left" rowspan="1" colspan="1">5</td><td align="left" rowspan="1" colspan="1">4</td></tr><tr><td align="left" rowspan="1" colspan="1">GRK2 (2)</td><td align="left" rowspan="1" colspan="1">90%±7%, <italic>n</italic> = 5</td><td align="left" rowspan="1" colspan="1">40%±6%, <italic>n</italic> = 6</td><td align="left" rowspan="1" colspan="1">40.0%</td><td align="left" rowspan="1" colspan="1">22.0%</td><td align="left" rowspan="1" colspan="1">58.0%</td><td align="left" rowspan="1" colspan="1">2</td><td align="left" rowspan="1" colspan="1">12</td><td align="left" rowspan="1" colspan="1">2</td><td align="left" rowspan="1" colspan="1">3</td><td align="left" rowspan="1" colspan="1">1</td><td align="left" rowspan="1" colspan="1">5</td></tr><tr><td align="left" rowspan="1" colspan="1">Gai2 (3)</td><td align="left" rowspan="1" colspan="1">83%±5%, <italic>n</italic> = 4</td><td align="left" rowspan="1" colspan="1">73%±6%, <italic>n</italic> = 5</td><td align="left" rowspan="1" colspan="1">73.0%</td><td align="left" rowspan="1" colspan="1">55.0%</td><td align="left" rowspan="1" colspan="1">91.0%</td><td align="left" rowspan="1" colspan="1">–</td><td align="left" rowspan="1" colspan="1">5</td><td align="left" rowspan="1" colspan="1">–</td><td align="left" rowspan="1" colspan="1">5</td><td align="left" rowspan="1" colspan="1">–</td><td align="left" rowspan="1" colspan="1">7</td></tr><tr><td align="left" rowspan="1" colspan="1">Gaq (3)</td><td align="left" rowspan="1" colspan="1">70%±8%, <italic>n</italic> = 7</td><td align="left" rowspan="1" colspan="1">66%±23%, <italic>n</italic> = 2</td><td align="left" rowspan="1" colspan="1">66.0%</td><td align="left" rowspan="1" colspan="1">0.0%</td><td align="left" rowspan="1" colspan="1">95.0%</td><td align="left" rowspan="1" colspan="1">–</td><td align="left" rowspan="1" colspan="1">3</td><td align="left" rowspan="1" colspan="1">–</td><td align="left" rowspan="1" colspan="1">1</td><td align="left" rowspan="1" colspan="1">–</td><td align="left" rowspan="1" colspan="1">3</td></tr><tr><td align="left" rowspan="1" colspan="1">PLCβ3 (1)</td><td align="left" rowspan="1" colspan="1">–</td><td align="left" rowspan="1" colspan="1">83%±15%, <italic>n</italic> = 3</td><td align="left" rowspan="1" colspan="1">83.0%</td><td align="left" rowspan="1" colspan="1">38.0%</td><td align="left" rowspan="1" colspan="1">100.0%</td><td align="left" rowspan="1" colspan="1">–</td><td align="left" rowspan="1" colspan="1">3</td><td align="left" rowspan="1" colspan="1">–</td><td align="left" rowspan="1" colspan="1">–</td><td align="left" rowspan="1" colspan="1">–</td><td align="left" rowspan="1" colspan="1">3</td></tr><tr><td align="left" rowspan="1" colspan="1">PLCβ4 (1)</td><td align="left" rowspan="1" colspan="1">87%±6%, <italic>n</italic> = 5</td><td align="left" rowspan="1" colspan="1">–</td><td align="left" rowspan="1" colspan="1">87.0%</td><td align="left" rowspan="1" colspan="1">69.0%</td><td align="left" rowspan="1" colspan="1">100.0%</td><td align="left" rowspan="1" colspan="1">–</td><td align="left" rowspan="1" colspan="1">4</td><td align="left" rowspan="1" colspan="1">–</td><td align="left" rowspan="1" colspan="1">4</td><td align="left" rowspan="1" colspan="1">–</td><td align="left" rowspan="1" colspan="1">4</td></tr></tbody></table></alternatives></table-wrap>
<disp-formula></disp-formula> <disp-formula></disp-formula> <inline-formula></inline-formula> <inline-formula></inline-formula>
<supplementary-material content-type="local-data" id="pcbi.1000185.s001"><label>Dataset S1</label><caption><p>(0.11 MB DOC)</p></caption></supplementary-material> <supplementary-material content-type="local-data" id="pcbi.1000185.s002"><label>Figure S1</label><caption><p>This figure shows exemplar MCMC realizations for parameter k109f (the UDP+P2YR forward binding rate) from three independent chains. The chains have converged to the stationary distribution which is the posterior distribution as measured by the PSRF (see <xref ref-type="sec" rid="s4">Materials and Methods</xref>).</p><p>(0.20 MB DOC)</p></caption></supplementary-material> <supplementary-material content-type="local-data" id="pcbi.1000185.s003"><label>Figure S2</label><caption><p>Posterior distributions and correlations The first figure shows that the pairwise marginal posterior distributions for the ligand binding reactions for P2YR and C5aR. The posterior distributions show the dissociation constants for the reactions are tightly constrained by the data, while the values of the forward and reverse rates that make up the ratio are not as well constrained by the data. Additionally, the UDP binding rates are not correlated with the C5a binding rates. k108f and k108r are the P2YR forward and reverse rates and k101f and k101r are the C5aR rates. The next two figures show the one-way marginal posterior density estimates from three independent MCMC chains with approximately 30,000 samples. The 20 estimates parameters are along the rows and the independent chains are along the columns. In each plot, the light blue density is the prior density and the green, purple and orange densities are the posterior densities. The vertical line shows the parameter value used in the model simulations in the paper and listed in ##SUPPL##13##Table S3##. All of the densities are plotted on a log scale. Each marginal posterior distribution estimate is constructed from independent MCMC chains. The results from each chain (three of them) are shown in the columns of the second figure below. In some cases the algorithm sampled heavily from one mode that was not explored as heavily by another chain. However, the PSRF criterion used to assay convergence and a visual inspection of overall posterior density correspondence do indicate that the posterior distributions are sufficiently sampled by all three chains in aggregate. Furthermore, the fit of the model to the data as shown in ##SUPPL##3##Figure S3## shows that the model point estimates are effective in fitting the actual calcium measurements.</p><p>(0.55 MB DOC)</p></caption></supplementary-material> <supplementary-material content-type="local-data" id="pcbi.1000185.s004"><label>Figure S3</label><caption><p>Peak height dose response. This figure shows the single ligand calcium dose responses for C5a and UDP stimulation.</p><p>(0.21 MB DOC)</p></caption></supplementary-material> <supplementary-material content-type="local-data" id="pcbi.1000185.s005"><label>Figure S4</label><caption><p>Knockdown simulations. This figure shows representative simulations and data for each knockdown experiment. A complete set of all 96 experiments is provided in a supplementary folder.</p><p>(0.69 MB DOC)</p></caption></supplementary-material> <supplementary-material content-type="local-data" id="pcbi.1000185.s006"><label>Figure S5</label><caption><p>Input model fit. This figure shows the input model (described in <xref ref-type="sec" rid="s4">Materials and Methods</xref>) fit to the FITC measurements. The ligand concentration that the cell sees does not transit instantaneously from 0 to the final concentration. The ligand concentration is expected to take an amount of time that is significant on the scale of the measurements made for this study to reach the final concentration.</p><p>(0.12 MB DOC)</p></caption></supplementary-material> <supplementary-material content-type="local-data" id="pcbi.1000185.s007"><label>Figure S6</label><caption><p>Large pathway diagram.</p><p>(0.18 MB DOC)</p></caption></supplementary-material> <supplementary-material content-type="local-data" id="pcbi.1000185.s008"><label>Figure S7</label><caption><p>System of differential equations. This figure shows the complete set of differential equations used to simulate the model. These equations are also available in the source c code for the model supplied. This system of equations with the initial conditions and nominal parameter values reported in ##SUPPL##11##Table S1## and ##SUPPL##12##Table S2##, respectively, completely define the model and allow for the reproduction of the simulations used in this paper on any platform.</p><p>(1.62 MB DOC)</p></caption></supplementary-material> <supplementary-material content-type="local-data" id="pcbi.1000185.s009"><label>Figure S8</label><caption><p>Hill function self-synergy. Consider a Hill function, .  is a dimensionless critical concentration <italic>y</italic>*, below which self-synergy will occur. Based on the analysis, we conclude that: (i) <italic>n</italic> must be greater than 1 for self-synergy to occur, (ii) self synergy never occurs if the concentration <italic>x</italic> exceeds equilibrium constant <italic>K</italic> (<italic>y</italic>&gt;1), and (iii) for <italic>n</italic>&gt;2, there is a large range of concentration for self-synergy. In the G protein model, <italic>x</italic>, is the concentration of IP3-IP3R, <italic>H</italic>(<italic>x</italic>) is the rate of change in cytosolic calcium concentration and <italic>n</italic> = 4. We have tested the validity of this self synergy hypothesis by stimulating the cells with both 20 nM UDP and 40 nM UDP (data not shown). Though at such low ligand concentrations, the measurement variability is high, we observed that the synergy ratio, on average was 1.17 compared to a value of 1.25 predicted by the model.</p><p>(0.02 MB DOC)</p></caption></supplementary-material> <supplementary-material content-type="local-data" id="pcbi.1000185.s010"><label>Figure S9</label><caption><p>Parameter sensitivity analysis. The parameter of interest is varied by 10% while all other parameters are kept constant. The parameters are grouped according to their functionalities. The sensitivity coefficient is the ratio of the relative change in the peak height to the relative change in the parameter value. The four most sensitive parameters (sensitivity coefficient &gt;2) in the Cacyt category are Vqssk50 (IP3+IP3K_a-&gt;IP4+IP3K_a (Vmax)), Kqssk50 (IP3+IP3K_a-&gt;IP4+IP3K_a (Km)), a1 (Ca leak into the cell from outside), and Kex (Na/Ca exchange activation const). The top 3 most sensitive parameters in the PLCb3 category are: k21bf* (PLCb3_Ca_Gbg_PIP2-&gt;PLCb3_Ca_Gbg+IP3+DAG), k20f (Gbg+PLCb3_Ca-&gt;PLCb3_Ca_Gbg), k21af* (PLCb3_Ca_Gbg+PIP2-&gt;PLCb3_Ca_Gbg_PIP2). A star next to the parameter name indicates it was estimated.</p><p>(0.04 MB DOC)</p></caption></supplementary-material> <supplementary-material content-type="local-data" id="pcbi.1000185.s011"><label>Protocol S1</label><caption><p>FITC protocol.</p><p>(0.03 MB DOC)</p></caption></supplementary-material> <supplementary-material content-type="local-data" id="pcbi.1000185.s012"><label>Table S1</label><caption><p>Model initial conditions. This table shows the initial conditions used for the model. The model was run for sufficient time for the species states in the model to reach equilibrium before ligand stimulation was added. The number of molecules was calculated using a cell volume of 1 pL.</p><p>(0.05 MB DOC)</p></caption></supplementary-material> <supplementary-material content-type="local-data" id="pcbi.1000185.s013"><label>Table S2</label><caption><p>Model parameters. This table shows the nominal parameters used for the model. Parameter distributions that were estimated are shown as shaded rows and with a star next to the parameter name in the table. The prior distribution for each parameter is as described in the <xref ref-type="sec" rid="s4">Materials and Methods</xref> section with mean value specified by the column labeled “prior”.</p><p>(0.14 MB DOC)</p></caption></supplementary-material> <supplementary-material content-type="local-data" id="pcbi.1000185.s014"><label>Table S3</label><caption><p>Parameter posterior uncertainty and references. This table shows the HPD intervals as computed by the R CODA library function “hpdinterval”. HPD intervals for each of the three MCMC chains were calculated and the union of those intervals is reported for each parameter in this table. The prior value reported in ##SUPPL##12##Table S2## was set using information from references listed in the appropriate column. The references used to form the basis of the parameter estimates are shown in the last column.</p><p>(0.37 MB DOC)</p></caption></supplementary-material> <supplementary-material content-type="local-data" id="pcbi.1000185.s015"><label>Table S4</label><caption><p>Goodness of fit evaluation. We use the mean squared error criterion to evaluate the goodness of our model fit to the data. We have used this data in the estimation procedure and thus does not constitute a true validation. However, we show that in general our model fits the bulk of the data. Those areas of lack-of-fit are usually due to extraordinary experiment-to-experiment variation and in some cases point to unaccounted mechanisms. We elaborate on one such mechanism (multiple GRK isoforms) in the text of the article.</p><p>(0.27 MB DOC)</p></caption></supplementary-material>
<table-wrap-foot><fn id="nt101"><label/><p>Five different cell lines that have a perturbation in the level of a key signal transduction protein were constructed by shRNAi lentiviral infection. The calcium response from these cell lines in addition to the wild-type cell line were used to fit relevant parameters in the model. Because shRNAi does not entirely remove the protein product, the fraction knockdown was estimated by qRT-PCR and by Western blot analysis. The standard error (se) was computed for each estimate and the upper and lower confidence intervals were computed as ±3·se. The knockdown confidence intervals are used in the GPCR model to construct prediction confidence intervals for the calcium response. Where several cell lines were constructed for each knockdown, the best was selected and reported in parenthesis. The sample size for each knockdown-ligand dose combination is shown in the last 6 columns.</p></fn></table-wrap-foot> <fn-group><fn fn-type="COI-statement"><p>The authors have declared that no competing interests exist.</p></fn><fn fn-type="financial-disclosure"><p>MLR was supported by a Department of Energy Computational Science Graduate Fellowship (DE-FG02-97ER25308). RAR and TIR were supported by National Institutes of Health (NIH) Grant GM62114. PF, TD, and APA were supported by the Alliance for Cellular Signaling (NIH Grant U54 GM62114). APA would also like to acknowledge support of the Howard Hughes Medical Institute. MIJ was supported by a grant from Microsoft Research and a grant from Yahoo! Research.</p></fn></fn-group>
<graphic xlink:href="pcbi.1000185.g001"/> <graphic id="pcbi-1000185-t001-1" xlink:href="pcbi.1000185.t001"/> <graphic xlink:href="pcbi.1000185.g002"/> <graphic xlink:href="pcbi.1000185.g003"/> <graphic xlink:href="pcbi.1000185.g004"/> <graphic xlink:href="pcbi.1000185.g005"/> <graphic xlink:href="pcbi.1000185.e001.jpg" mimetype="image" position="float"/> <graphic xlink:href="pcbi.1000185.e002.jpg" mimetype="image" position="float"/> <inline-graphic xlink:href="pcbi.1000185.e003.jpg" mimetype="image"/> <inline-graphic xlink:href="pcbi.1000185.e004.jpg" mimetype="image"/>
<media xlink:href="pcbi.1000185.s001.doc"><caption><p>Click here for additional data file.</p></caption></media> <media xlink:href="pcbi.1000185.s002.doc"><caption><p>Click here for additional data file.</p></caption></media> <media xlink:href="pcbi.1000185.s003.doc"><caption><p>Click here for additional data file.</p></caption></media> <media xlink:href="pcbi.1000185.s004.doc"><caption><p>Click here for additional data file.</p></caption></media> <media xlink:href="pcbi.1000185.s005.doc"><caption><p>Click here for additional data file.</p></caption></media> <media xlink:href="pcbi.1000185.s006.doc"><caption><p>Click here for additional data file.</p></caption></media> <media xlink:href="pcbi.1000185.s007.doc"><caption><p>Click here for additional data file.</p></caption></media> <media xlink:href="pcbi.1000185.s008.doc"><caption><p>Click here for additional data file.</p></caption></media> <media xlink:href="pcbi.1000185.s009.doc"><caption><p>Click here for additional data file.</p></caption></media> <media xlink:href="pcbi.1000185.s010.doc"><caption><p>Click here for additional data file.</p></caption></media> <media xlink:href="pcbi.1000185.s011.doc"><caption><p>Click here for additional data file.</p></caption></media> <media xlink:href="pcbi.1000185.s012.doc"><caption><p>Click here for additional data file.</p></caption></media> <media xlink:href="pcbi.1000185.s013.doc"><caption><p>Click here for additional data file.</p></caption></media> <media xlink:href="pcbi.1000185.s014.doc"><caption><p>Click here for additional data file.</p></caption></media> <media xlink:href="pcbi.1000185.s015.doc"><caption><p>Click here for additional data file.</p></caption></media>
[{"label": ["10"], "element-citation": ["\n"], "surname": ["Berg", "Tymoczko", "Stryer"], "given-names": ["JM", "JL", "L"], "year": ["2002"], "source": ["Biochemistry. 5th ed"], "publisher-loc": ["Bethesda (Maryland)"], "publisher-name": ["W.H. Freeman; National Center for Biotechnology Information"]}, {"label": ["42"], "element-citation": ["\n"], "surname": ["Hindmarsh", "Brown", "Grant", "Lee", "Serban"], "given-names": ["AC", "PN", "KE", "SL", "R"], "year": ["2005"], "article-title": ["SUNDIALS: Suite of Nonlinear and Differential/Algebraic Equation Solvers."], "source": ["ACM Trans Math Softw"], "volume": ["31"], "fpage": ["363"], "lpage": ["396"]}, {"label": ["45"], "element-citation": ["\n"], "surname": ["Robert", "Casella"], "given-names": ["CP", "G"], "year": ["2004"], "source": ["Monte Carlo Statistical Methods"], "publisher-loc": ["New York"], "publisher-name": ["Springer"]}, {"label": ["46"], "element-citation": ["\n"], "surname": ["Gelman", "Rubin"], "given-names": ["A", "DB"], "year": ["1992"], "article-title": ["Inference from iterative simulation using multiple sequences."], "source": ["Stat Sci"], "volume": ["7"], "fpage": ["457"], "lpage": ["472"]}]
PLoS Comput Biol. 2008 Sep 26; 4(9):e1000185