Modeling of G-protein-coupled Receptor Signaling Pathways*

G-protein-coupled receptors (GPCRs)2 are the largest family of cell membrane receptors. An estimated 50% of current pharmaceuticals target GPCRs (1), suggesting that further increases in our understanding ofGPCRs and the signaling pathways they initiate will lead to new drug targets. Mathematical and computationalmodeling (here, simply “modeling”) has a substantial history in modern biology and pharmacology (2, 3) and offers a powerful tool for examining GPCR pathways. Such models can be used to better understand hypothesized mechanisms, run virtual (in silico) experiments, interpret data, suggest new drug targets, motivate experiments, and offer new explanations for observed phenomena.

G-protein-coupled receptors (GPCRs) 2 are the largest family of cell membrane receptors. An estimated 50% of current pharmaceuticals target GPCRs (1), suggesting that further increases in our understanding of GPCRs and the signaling pathways they initiate will lead to new drug targets. Mathematical and computational modeling (here, simply "modeling") has a substantial history in modern biology and pharmacology (2,3) and offers a powerful tool for examining GPCR pathways. Such models can be used to better understand hypothesized mechanisms, run virtual (in silico) experiments, interpret data, suggest new drug targets, motivate experiments, and offer new explanations for observed phenomena.

Many Simultaneous Kinetic Processes
The more we learn about GPCRs and the pathways they activate, the more complicated our picture of GPCR signaling becomes. At the subsecond to minute time scale, ligand binding, interactions of receptors with G-proteins, G-protein activation/deactivation, and the action of RGS (regulator of G-protein signaling) proteins in GAP or non-GAP roles occur. Depending on the particular G-protein subunit, many downstream signaling molecules (e.g. adenylyl cyclase, cAMP, phospholipase C, Ca 2ϩ , membrane channels) are transiently modulated, locally or over the entire cell. At a slightly longer time scale, receptor phosphorylation, arrestin binding, and activation of non-G-protein-dependent signaling pathways occur (4,5). Receptor trafficking events, i.e. internalization, recycling, routing to lysosomes, and up-regulation, as well as new receptor synthesis and regulation of gene expression, also occur over the time scale of minutes to hours. Many of the receptor and membrane level events are shown schematically in Fig. 1A. Noble attempts to consolidate information on intracellular signaling pathways initiated by GPCR binding are available (and evolving) (e.g. Database of Science Signaling) (6,7).
Put simply, it is difficult to intuit the net result of so many simultaneous kinetic processes. Add nonlinearities such as feedback, time-varying sequestering of molecules via scaffolds and other mechanisms, and multiple receptor and G-protein species to the temporal and spatial variations already mentioned, and it is virtually impossible. Modeling helps: putting in hypothesized mechanisms and numbers (rates, concentrations) allows both qualitative and quantitative insights. With models one can ask questions such as, which of several simultaneous pathways plays a larger role? Do these events happen fast enough to be important in signaling? What mechanisms might allow modulation of signaling, desensitization, or receptor cross-talk? What factors influence ligand efficacy? Interruption of processes at which points (i.e. drug targets) is most effective? Which experimental protocol is most likely to emphasize a particular mechanism?

