A Dual Receptor Cross-talk Model of G protein-coupled Signal
Transduction

Running Head: Dual Receptor GPCR Crosstalk Model

Patrick Flaherty^1^, Mala L. Radhakrishnan^2^, Tuan Dinh^3^, Michael I.
Jordan^4^, Adam P. Arkin^5\*^

^1^Department of Electrical Engineering and Computer Sciences, UC
Berkeley, Berkeley, CA, 94720

^2^Department of Chemistry, Computer Science and Artificial Intelligence
Laboratory, Massachusetts Institute of Technology, Cambridge, MA 02139

^3^Department of Bioengineering, UC Berkeley, Berkeley, CA, 94720

^4^Department of Statistics and Computer Science Division, UC Berkeley,
Berkeley, CA 94720

^5^ Physical Biosciences Division, Lawrence Berkeley Laboratory, Howard
Hughes Medical Institute, Department of Bioengineering, UC Berkeley,
Berkeley, CA 94720

^\*^Corresponding Author: Adam P. Arkin, 717 Potter Street, Berkeley,
CA, 94710,

V: (510) 495-2366, F: (510) 486-6219, aparkin@lbl.gov

Number of text pages: **16**

Number of figures: **5**

Number of tables: **1**

Number of words in the abstract: **216**

Number of characters in the paper: **44,725**

**Abstract**

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
factor 5a (C5a) and uridine diphosphate (UDP) on the intracellular
second messenger calcium to elucidate the principles that govern the
processing of multiple signals by G protein-coupled receptors. 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.

# Author Summary

# 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 receptor. 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 Bayesian statistical methods to integrate a diverse set of experimental data into our model and quantify confidence intervals on our model predictions. We used this model, trained on single receptor data, to predict the non-additive signal processing of two GPCR signals. Validation experiments support our hypothesized mechanism for dual receptor signal processing and the predictions of the model. 

# Introduction

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 (Kroeze et al.).
Though the system is physiologically and pharmacologically important,
the mechanism by which the system integrates multiple signals is not
well understood (Werry et al.).

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) (Casey and
Gilman; Wu et al.). In addition to its catalytic activity, PLCβ acts as
a GTPase for Gα-GTP (Mukhopadhyay and Ross 1999). IP3 binds to specific
receptor-channels on the membrane of the ER to release Ca^2+^ into the
cytosol (Patterson et al. 2004). DAG and Ca^2+^ bind to and activate
protein kinase C (PKC) which may phosphorylate and inactivate specific
PLCβ isoforms (Ananthanarayanan et al.). G protein receptor kinase (GRK)
is activated once it is phosphorylated by PKC (Penela et al.) and is
localized to the plasma membrane by Gβγ (Pitcher et al.). 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 (Penela et
al.). Activated GRK can then phosphorylate specific GPCRs which leads to
receptor inactivation -- perhaps directly or by arrestin activity
(Penela et al.). 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
(Berg et al. 2002).

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 (Allegretti et al.).
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 (Warny et al.). 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(Yoshioka et al. 2001).

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 (Lukas 2004). Mishra & Bhalla built a model to
investigate the role of IP4 as a signal coincidence detector in the GPCR
pathway (Mishra and Bhalla). The model by Lemon et al. predicts the
calcium response to UTP stimulation and is the closest in focus to our
model (Lemon et al.). 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 (Maurya and Subramaniam 2007a, 2007b).

Several hypotheses for the mechanism of crosstalk and synergy among
GPCR-mediated pathways have been proposed. Cross-talk 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 (Quitterer
and Lohse 1999). Zhu et al. speculated that PLC is under either
conditional or dual regulation of Gβγ and Gα (Zhu and Birnbaumer 1996).
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 (Werry et al. 2003).

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.

**Results**

There are two main features of the structure of our model, shown in
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.

