Both Ligand- and Cell-Specific Parameters Control Ligand Agonism in a Kinetic Model of G Protein–Coupled Receptor Signaling

G protein–coupled receptors (GPCRs) exist in multiple dynamic states (e.g., ligand-bound, inactive, G protein–coupled) that influence G protein activation and ultimately response generation. In quantitative models of GPCR signaling that incorporate these varied states, parameter values are often uncharacterized or varied over large ranges, making identification of important parameters and signaling outcomes difficult to intuit. Here we identify the ligand- and cell-specific parameters that are important determinants of cell-response behavior in a dynamic model of GPCR signaling using parameter variation and sensitivity analysis. The character of response (i.e., positive/neutral/inverse agonism) is, not surprisingly, significantly influenced by a ligand's ability to bias the receptor into an active conformation. We also find that several cell-specific parameters, including the ratio of active to inactive receptor species, the rate constant for G protein activation, and expression levels of receptors and G proteins also dramatically influence agonism. Expressing either receptor or G protein in numbers several fold above or below endogenous levels may result in system behavior inconsistent with that measured in endogenous systems. Finally, small variations in cell-specific parameters identified by sensitivity analysis as significant determinants of response behavior are found to change ligand-induced responses from positive to negative, a phenomenon termed protean agonism. Our findings offer an explanation for protean agonism reported in β2--adrenergic and α2A-adrenergic receptor systems.