Model Development
How do you create a model? Mathematical/computational modeling that is mechanistic (as opposed to empirical models, which approximate the shape of a relationship without any mechanistic basis, or statistical models) generally incorporates three steps. In the first, you need to turn the biology into rules (e.g."if two receptors collide, there is a 50% probability of dimerization") or equations. This means that you have to decide what is important and what is not, at least for the question(s) you are asking. Models may well leave out something you know happens, perhaps because it is not significant at the time or length scale you are examining or is something you would rather not focus on now. In the second step, you solve the equations or run the simulation governed by these rules. The question being asked is simply, "If all these things happen as I have written, i.e. with these particular rates or rules, what is the outcome?" In the third step, you use the solution of the model to say something meaningful about the biology. If the equations/rules give behavior consistent with the outcome of an experiment, this suggests that the model captures the relevant biological mechanism(s) behind the observation. You may be able to fit key model parameters, offering quantitative insights as to the rates of particular processes or the relative importance of a particular pathway.
Models can vary significantly in the level of detail used to describe molecular interactions. For example, models of GPCR signaling may include vast detail in the G-protein activation/ deactivation cycle (8) or may have no explicit inclusion of G-proteins at all but rather lump their effect into what happens to a downstream component. In other words, there is a choice as to how fine-or coarse-grained you make your model. More abstract models, ones that perhaps only identify molecules as signaling components A, B, and C and attempt to understand how the rules governing their behavior influence responses, are also possible (9). The appropriate form depends on the question you are asking and also on the data you have for model validation and testing. Modeling is meant to be an iterative process with experimentation, each one driving the other.
Typically, a minimal model is constructed and then grows in complexity, driven by new hypotheses and new data. Simple is good. There is no glory to be had by constructing a complicated model when a simple one can elucidate the key features of the biology; in fact, a complicated model may obscure a simple * This is the first article of six in the Thematic Minireview Series on Computational Biochemistry: Systems Biology. This minireview will be reprinted in the 2009 Minireview Compendium, which will be available in January, 2010. 1 To whom correspondence should be addressed. E-mail: linderma@ umich.edu. 2 The abbreviations used are: GPCRs, G-protein-coupled receptors; ODEs, ordinary differential equations; PDEs, partial differential equations; FRET, fluorescence resonance energy transfer; GRKs, GPCR kinases; PKA, cAMPdependent protein kinase.
result. (Indeed, with a complicated model, one can use reducedorder modeling to figure out which pathways are the most important, essentially approaching the simplicity goal from the other side (10).) Although models typically involve numbers, often the numbers themselves are less important than the qualitative insight (e.g. how RGS proteins shift the dose-response curve).
Determining model parameters such as rate constants is an important and often time-consuming step. Different approaches are used and are not mutually exclusive. Focus on a particular time scale in experiments may remove some parameters from consideration to allow remaining parameters to be more accurately determined (e.g. short time scale experiments will not be confounded by receptor recycling); inhibitors of particular pathways are also used. A complementary approach is to fit all parameters at once, e.g. using regression analysis or Monte Carlo techniques (11). Techniques for discrimination between models that appear to fit the data reasonably well and for systematic assessment of the impact of parameter uncertainties are available (10 -12) and will be increasingly important as more components of signaling pathways are identified and more detailed models become possible.

Types of GPCR Pathway Models
The pharmacological literature has a rich history in modeling the equilibrium states of GPCRs, models composed of algebraic equations that describe how the addition of ligand and/or G-protein will change the distribution of receptor states (13). Rates of processes are not included. However, equilibrium models have limited applicability in understanding the inherently kinetic process of signaling, especially downstream of G-protein activation, and predictions of kinetic and equilib-rium models with similar parameter values can be markedly different (14).
Kinetic models of GPCR signaling aim to link the time course of GPCR binding and other receptor level events with the kinetics of early (e.g. G-protein activation) and later (e.g. cAMP and Ca 2ϩ dynamics, desensitization, mitogen-activated protein kinase (MAPK) activation) downstream events. Many of these models are formulated as ordinary differential equations (ODEs), describing how the concentration (or number) of each particular species evolves with time. Examples of how to write such equations for receptor systems are given in Ref. 15 and typically involve mass action kinetics.
Because signaling propagates from GPCRs at the cell membrane into the cell and because molecules such as Ca 2ϩ may be released from discrete sources within the cell, there are likely to be spatial gradients of signaling molecules. To follow both spatial and temporal information during signaling, partial differential equations (PDEs), which arise naturally from the diffusion equation, can be used. Readily available software to solve ODEs and/or PDEs includes Mathematica (Wolfram Research, Inc.), MATLAB (The MathWorks, Inc.), and Berkeley Madonna TM . For solving equations in a more cell-specific framework, one excellent tool is Virtual Cell (16).
Stochastic models allow events to be described by probabilities; unlike the (deterministic) models above, given the same input, a stochastic model will give (somewhat) different results each time. Such models (equation-or rule-based) of GPCR signaling are used when there are small numbers of molecules involved or when it is useful to track likely movements of individual molecules (17,18). A disadvantage of stochastic models is that they are more difficult to solve, and simulations must be FIGURE 1. A, many of the receptor and membrane level events in GPCR signaling. RK, receptor kinase. B, relative importance of signal attenuation via GRK-mediated phosphorylation/arrestin binding or PKA-mediated phosphorylation following isoproterenol binding to endogenous ␤ 2 -adrenergic receptors on HEK-293 cells (27). C, cubic ternary complex model, an equilibrium model incorporating differing abilities of ligands to stabilize the active receptor state (R*) (13). D, depletion of "activate-able" G-proteins (empty ovals) from the vicinity of long-lived receptor-ligand complexes (44,45). Activated G-proteins are shown as filled ovals. E, kinetic scaffolding. Increased GTP hydrolysis gives a spatial focusing of the G␣ GTP signal (52). F, clusters of GPCRs resulting from dimerization and diffusion kinetics (18). Receptors are shown as filled (dimers) or empty (monomers) circles.
run multiple times to gather statistics on the range of possible outcomes. Tolkovsky and Levitzki (19) offered one of the earliest kinetic models of GPCR signaling, the collision coupling model, suggesting that ligand-bound receptors and enzymes (here, G-proteins) collide and couple transiently to produce enzyme activation. Model behavior was qualitatively consistent with their data on cAMP accumulation during binding of epinephrine to ␤-adrenergic receptors. Many of the models described below rely on collision coupling. Several different types of models suggest that diffusion is fast enough to allow a collision coupling mechanism to produce observed rates of G-protein activation and/or downstream signaling (15, 20 -22), although mechanisms to limit contact between molecules (e.g. lipid rafts) are not ruled out. GTP hydrolysis may be rate-limiting as compared with the steps of GTP binding and GDP release, and so those latter steps are not included in many models (23,24).