**\[FIGURE 1 HERE\]**

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^2+^/Ca^2+^ exchanger, a
sarco(endo)plasmic reticulum Ca^2+^-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 Figure S7 (supplementary
information). The parameters and initial conditions are in Table S2 and
Table S1 respectively (supplementary information). Where available, we
have relied on *in-vitro* or *in-vivo* biochemical experiments for the
reactions and parameter values (see supplementary 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 Methods 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.

We briefly discuss the reactions involving the Na^2+^/Ca^2+^ 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 (Kehrl 1998); 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 (Ross and Wilkie 2000; Cunningham et al. 2001; Kehrl
and Sinnarajah 2002). 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^2+^ concentration back to the
resting level after stimulation. We have modeled the SERCA pump as in
the Keizer and DeYoung model (Keizer and De Young 1992). The IP3
dependent opening of ER calcium channels was found to be cooperative
(Meyer et al. 1988) and we have used the Meyer and Stryer model for the
IP3-gated channel with a Hill coefficient of four (Meyer and Stryer
1988; Keizer and De Young 1992). Finally, many other proteins such as
calmodulin and the fluorescent indicator Fura-2 bind Ca^2+^. Because our
measurements reflect these effects, we have included a general buffer
for cytosolic calcium.

### Isoform Specificity

Complement factor 5a activates the C5a receptor which is a Gαi-coupled
receptor (Jiang et al. 1996). 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 *in-vitro* 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 (Park et al. 1993; Wu et al. 1993).
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).

UDP stimulates the P2Y6 receptor and the associated Gαq-GTP activates
both PLCβ3 (Smrcka and Sternweis 1993) and PLCβ4 (Lee et al. 1994). The
GTPase rate of Gαq is increased 1000-fold when bound to
PLCβ (Mukhopadhyay and Ross 1999). 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 (Smrcka and
Sternweis 1993), but does not activate PLCβ4 (Jiang et al. 1994).

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 (Runnels and Scarlata 1999). 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.

### Calcium-dependent Feedback

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

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 (Rhee
2001). Because the dissociation constant for PLCβ-Ca^2+^ in our model is
larger than the basal concentration of cytosolic calcium, as more Ca^2+^
is released from the ER, more PLCβ-Ca^2+^ becomes available to bind Gαq
or Gβγ. This positive feedback mechanism accelerates the release of
Ca^2+^.

In our model, Ca^2+^ and IP3 cooperatively open the channel between the
ER and the cytosol. It is believed that Ca^2+^ initially stimulates the
IP3 receptor with maximal stimulatory effect at 100-300nM (Patterson et
al. 2004). At higher concentrations, Ca^2+^ has an inhibitory effect. We
use the IP3 receptor model structure in the Keizer and DeYoung model for
this component (Keizer and De Young 1992).

Protein kinase C (PKC) has been shown to phosphorylate PLCβ3 which
inhibits PLCβ3 activation due to Gαq and Gβγ (Yue et al.; Litosch 2002).
PKC is activated when bound to DAG and Ca^2+^ (Ananthanarayanan et al.
2003; Spitaler and Cantrell 2004). Because the preferred order of
binding is not entirely known, PKC, DAG and Ca^2+^ form a thermodynamic
cycle of reversible reaction with only the PKC-DAG-Ca^2+^ form active.
In our model, the dissociation constant of PKC and Ca^2+^ is much
greater than the basal Ca^2+^ concentration, and upon binding DAG, the
PKC-DAG complex has a higher affinity for Ca^2+^ making the order of
binding preferentially PKC to DAG then PKC-DAG to Ca^2+^. 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.

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 (Penela et al. 2003). When properly localized, GRK2 may bind
Gβγ and then phosphorylate the C5a-C5a receptor complex to inactivate it
(Langkabel et al. 1999). 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 *in-vivo*.

# *Single Ligand Experiments*

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
data set composed of 96 Fura-2 time series measurements as described in
the Materials and Methods 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 data set. The data set also contains
calcium measurements on 5 different shRNAi knockdown cell lines
constructed by lentiviral infection. The time interval between samples
is approximately 3-4 seconds and each time series is approximately
100-300 seconds of post-stimulation data. Table 1 shows a summary of the
knockdown data used for statistical parameter estimation for this model
in addition to the wild-type experiments.

**\[TABLE 1 HERE\]**

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 data set and the fit of the model to
each single ligand experiment are available in the supplementary
information. We briefly discuss some issues relating to goodness of fit
and the Bayesian parameter estimation here.

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. 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 5nM and 250nM 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
Table S3. Those parameters with large HPD intervals are not well
informed by the measurements and are candidates for directed biochemical
experiments.