Introduction
G protein-coupled receptors (GPCRs) are the largest class of cell membrane receptors with almost 2,000 members identified [1]. It is estimated that more than 50% of pharmaceuticals target GPCRs [2]. While the majority of pharmacologic research has focused on ligand-specific properties that influence cell behavior, relatively few studies focus on cell-specific parameters that may also determine cell responses [3][4][5]. Studying the effect of changes in both ligandand cell-specific parameters on cellular behavior is complicated by the large number of interactions and feedback mechanisms inherent in cellular signaling. Thus it becomes necessary to use quantitative models to aid in the analysis of these systems.
Typical models of GPCR signaling are termed ternary complex models or TCMs (reviewed in [6][7][8]). These models feature ligand (L) binding to receptor (R) to form a ligandreceptor complex (LR), and LR interaction with G protein (G) to form the ternary ligand-receptor-G protein (LRG) complex. Subsequent equilibrium models of GPCR signaling have remained true to this paradigm while incorporating additional receptor (e.g., active receptor R* or inactive receptor R) or G protein states or other effectors to account for experimental findings [9][10][11][12][13][14][15]. A key feature of these models is that active receptors can associate with G protein in the absence or presence of ligand to form R*G and LR*G, respectively, and both these complexes can signal. Kenakin and colleagues have proposed a thermodynamically complete representation of ligand, receptor, and G protein interac-tions termed the cubic ternary complex model (cTCM) and shown in Figure 1 [11]. In this model, inactive receptor can both bind ligand (to form LR) and associate with G protein (to form RG or LRG).
TCMs are typically equilibrium models and while they have been widely used, it is well known that kinetic models are better able to replicate the intrinsic dynamics of signal transduction, as has been discussed previously (see [14,[16][17][18][19][20]). Furthermore, predictions of kinetic and equilibrium models with similar parameter values can be markedly different [21,22]. Indeed, a number of groups have discussed the importance of kinetics in analyzing GPCR systems [8,20,[23][24][25][26][27]. Our group has thus developed a kinetic version of the cTCM termed the cubic ternary complex activation model (cTCAM, Figure 1) [22]. The cTCAM incorporates a G protein activation feedback loop whereby G proteins (G) couple to and are activated by active receptors (R* and LR*), allowing for GTP binding and uncoupling of the G protein heterotrimer into a and bc subunits. The GTP on active G proteins (Ga GTP ) is hydrolyzed by the intrinsic GTPase activity of the alpha subunit (with or without participation of regulator of G protein signaling (RGS) proteins) to form inactive Ga GDP subunits. The feedback loop is completed when the Ga GDP and Gbc subunits couple to reform inactive G protein (G). We found that the predictions of the kinetic model (cTCAM) can be strikingly different than those of the equilibrium model (cTCM) in terms of the character of the response (positive/neutral/inverse agonism), suggesting the importance of using the more realistic dynamic model [22].
In this work we use the cTCAM as a dynamic model of the initial events in GPCR signaling to identify both the ligandand cell-specific parameters that may be key determinants of the character of cellular response. One drawback of the kinetic cTCAM, and indeed for most if not all signal transduction models, is that many parameter values have not been measured and others may vary over several orders of magnitude [10,12,22,25,28]. Additionally, the large number of parameters and incorporation of the G protein activation feedback loop make intuition of model behavior difficult. Thus it becomes necessary to introduce techniques for efficient sampling of the input parameter space and quantification of model output. In the risk analysis and environmental engineering fields, uncertainty and sensitivity analysis have been routinely used to sample input parameter space and identify key parameters [29,30]. These techniques have recently been introduced into the biological sciences. In particular, Latin hypercube sampling (LHS) with partial rank correlation coefficients (PRCC) has been used to perform uncertainty and sensitivity analysis in epidemiological studies of HIV and tuberculosis [31,32], to understand the dynamics of tuberculosis infection and immunity in host-pathogen models [33,34], and to analyze parameter sensitivity in a T cell receptor activated Erk-MAPK signaling pathway [35]. Additionally, several studies have used genetic algorithm-based search methods to perform parameter fitting and sensitivity analysis for models of T cell receptor and GPCR activation [12,35,36]. Rundell and colleagues found that these global analysis methods (LHS/PRCC and genetic algorithm approaches) and others (Sobol's method and Fournier amplitude sensitivity test (FAST)) give very similar results [35].
LHS has been shown to be a computationally efficient method for sampling parameter ranges. It is more than one order of magnitude more efficient than random sampling methods [29,30,37]. Additionally, statistical techniques can be used to identify parameters that are most important in determining output variables. Correlation coefficients can be readily calculated to identify parameters whose variation is strongly correlated with variations in an output parameter of interest. For nonlinear monotonic systems such as the cTCAM, PRCC is known to be the most appropriate [29,38,39]. PRCC values can be calculated at each time point of the simulation, and the relative importance of the parameters can be tracked over time. Here we use LHS and PRCC to identify parameters that are important determinants of G protein activation in a general model of GPCR signaling (cTCAM, Figure 1). In particular, we are interested in how small variations in parameter values might give rise to large differences in the character (positive/neutral/inverse agonism) of ligand-induced responses.
Different ligands, while binding to the same receptors, are often able to transmit different levels of signal per bound receptor and thus have different levels of response. Recent studies have shown that the same ligand may not only induce varying levels of response but also both positive and negative responses in some GPCR systems. Here we refer to the ability of a ligand to act both as a positive agonist and an inverse agonist as protean agonism (a term previously introduced by Kenakin [40]). Protean agonists (ligands that are able to induce protean agonism) have been identified in several GPCR systems, including b 2 -adrenergic, a 2A -adrenergic, opioid, and histamine H3 receptors [41][42][43][44]. After identification of ligand-and cell-specific parameters that play critical roles in determining the character of a response, we then focus on two particularly interesting studies of b 2 -adrenergic and a 2A -adrenergic receptors in which both positive and inverse agonism are produced by introduction of the same ligand to the system. Chidiac and colleagues reported that the b 2 -adrenergic receptor partial agonist dichloroisoproterenol (DCI) acted as both a positive and an inverse agonist in Sf9 cells overexpressing b 2 -adrenergic receptor despite the fact that the treatment of the cells was similar in both cases [41]. Jansson and colleagues reported that the a 2A -adrenergic ligand levomedetomidine (levomed) had opposite effects on cAMP production in different cell lines, activating the receptor in several systems (S115 and PC10 cells) while inhibiting constitutive activity of the endogenous receptor in HEL 92.1.7 cells [42,45,46]. We find that changes in cellspecific parameters identified as key determinants of response behavior by sensitivity analysis are consistent with these seemingly contradictory behaviors in the b 2 -and a 2Aadrenergic receptor systems.

Author Summary
G protein-coupled receptors (GPCRs) are transmembrane proteins involved in physiological functions ranging from vasodilation and immune response to memory. The binding of both endogenous ligands (e.g., hormones, neurotransmitters) and exogenous ligands (e.g., pharmaceuticals) to these receptors initiates intracellular events that ultimately lead to cell responses. We describe a dynamic model for G protein activation, an immediate outcome of GPCR signaling, and use it together with efficient parameter variation and sensitivity analysis techniques to identify the key cell-and ligandspecific parameters that influence G protein activation. Our results show that although ligand-specific parameters do strongly influence cell response (either causing increases or decreases in G protein activation), cellular parameters may also dictate the magnitude and direction of G protein activation. We apply our findings to describe how protean agonism, a phenomenon in which the same ligand may induce both positive and negative responses, may result from changes in cell-specific parameters. These findings may be used to understand the molecular basis of different responses of cell types and tissues to pharmacological treatment. In addition, these methods may be applied generally to models of cellular signaling and will help guide experimental resources toward further characterization of the key parameters in these networks.

Cell-Specific Parameters Are Highly Correlated with Variation in Response Characteristics
LHS was used to efficiently sample parameter values from the ranges listed in Table 1. One LHS simulation typically sampled each of the 16 parameters 1,000 times, producing 1,000 solutions to the model equations (Text S1). The formation of Ga GTP is computed (see Equation S.27 in Text S1). Eight example solutions of the time course of Ga GTP formation are shown in Figure 2A. As expected, variations in parameter values caused large variations in response behavior. For example, for the cases shown in Figure 2A, at steady state prior to ligand addition (time zero), Ga GTP values varied between approximately two and 240 per cell. Once ligand is added to the system, it binds to receptors which in turn bind to G protein, changing the distribution of receptors and G proteins among their various states and causing increases or decreases in G protein activation (as seen by changes in Ga GTP formation, Figure 2A). In some instances, G protein activation increases rapidly, and then falls due to GTP hydrolysis. Ligand-induced responses occur on the order of 5 s to 30 s at this ligand concentration, which is approximately the timescale over which G protein activation is known to occur [47][48][49]. Calculating the percent change of Ga GTP upon ligand addition (designated %OverBasal) allowed for easier observation of decreases in G protein activity upon ligand binding; in Figure 2B these are seen as negative %OverBasal values. In this way, responses were directly related to the pharmacological classifications of ligand efficacy-positive (increase in %OverBasal), neutral (no change in %OverBasal), and inverse agonism (decrease in %OverBasal)-and thus can be compared with experimental data that normalize to control conditions.
To determine the correlation between parameter values and levels of G protein activation, PRCC values were calculated at 0.25-s intervals for varying ligand concentrations (0.1 nM to 0.1 mM). Correlations (PRCC values) for each parameter listed in Table 1 were calculated with respect to the two different measurements (outputs) of G protein activation discussed above, the number of Ga GTP and %OverBasal. Table 2 shows the rank order of PRCC values for the two responses at two time points, 5 s and 2.5 min, and two ligand concentrations, 1 nM and 10 lM, representing pre-and post-steady state time points and sub-and suprasaturating ligand conditions, respectively.  The cTCM (black) is a thermodynamically complete equilibrium representation of ligand (L), receptor (R), and G protein (G) interactions [11]. Association and dissociation of L and R is represented here from top to bottom, R and G interactions from front to back, and the interconversion of inactive and active R states from left to right of the cube. Based on the cTCM, the cTCAM (black and red) incorporates the dynamics of activation and recycling of G protein (dashed lines) into a kinetic model of LRG interactions [22]. A brief summary of model parameters is found in Table 1. Description of model parameters, assumptions, and equations are given in Text S1. doi:10.1371/journal.pcbi.0030006.g001 Previous analysis of the equilibrium cTCM in our group found that K act , G total , a, d, c, and K g all played a role in determining the character (positive/neutral/inverse agonism) of the response [20]. By using the cTCAM to include the necessary activation events, however, a somewhat different set of parameters was found to be highly correlated with response generation as summarized in Table 2. Not surprisingly, the ligand-specific parameter most correlated with both measures of response generation (instantaneous number of Ga GTP and %OverBasal) was the effectiveness with which the ligand induces an active receptor conformation (a). Indeed, this was the only ligand-specific parameter found to be highly correlated with the Ga GTP response. The character of responses was influenced by several cell-specific parameters including receptor and G protein expression (R total and G total ), parameters involved in the G protein activation loop (k Gact , k GTP , and k G ), the equilibrium ratio of active to inactive receptor (K act ), and the rate and efficiency of receptor-G protein coupling (b and k 11 ). Although generally the parameters that were highly correlated did not differ between the two response measures (Ga GTP and %OverBasal), the rank order of the PRCC values did vary. Similar results were found when the microscopic reversibility assumption of the model (discussed in the Methods section) was relaxed and all forward and reverse rate constants (shown in Figure S1) were varied (unpublished data).
The rank order of the parameters that were highly correlated with Ga GTP did not vary significantly with time, and thus PRCC values for the parameters at 2.5 min but not 5 s are shown in Table 2. Additionally, the PRCC values did not vary significantly at varying ligand concentrations, with the exception that k G and a were found to be lower in rank order at the lower ligand concentrations.
For the %OverBasal response, significant differences were seen between low and high ligand concentrations. At low ligand concentration (1 nM in Table 2), the parameters for reversible ligand binding (k 3 and K a ) were highly correlated with %OverBasal. Additionally, the PRCC values of both k 3 and K a were found to vary with time, and thus the PRCC values for all parameters at both 5 s and 2.5 min at this ligand concentration are listed in Table 2. As shown in Figure 3, the ligand association rate constant k 3 was found to be highly correlated at 5 s, but was not significantly correlated after 2.5 min of ligand binding. In contrast, the PRCC value for K a rapidly increased over the first 30 s of the simulation and was significantly correlated with %OverBasal at longer times (greater than 1 min). Thus at sub-saturating ligand concentrations in the first few seconds of ligand binding when response characteristics are likely to be determined, changes in response are highly sensitive to the ligand association rate constant. At high ligand concentration (10 lM in Table 2), k 3 and K a were not highly correlated with %OverBasal, nor was there a significant difference between PRCC values at the two times, and thus only PRCC values at 2.5 min are listed in Table 2.

Receptor and G Protein Levels Independently Affect Response Characteristics
The levels of receptor and G protein expression were both identified as important parameters in response generation ( Table 2). To investigate the impact that variations in these parameters have on the character of response generation, simulations varying R total and G total over their physiologic range were performed. To simulate (but not replicate) experiments that measure the accumulation of a signal over time, such as radioligand assays, the integral of Ga GTP over the first 10 s of ligand binding was calculated and normalized to basal values with no ligand present and is designated %Accum according to Equation 2. Qualitatively similar results are obtained using the instantaneous number of Ga GTP at either 5 s or 10 s after ligand binding (unpublished data). For certain sets of parameters values, particularly when a, d, and c 6 ¼ 1, the %Accum upon ligand addition may vary from positive to negative values. The explanation for this protean agonist behavior is as follows. In conditions of constitutive activity (L ¼ 0, Kact 6 ¼ 0), there exist two receptor-G protein complexes (RG and R*G) of which only R*G can activate G protein. Upon ligand addition, the number of receptor-G protein complexes increases to four (RG, R*G, LRG, and LR*G), two of which (R*G and LR*G) can activate G protein. At low values of R total and G total , the addition of ligand to the system can actually reduce G protein activation because of the redistribution of receptors and G proteins among their various states. When this occurs, the ligand behaves as an inverse agonist; a decrease in the %Accum upon ligand addition is then seen (Figure 4A and 4B). As the total numbers of receptors and G proteins are increased, the total number of R*G and LR*G generated upon ligand addition greatly outnumber those of R*G prior to ligand addition, and the ligand behaves as a positive agonist [50]. Figure 4A shows the behavior of %Accum when R total is held constant at 5,000/cell and G total is varied. Note that as G protein expression is increased, the same ligand is predicted to change from an inverse to a neutral to a positive agonist. Similar trends are seen when G total is held constant at 10,000/cell and R total is varied ( Figure 4B).
Previous studies have proposed that the ratio of receptors to G proteins may play an important role in determining the efficacy of a ligand [40,51]. Furthermore, previous analysis of the equilibrium cTCM has shown that relatively large changes (.25-fold) in K act and G total may change the efficacy of a ligand, from positive to neutral or neutral to negative responses [50]. To investigate the role of the R/G ratio in our more physiological model, %Accum at a saturating ligand concentration (10 lM) from Figure 4A and 4B was plotted versus the ratio of R total /G total in Figure 5. As seen in Figure 4, the level of response increased as either receptor or G protein levels were increased. However, because increasing G total decreases the ratio R/G, the response generated in simulations of constant R total and varying G total (dashed line in Figure 5) increases as R total /G total decreases. In contrast, as receptor expression is increased, the level of response and the ratio R total /G total increase as indicated by a positive slope in the curve (solid line in Figure 5). Thus, there is not a clear relationship between the ratio R total /G total , and this measure cannot be used as a straightforward predictor of response efficacy. These results indicate that the individual levels of receptor and G protein expression are important in determining ligand efficacy.
Our Results Offer New Explanations for Protean Agonism in the a 2A -and b 2 -Adrenergic Receptor Systems a 2A -Adrenergic receptors couple to Ga i proteins, activating the G proteins, which in turn inhibit adenylyl cyclase activation. Jansson and colleagues have reported that the a 2A -adrenergic ligand levomed is a positive agonist in PC10 cells, causing an inhibition of cAMP production as shown in Figure 6A [45]. However, the same group has reported that levomed also acts as an inverse agonist in HEL 92.1.7 cells, causing an increase in cAMP production ( Figure 6A) [42]. These results suggest that a parameter (or parameters) different between the two cell types may play a critical role in determining the character of the response, and cause protean agonism.
Guided by the results from our uncertainty and sensitivity analysis (Table 2), we tested whether changes in a single cellspecific parameter, a parameter that might differ between the PC10 and HEL 92.1.7 cells, could produce protean agonism. Small (less than or equal to half an order of magnitude) changes in any of three parameters, G total , k Gact , or K act , were all able to produce protean agonism, as shown with representative simulations in Figure 6B-6D. All three changes represent physiologically reasonable explanations for the protean agonism. Differences in G protein expression (G total ) between two cell lines and indeed within a cell line are common, although actual expression levels are rarely quantified [52][53][54]. Differences in the rate constant of G protein activation (k Gact ) between cell lines then may be due to differences in G protein isoform expression profiles [55]. Further, differences in the equilibrium ratio of active to inactive receptors (K act ) between cells expressing endogenous receptors and those with transfected receptors, in this case, between the transfected PC10 and endogenous HEL 92.1.7 cells studied by Jansson and colleagues, would also seem quite plausible.
Interestingly, although differences in a 2A -adrenergic receptor expression between HEL 92.1.7 cells (2,900-4,100 receptors at the cell membrane; [56]) and PC12 cells (about 5-  fold greater; [57]) have been reported, in our simulations these modest changes in R total did not produce protean agonism. Similarly, although the GTPase activity of the Ga protein subunit (k GTP ) can vary widely between different cell types depending on the expression of RGS proteins [3,25,58], variation of k GTP did not produce protean agonism in our model. Other cellular parameters identified by uncertainty and sensitivity analysis ( Table 2) also did not produce protean agonism, at least not for moderate (less than an order of magnitude) changes in their values and for physiologically reasonable values of the remaining parameters.
As a second example of protean agonism, Chidiac and colleagues have reported that the b 2 -adrenergic ligand DCI produces both stimulatory and inhibitory responses in Sf9 cells despite similar treatment of the cells ( Figure 7A) [41]. The parameter with the most likely variation in this system is G protein expression (G total ). While quantitative measurement of Ga s -like proteins in Sf9 cells has not to our knowledge been reported, several studies have shown that there may be wide variations in G protein expression in this system. Seifert et al. [52] and several references therein report that endogenous Ga s -like proteins could not be detected with immunoblot assays of Sf9 membranes, while Kleymann et al. [59] and Leopoldt et al. [54] have detected endogenous Ga slike proteins in immunoblots of membranes. Additionally, infection of Sf9 cells with baculoviruses has been shown to downregulate the expression of endogenous G proteins [54]. Simulations at varying G total showed that as little as a 5-fold change in G total was able to produce protean agonism in our model, as shown in Figure 7B.
Chidiac and colleagues further report that after treatment with isoproterenol (a b 2 -adrenergic agonist) DCI acts as an inverse agonist only ( Figure 7C) [41]. In other words, protean agonism by DCI is not seen after receptor desensitization. In another cell type (S49 lymphoma cells), Insel and colleagues have found that redistribution of Ga protein subunit between the membrane and cytosolic compartments occurs after prolonged treatment with isoproterenol [60]. The possibility that this loss of G proteins contributes to inverse agonism in Sf9 cells was tested by decreasing the number of available G proteins (G total ). Ligand stimulation of these desensitized cells (parameters as in Figure 7B but with half the number of G proteins) produced only inverse agonism in our model ( Figure 7D). Thus our analysis suggests the observation of positive or inverse agonism in these cells is sensitive to G protein number and that small changes in G protein number can account for observations of both protean agonism (before desensitization) and inverse agonism (after desensitization).

Discussion
As is common in large signal transduction networks, there are significant uncertainties in model parameter values, and it is difficult to intuit model behavior. Using efficient (A) Total G protein (G total ) is varied from 1,000,000/cell to 3,000 per cell, and R total ¼ 5,000/cell. (B) Total receptor (R total ) is varied from 30,000 to 1,000 per cell, and G total ¼ 10,000/cell. Parameter values: k 1 ¼ 1 s À1 , k 3 ¼ 1 3 10 7 M À1 s À1 , k 11 ¼1 3 10 À4 (number/cell) À1 s À1 , k Gact ¼5 s À1 , k GTP ¼1 s À1 , k G ¼1310 À4 (number/cell) À1 s À1 , doi:10.1371/journal.pcbi.0030006.g004 Figure 5. The Effect of the Ratio of Receptor to G Protein on G Protein Activation R total was set at 5,000/cell and G total varied (dashed line). As G protein expression increases, the response changes from negative to positive agonism. G total was set at 100,000 per cell and R total varied (solid line). As receptor expression increases, the response changes from negative to positive agonism. doi:10.1371/journal.pcbi.0030006.g005 parameter variation and sensitivity analysis methods, we identify parameters in a dynamic G protein activation model that play key roles in determining the level and indeed the character (positive/neutral/inverse agonism) of the response. This type of analysis allows for quick and efficient identification of important parameters. Drug design and development typically focus on altering ligand-specific parameters to produce desired cell response. Significantly, however, we find that not only ligand-specific parameters, but also cell-specific parameters, play critical roles in determining response behavior. To demonstrate how changes in cellular specific parameters can dramatically change cell response, we apply our findings towards studies of protean agonism in the a 2Aand b 2 -adrenergic receptor systems.
Our analysis shows that several cell-specific parameters not previously identified contribute significantly to the character of ligand-induced responses. These include the total receptor concentration (R total ), the GTP hydrolysis rate constant (k GTP ), and the G protein activation rate constant (k Gact ) (see Table 2 for a full list). Receptor number is known to vary widely in different cell types, and is regularly manipulated using transfection technologies. This is one of the few parameters in our model that is routinely quantified. It is well known that the GTPase activity of the Ga protein subunit may vary and is largely dependent on the presence of RGS proteins in the system [61,62]. Therefore, it is not surprising that the rate constant for GTP hydrolysis (k GTP ) would be found to be important in response generation; indeed much attention has been given to RGS proteins as new therapeutic treatments [3,63,64]. Our findings are also consistent with an elegant study by Bornheimer and colleagues that analyzed G protein activation by active receptor and deactivation by GTPase activation proteins (GAPs) [12]. They found that local concentrations of receptors and GTPase activation proteins mediate various regimes of response behavior by kinetically controlling G protein activity. Following its identification by sensitivity analysis, we find that small changes in k Gact (the rate constant for G protein activation) may result in protean agonism such as reported by Jansson and colleagues for the a 2A -adrenergic system [42,45].
G protein concentration and K act , the ratio of active to inactive receptor states, are two cell-dependent parameters identified by our study (using a dynamic model of G protein activation) and two previous studies (using equilibrium models) as parameters that are key to determining the character of response and, when varied, could cause protean agonism (Table 2 and Figures 6B, 6D, and 7B) [20,40,50]. Overexpression or underexpression of receptor and G protein or expression of different isoforms of G proteins may give very different results than those found in endogenous systems. Additionally, G protein expression may vary largely in cells depending on a variety of factors, including cell type, state of cell development [65], prior treatments to the cells [54,60,66,67], and disease state [68][69][70]. Although K act has not been rigorously quantified, receptor mutation and fluorescent imaging studies are lending valuable insight into how conformational changes in the receptor confer activity with and without ligand stimulus [71][72][73][74][75].
Our findings on the sensitivity of G protein activation to these parameters are corroborated by experimental studies implicating changes in G protein and receptor expression and receptor activity in cardiac and Alzheimer disease. For example, it has been shown that both b 1 -and b 2 -adrenergic receptors are expressed at decreased levels in studies of cardiac failure, while Ga i subunits are increased [76]. As a second example, although muscarinic receptor density is not changed in brains of Alzheimer patients, the functionality of the receptors has been shown to be compromised (suggesting an inactive receptor state and changes in K act ), and decreases in the function of Ga q protein have also been reported [77]. Thus, both experimental reports and our modeling results suggest that small changes in these key parameters may disrupt normal signaling and lead to disease states.
To account for day-to-day and cell-to-cell variability, experimental results are almost always presented normalized to basal levels. However, clearly some information is lost when normalizing, both in experimental data and in modeling studies. For example, in experimental systems normalizing washes out variations in the ''basal state'' of the cells that may be an indicator of the current state of the cell or of future response characteristics. In the context of our model, un-normalized cell population or single cell data would allow for a more in-depth study of how parameters contribute to both basal and post-ligand treatment responses.
Our findings represent what may be only a small sampling of cellular parameters that influence ligand efficacy. Indeed, small variations in combinations of parameters could also lead to observations of protean agonism. For the sake of conciseness, these are not analyzed here. However, the methods introduced here provide a computationally efficient mechanism to begin to explore all possibilities. One limitation of this type of analysis, particularly PRCC analysis, is that it can only identify trends in certain directions within the given parameter space. Other analysis techniques for LHS sampling, such as subjective, differential sensitivity analysis, one-at-a-time design, and the adjoint method have been used, although these also have significant limitations [39]. Another limitation for using PRCC is that the system must be monotonic; however, this is easily checked by monitoring scatterplots, and there exist other global analysis methods by which to analyze the impact of parameter variation on model output. Recent work has compared multiple global analysis methods and found that LHS/PRCC, genetic algorithm approaches, Sobol's method, and Fournier amplitude sensitivity test (FAST) give very similar results [35]. Finally, both kinetic measurements and modeling will be important to our progress in understanding and manipulation of G protein activation. New technologies, such as those reported in yeast and HEK 293 cells using fluorescence resonance energy transfer (FRET) and bioluminescence resonance energy transfer (BRET) to monitor receptor-G protein and G protein subunit interactions [48,49] have the potential to quantify these interactions in real time. Modeling studies by our group and others [12,33,[78][79][80][81][82][83] facilitate our understanding of the complexities involved in cell signaling and identify pathway interactions that are key to describing normal and pathological cell functions. The dynamic model that we present here, a general model of GPCR and G protein activation, can be easily expanded to include details of particular GPCR systems and phenomena such as receptor desensitization and internalization, processes that do not operate under steady state conditions. Further, as the components involved in signal transduction become better described and incorporated into signaling databases (e.g., SigPath [84], Alliance for Cell signaling [85], National Center for Genome Resources's PathDB [86]) and as more complex models of signaling pathways become possible, it will become increasingly important to systematically assess parameter uncertainty and quantify regimes of model behavior.

Materials and Methods
The cubic ternary complex activation model. The cubic ternary complex activation model (cTCAM) (Figure 1) is a kinetic extension of the (equilibrium) cTCM integrating a feedback loop of G protein activation and recycling with dynamic LRG interactions [22]. Each reaction in the model is governed by mass action kinetics. Reaction rate constants describe the association and dissociation of ligand and receptor (top to bottom of the cube) and receptor and G protein (front to back of the cube) as seen in Figure 1. The interconversion of inactive and active receptor states is described by forward and backward rate constants as viewed from left to right of the cube.
The cTCAM has a total of 27 kinetic rate constants, 24 describing the binding and dissociation reactions between ligand, receptor, and G protein, and three describing the G protein activation cycle. Additionally, the total concentrations of ligand, receptor, and G protein are required, bringing the total number of model parameters to 30. As described in Text S1, the model can be reduced to 16 input parameters (plus ligand concentration). A description of all the parameters used in the cTCAM is presented in Table 1. Briefly, four of these are ligand-specific parameters that control the binding of ligand to the receptor and the ability of ligand to bias receptor activation and G protein association (K a , k 3 , a, and c, respectively). The eleven cell-specific parameters in the model determine the strength of association of the inactive and active receptor states with G protein (K g , k 11 , g, and b), the ratio of active to inactive receptor (k 1 , K act ), the kinetics of G protein activation, deactivation, and recombination (k Gact , k GTP, and k G ), and the numbers of G proteins and receptors (G total , R total ). The remaining constant, d, is both a ligand-and cell-specific parameter that governs the synergistic effects between ligand binding, receptor activation, and G protein association [11].
The equations describing the cTCAM are presented in Text S1. To calculate the initial conditions for each species in the model under varying parameter values, the total receptor and G protein concentrations (R total and G total ) were set as the initial values of R and G while the initial values of all other species were set to zero. The system was allowed to come to steady state (typically less than 100 s) at which point a bolus of ligand (L) was added to the system (denoted as time ¼ 0), and the formation of each species was tracked over time.
Model analysis. The level of Ga GTP is the key output of the model. However, experimental data are routinely normalized to basal values (i.e., the amount of activated G protein in the absence of ligand). Therefore, for the model G protein activation was quantified as the percent change in Ga GTP from basal (L ¼ 0) as given by Quantifying response in this way is directly related to the pharmacological classifications of ligand efficacy-positive (increase in %OverBasal), neutral (no change in %OverBasal), and inverse agonism (decrease in %OverBasal)-and is analogous to that previously used to analyze response activation in the cTCM and cTCAM models [20,22,81]. Dose response curves were generated by calculating the integral of Ga GTP ( R Ga GTP [t] dt) over the first 10 s of ligand stimulation at varying ligand concentrations. These dose-response curves are used to look at how responses may change upon varying receptor and G protein totals (Figure 4) and for comparison with experiments that measure the accumulation of radioligand, for example [ 3 H]cAMP accumulation in studies of protean agonism in b 2 -adrenergic and a 2Aadrenergic receptors (Figures 6 and 7). Values of the integral were normalized to basal values according to %Accum ¼ 100 To distinguish this measurement from %OverBasal, this analysis is termed ''%Accum.'' Uncertainty and sensitivity analysis using Latin hypercube sampling and partial rank correlation coefficient. To efficiently sample the ranges over which input parameters (x i , i ¼ 1,2,..,X) may vary, LHS was implemented for the cTCAM using methods described previously [29,30,87,88]. Briefly, each input parameter is assigned a range according to values found in the literature (Table 1). Each parameter's value range was divided into N equal probable segments according to a specified probability distribution function for that parameter. For a typical simulation, N ¼ 1,000. Distributions for most of these parameters are not known and therefore uniform distributions were used for each parameter. A random value was then chosen from each segment, so that each parameter became a vector of N values. The N values of the parameter vectors were then randomly paired to generate an N by X input matrix where X was the number of parameters to be varied (X ¼ 16 for the cTCAM). The differential equations describing the model (Text S1) were then solved, generating a vector of N solutions.
PRCCs were calculated to quantify the relative importance of each parameter in generating a desired output (measures of Ga GTP described above). Partial correlation measures the strength of the linear relationship between the output and an input variable (x i ) after the effect of all other elements of x have been removed [38], while the rank transformation is used to linearize the nonlinear monotonic relationship between the input parameters and the output [29]. PRCC values vary between À1 (perfect negative correlation) and 1 (perfect positive correlation). PRCC values were calculated as described previously [29,30,87,88]. Briefly, solutions of Ga GTP at desired time points were added to the LHS input matrix to generate an N by X þ 1 matrix. The values for each of the X þ 1 parameters were then ranked from 1 to N and the resulting matrix was used to calculate a partialized matrix in which the linearized effects of the other parameters are taken out of each parameter. Correlation coefficients were then calculated from the partialized matrix. Scatterplots were generated to assure that the monotonicity assumption applies [29]. The calculated PRCC values were differentiated based on p-values derived from a Student's t test and were then ranked according to their absolute value. Figure S1. The Cubic Ternary Complex Activation Model with Rate Constants Found at doi:10.1371/journal.pcbi.0030006.sg001 (77 KB TIF).