Insights into GPCR Pathways, Dynamics, and Rate-limiting Steps
To what extent do G-protein activation/deactivation dynamics modulate downstream responses? In the mating pheromone signaling pathway in Saccharomyces cerevisiae, arguably the most biochemically well characterized eukaryotic signaling pathway, two recent modeling efforts support a central role for G-protein dynamics. Using an ODE model of the receptor/ligand binding dynamics and the G-protein cycle, Hao et al. (24) argued for pheromone-dependent transcriptional induction of Sst2 (which functions as an RGS protein) as a negative feedback mechanism that leads to desensitization as well as a positive feedback mechanism that leads to pheromone-dependent loss of Sst2 protein. Experiments driven by that model indeed showed that Sst2 is ubiquitinated and degraded. Examining the same signaling pathway, Yi et al. (25) used an ODE model together with fluorescence resonance energy transfer (FRET) measurements reporting the association state of the G-protein heterotrimer following ligand binding. Changing the rate of G-protein deactivation via deletion of Sst2 was found to shift the dose-response curve to higher sensitivity; Sst2 was determined to increase the rate constant for G-protein deactivation by ϳ25-fold. These results suggest that RGS proteins and proteins that down-regulate receptor activity (e.g. GPCR kinases (GRKs)) should be useful drug targets (26). Bornheimer et al. (8) examined the role of RGS proteins in a kinetic network in more detail, allowing for a ternary complex of receptor, G-protein, and RGS protein and fitting GTP hydrolysis data from vesicle preparations containing the m1 muscarinic acetylcholine receptor, G q , and RGS4. The model demonstrated different behaviors or "signaling phenotypes" that might be observed depending on local concentrations of GPCR, G-protein, and RGS. One of these phenotypes, in which G-proteins and receptors are clustered, offers a mechanistic hypothesis for the previously unexplained and unexpected observation that, in some systems, RGS proteins do not change the maximal response, only the kinetics.
Receptor dynamics, including synthesis and desensitization kinetics, also play important roles in determining response dynamics in GPCR pathways. Yi et al. (25) found that both receptor synthesis and ligand-induced receptor endocytosis contributed to longer term features in the yeast mating response. Violin et al. (27) used a novel live cell cAMP sensor to follow the dynamic response to isoproterenol binding to endogenous ␤ 2 -adrenergic receptors on HEK-293 cells and fit the data with an ODE model. They determined that, at least for the conditions/system tested, receptor inactivation due to GRKmediated phosphorylation followed by binding to ␤-arrestin is much more significant than receptor inactivation due to cAMPdependent protein kinase (PKA)-mediated phosphorylation (Fig. 1B); enhanced cAMP clearance by PKA activation of PDEs also played an important role in limiting the response. The model suggests the next round of experiments, i.e. whether the relative importance of the different desensitization pathways changes with receptor concentration (because the GRK/arrestin pathway should saturate) or the particular agonist or agonist concentration used, and may give clues as to relevant drug targets to affect desensitization.
When are spatial gradients in concentrations important? Models of chemotaxis suggest that such gradients are very important, and a central hypothesis in the field (local excitation, global inhibition) requires that cell components are not "well mixed" (28). Narang (29) developed a model for eukaryotic gradient sensing that relies on spatial segregation of distinct signaling molecules at the front and rear of the cell. G i activation plays a key role for the "frontness" pathway (the actin-rich front of the cell), whereas G 12/13 activation generates "backness" signals (controlling the morphology of the cell rear). The two pathways inhibit each other. Alternatively, Levine et al. (9) suggested a model involving a membrane-bound activator and diffusing inhibitor, which they postulated could be G␣ and G␤␥, respectively. These models make predictions for the distributions of particular signaling molecules within the cell as a function of precise ligand (chemoattractant) gradients that may now be achievable using microfluidic techniques (30).
Spatial gradients also play roles in the visual phototransduction cascade, the most quantitatively characterized G-protein signaling system. Building on early insights in phototransduction modeling (20), Bisegna et al. (31) developed a detailed model incorporating structural features of the rod outer segment of vertebrates as well as spatial and temporal modeling of second messenger generation. Diffusion of both cGMP and Ca 2ϩ during signaling is predicted to be important for the suppression of variability in the single photon response. The need for spatial detail in models is also argued by Saucerman et al. (32) and Iancu et al. (33) for cardiac myocytes.
Models such as these that quantitatively describe the dynamics of GPCR signaling pathways have a tremendous potential in terms of analyzing alternative therapeutic agents to affect activation and/or desensitization. For example, Saucerman et al. (34) used an integrated ODE model of ␤-adrenergic receptor signaling and excitation-contraction coupling in cardiac myocytes to compare potential experimental and pharmacological treatments (e.g. receptor or adenylyl cyclase overexpression for the ability to maximize isoproterenol response while maintaining physiological basal cAMP levels).
Tight integration of mathematical/computational modeling with live cell, high time resolution imaging of ligand binding and signaling pathways will hasten the evolution of more complete, accurate, and useful models. For GPCR pathways, recently developed probes include a fluorescent reporter of PKA-mediated phosphorylation (32), a FRET-based sensor for cAMP (27), FRET-based sensors that report association of G␣ and G␤␥ (35), and a bioluminescence resonance energy transfer-based assay to monitor arrestin-adaptin interactions (36). The integration of GPCR pathway models and experiments in describing calcium dynamics (for which probes have been available for some time) has also been described (37)(38)(39).