**\[FIGURE 2 HERE\]**

### Wild-type Experiments

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. 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 data set 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.

**\[FIGURE 3 HERE\]**

### Knockdown Experiments

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 (Shin et al. 2006). 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.

**\[FIGURE 4 HERE\]**

Figure 4 shows simulations and experimental data for three
representative knockdown experiments. The upper-left panel of Figure 4
shows a GRK knockdown line stimulated with 250nM 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.

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 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 (p=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.

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 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 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.

Because this data set 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.

# *Double Ligand Experiments*

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 (Natarajan et al.). To
quantify the amount of synergy or non-additivity that is present in the
calcium response, a *synergy ratio* 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 (Natarajan et al.).

**\[FIGURE 5 HERE\]**

The left panel of 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 Figure 5. A χ^2^ 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 (p-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
(Natarajan et al.) and has been observed by us in this cell line.

The right panel of 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.

**Discussion**

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
(Werry et al. 2003). In particular, simultaneous Gβγ and Gαq binding to
PLCβ (Zhu and Birnbaumer 1996) and Gβγ exchange between Gαi and
Gαq-coupled receptors have been proposed as potential mechanisms
(Quitterer and Lohse 1999). 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.

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 (Zhu and Birnbaumer) accounts
for the specificity of synergy, but not the robustness pattern of the
synergy.

We observe in the model that if the Gαq-PLCβ3-Ca^2+^ and
Gαq-PLCβ4-Ca^2+^ 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.

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 (De Young and Keizer 1992). Therefore, for a small range of
IP3 concentration, the amount of Ca^2+^ released is more than additive
(see Figure S8, supplementary information). We conclude from our
analysis of the model that the synergy ridge in 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 Figure 5).

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 *in-vitro* and other biochemical
experiments. The Bayesian update rule allowed us to estimate parameters
with our best current data set 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

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^2+^ response due to simultaneous
stimulation by C5a and UDP. We have shown that our representative model
is consistent with this experimental data set 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 data set (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.

**Materials and Methods**

The model equations are given in Figure S7 of the supplementary
information. The initial conditions and parameter values are in Tables
S1 and S2 respectively. All the data used in this work and a stand-alone
implementation of the model is provided at
http://genomics.lbl.gov/supplemental/flaherty-gpcr/. The model was
simulated using CVODE (Hindmarsh et al. 2005) and the GNU Scientific
Library. Further details on materials and methods are available in
supplementary information.

## Experimental Methods

## Intracellular free calcium in cultured adherent RAW264.7 cells was measured in a 96-well plate format using the Ca^2+^-sensitive fluorescent dye Fura-2 (Grynkiewicz et al.; Tsien). 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 [***http://www.signaling-gateway.org***](http://www.signaling-gateway.org/)). 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 Protocol #1 in the supplementary information.

## Statistical Inference

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 data set were estimated and
are denoted with a star in 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.

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.

The likelihood is a function of the parameters (θ) and links the prior
distribution with the posterior distribution under Bayes rule

![](media/image1.pct)

where y denotes the observed data.

In our model, the likelihood function is a Gaussian distribution
according to the non-linear regression equation
$y = f(q) + e,e\text{\textasciitilde}N\left( 0,s^{2} \right)$, where
![](media/image2.pct) 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.

The Metropolis-Hastings algorithm (Robert and Casella 2004) was used to
estimate the posterior density of the parameters Pr(θ\|*y*). Three
independent chains were simulated from different initial parameter
values. To assess convergence of the posterior distribution estimate, we
used the Gelman-Rubin potential scale reduction factor (PSRF) (Gelman
and Rubin 1992). 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.

Posterior prediction confidence intervals were constructed using the
percentiles from the predictive distribution approximated with 2000
Monte Carlo samples from ![](media/image3.pct) at each of 100 simple
random samples from ![](media/image4.pct) obtained from

![](media/image5.pct)

where
$\text{Pr}\left( y_{\text{new}}|q_{i} \right)\text{\textasciitilde}N\left( f(q),s^{2} \right)$
and $s^{2}$ 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 $s^{2}$. A small factor of 1nM^2^ was added to each
variance estimate to bound variance estimates away from zero.

# Acknowledgements

PF, APA and MIJ designed the research. PF and MLR performed the
research. PF and TD analyzed the data. PF, MIJ and APA wrote the paper.
The authors would like to thank William Seaman, Robert Rebres, Tamara
Roach, Mel Simon and the Alliance for Cellular Signaling for the
experimental data used in this study and many informative and
interesting conversations.

Table 1:

The data set used for parameter estimation is shown in this table. 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.

Table 1:

<table>
<colgroup>
<col style="width: 10%" />
<col style="width: 9%" />
<col style="width: 1%" />
<col style="width: 11%" />
<col style="width: 8%" />
<col style="width: 6%" />
<col style="width: 7%" />
<col style="width: 7%" />
<col style="width: 7%" />
<col style="width: 8%" />
<col style="width: 7%" />
<col style="width: 6%" />
<col style="width: 8%" />
</colgroup>
<tbody>
<tr class="odd">
<td> </td>
<td><strong> </strong></td>
<td colspan="2"><strong> </strong></td>
<td><strong> </strong></td>
<td><strong> </strong></td>
<td><strong> </strong></td>
<td colspan="6"><strong>Sample Size</strong></td>
</tr>
<tr class="even">
<td></td>
<td colspan="3"><strong>Measured Fraction Knockdown</strong></td>
<td colspan="3"><strong>Model Value</strong></td>
<td colspan="3"><strong>C5a</strong></td>
<td colspan="3"><strong>UDP</strong></td>
</tr>
<tr class="odd">
<td>Cell Line</td>
<td colspan="2">qRT-PCR</td>
<td>Western</td>
<td>Nominal</td>
<td>Lower</td>
<td>Upper</td>
<td>&lt;10nM</td>
<td>10 - 100nM</td>
<td>&gt;100nM</td>
<td>&lt;1μM</td>
<td><p>1 –</p>
<p>10μM</p></td>
<td>&gt; 10μM</td>
</tr>
<tr class="even">
<td>Wild-type</td>
<td colspan="2">-</td>
<td>-</td>
<td>-</td>
<td>-</td>
<td>-</td>
<td>4</td>
<td>8</td>
<td>3</td>
<td>5</td>
<td>5</td>
<td>4</td>
</tr>
<tr class="odd">
<td>GRK2 (2)</td>
<td colspan="2">90% ± 7%, n=5</td>
<td>40% ± 6%, n=6</td>
<td>40.0%</td>
<td>22.0%</td>
<td>58.0%</td>
<td>2</td>
<td>12</td>
<td>2</td>
<td>3</td>
<td>1</td>
<td>5</td>
</tr>
<tr class="even">
<td>Gai2 (3)</td>
<td colspan="2">83% ± 5%, n=4</td>
<td>73% ± 6%, n=5</td>
<td>73.0%</td>
<td>55.0%</td>
<td>91.0%</td>
<td>-</td>
<td>5</td>
<td>-</td>
<td>5</td>
<td>-</td>
<td>7</td>
</tr>
<tr class="odd">
<td>Gaq (3)</td>
<td colspan="2">70% ± 8%, n=7</td>
<td>66% ± 23%, n=2</td>
<td>66.0%</td>
<td>0.0%</td>
<td>95.0%</td>
<td>-</td>
<td>3</td>
<td>-</td>
<td>1</td>
<td>-</td>
<td>3</td>
</tr>
<tr class="even">
<td>PLCb3 (1)</td>
<td colspan="2">-</td>
<td>83% ± 15%, n=3</td>
<td>83.0%</td>
<td>38.0%</td>
<td>100.0%</td>
<td>-</td>
<td>3</td>
<td>-</td>
<td>-</td>
<td>-</td>
<td>3</td>
</tr>
<tr class="odd">
<td>PLCb4 (1)</td>
<td colspan="2">87% ± 6%, n=5</td>
<td>-</td>
<td>87.0%</td>
<td>69.0%</td>
<td>100.0%</td>
<td>-</td>
<td>4</td>
<td>-</td>
<td>4</td>
<td>-</td>
<td>4</td>
</tr>
</tbody>
</table>

#  References