Factors That Influence Efficacy
All GPCR ligands for a particular receptor are not the same in terms of pathway activation. Even after correcting for differences in binding affinity, ligands have various abilities to induce responses (i.e. efficacy) and can be classified as full and partial agonists, antagonists, and inverse agonists. Certainly, the details of exactly how (and how fast) a ligand influences GPCR conformation matter (40,41), but how does this play out in terms of an analysis of the signaling pathway, and what role can modeling play in deciphering this?
One explanation is that ligands have differing abilities to stabilize active receptor states. Building on the ternary complex model, an equilibrium model that describes the reversible association of ligand, receptor, and G-protein, the extended ternary complex model and the more thermodynamically complete cubic ternary complex model include both active and inactive receptor states (R* and R) in terms of ability to activate G-proteins. Differences in ligand efficacy are postulated to be due to differences in the ability of the ligand to bias the receptor into the active state (Fig. 1C) (13). Ligands that favor the inactive state are inverse agonists; those that favor the active state are agonists. Although the models just mentioned are equilibrium models, a recent kinetic model by Kinzer-Ursem et al. (42) uses kinetic data on the binding of several N-formyl peptides to the N-formyl peptide receptor on neutrophils to demonstrate that a receptor state not discerned from kinetic binding experiments is necessary to account for the oxidant production and actin polymerization response; this state is interpreted as the active receptor state, and the model provides a way to measure the relative strength of a ligand to bias the receptor into this active conformation.
Using the cubic ternary complex model or a kinetic version of it, one can make an interesting prediction: the existence of ligands that can function as either positive agonists or inverse agonists depending on cell-and ligand-specific properties ("protean agonists") (12,43). These predictions are qualitatively consistent with the (scant) experimental data on protean agonism to date and hint at novel strategies for manipulating cell behavior.
Additional factors that influence efficacy have also been hypothesized and supported with models. When the collision coupling model is implemented in a framework that allows for spatial gradients of membrane components to be tracked, e.g. via simulations of individual receptors and G-proteins, one finds that when comparing cases of equal occupancy of receptors by agonist, increased movement of the agonist among receptors (shorter half-life of the receptor-agonist complex) results in increased G-protein activation, as receptors newly occupied by agonist have access to G-protein close by (Fig. 1D) (44,45). Such movement has also been predicted to partially decouple G-protein activation from receptor phosphorylation (46).
The recent elucidation of non-G-protein-dependent signaling by GPCRs adds a new dimension to understanding efficacy because ligands very capable of signaling via one pathway may be less capable at signaling via the other. Recently, Drake et al. (47) examined a broad range of ligands for the ␤ 2 -adrenergic receptor; many had similar efficacies for both pathways, but a few showed a bias for the ␤-arrestin-associated pathway. The development of models that account for these additional receptor states and pathways may aid in the development of novel therapeutics (48,49).

New Ideas on Cross-talk and Dimerization
Cross-talk can be broadly defined as the influence of activation of one receptor type on signaling through a second receptor type, i.e. a lack of signaling specificity. When multiple GPCR species can activate the same G-proteins, the collision coupling model predicts that there is competition between receptors for the same G-protein pool (50). Such competition would be one route to cross-talk. However, cross-talk in GPCR systems may in many cases be limited. This could be achieved by limiting contact of receptors and/or effectors through physical mechanisms (e.g. scaffolds, lipid rafts) (51). Several models also shed light on additional mechanisms that may limit cross-talk. Bornheimer et al. (8) suggested that clustering/coupling interactions between receptor, G-protein, and/or RGS proteins could offer such physical scaffolding. Alternatively, Zhong et al. (52) showed that the membrane level "spread" of G-protein signaling could be limited by a kinetic scaffolding mechanism in which RGS proteins, by their acceleration of GTP hydrolysis, reduce depletion of local (near the receptor) G␣ GTP levels enough to allow rapid recoupling of G-protein to receptor. This would maintain G-protein activation near the receptor but deplete it farther away (Fig. 1E), a spatial focusing of the signal.
Why do GPCRs dimerize? GPCRs may associate as homodimers and heterodimers and possibly as higher oligomers in the membrane (53,54). Such associations may be transient and altered by ligand binding. There is debate over whether such oligomerization is essential for signaling, although it has been recently demonstrated that both monomeric rhodopsin and ␤-adrenergic receptors incorporated into phospholipid particles efficiently activate G-proteins (55,56). The relevance of this finding for other GPCRs, for GPCRs in vivo, and for non-G-protein-dependent pathways is still an open question. Modeling has suggested another way to think about dimerization (18). When both association and dissociation of dimers are relatively fast (i.e. dimers are not high affinity) and diffusion is relatively slow, larger clusters of receptors can form (Fig. 1F). Such dimerization-induced clustering thus affects the overall organization of receptors in the membrane and is rapid enough to play a role in signaling. Because organization may affect access of receptors to signaling molecules (both positively, e.g. making it easier to "share" an effector, and negatively, e.g. making it more difficult for an effector to reach a receptor) and because a cluster of receptors could essentially act as a multivalent ligand for downstream players, this suggests that changes in dimerization may be used to modulate information flow through signaling pathways as well as cross-talk.
Finally, although not the focus of this review, a number of studies have explored the structure and/or kinetics of signaling pathways in a more general framework and predicted behaviors, including noise filtering, ultrasensitivity, oscillations, bistability, and robustness, that may have application for the modulation of signaling in GPCR pathways (57)(58)(59).

Future Prospects
Modeling is likely to become an increasingly useful tool for understanding GPCR pathways. The United States Food and Drug Administration Critical Path Initiative has recently identified model-based drug development, including drug and disease modeling, as an important goal (www.fda.gov/oc/ initiatives/criticalpath). Powerful and increasingly userfriendly software is available. Interesting new hypotheses generated by modelers are being published in prominent biochemistry, biology, and biophysics journals. In addition, new experimental techniques, particularly RNA silencing of pathway molecules and novel fluorescent probes that allow for single cell kinetic data, are available to both generate data for model building and allow testing of key model findings.