Modelling and simulation of biased agonism dynamics at a G protein-coupled receptor.

Theoretical models of G protein-coupled receptor (GPCR) concentration-response relationships often assume an agonist producing a single functional response via a single active state of the receptor. These models have largely been analysed assuming steady-state conditions. There is now much experimental evidence to suggest that many GPCRs can exist in multiple receptor conformations and elicit numerous functional responses, with ligands having the potential to activate different signalling pathways to varying extents-a concept referred to as biased agonism, functional selectivity or pluri-dimensional efficacy. Moreover, recent experimental results indicate a clear possibility for time-dependent bias, whereby an agonist's bias with respect to different pathways may vary dynamically. Efforts towards understanding the implications of temporal bias by characterising and quantifying ligand effects on multiple pathways will clearly be aided by extending current equilibrium binding and biased activation models to include G protein activation dynamics. Here, we present a new model of time-dependent biased agonism, based on ordinary differential equations for multiple cubic ternary complex activation models with G protein cycle dynamics. This model allows simulation and analysis of multi-pathway activation bias dynamics at a single receptor for the first time, at the level of active G protein (αGTP), towards the analysis of dynamic functional responses. The model is generally applicable to systems with NG G proteins and N* active receptor states. Numerical simulations for NG=N*=2 reveal new insights into the effects of system parameters (including cooperativities, and ligand and receptor concentrations) on bias dynamics, highlighting new phenomena including the dynamic inter-conversion of bias direction. Further, we fit this model to 'wet' experimental data for two competing G proteins (Gi and Gs) that become activated upon stimulation of the adenosine A1 receptor with adenosine derivative compounds. Finally, we show that our model can qualitatively describe the temporal dynamics of this competing G protein activation.

Introduction

Mathematical modelling and scientific computing are powerful tools for the analysis of cell signalling in pharmacology. “Analytical pharmacology”, which has its roots in classical receptor theory and largely focuses on equilibrium cell responses to drugs, provides a vital theoretical basis which underpins drug classification and prediction of drug mechanism of action (Kenakin and Christopoulos, 2011). Much of the analysis has centred on assumptions of a single ligand binding a monomeric G protein-coupled receptor (GPCR), activating a single active state and coupling a single G protein. Concepts like allosterism, inverse agonism, oligomerisation and “biased signalling” are now widely accepted and have enhanced receptor theory towards better understanding of drug-receptor interactions and informed drug discovery (Kenakin and Williams, 2014). GPCRs represent a target for perhaps up to half of all current drugs (Woodroffe et al., 2009), and as such, development of the theory for ligand-GPCR interactions and their consequences is key.

Biased agonism is now a widely accepted phenomenon whereby a ligand may activate multiple different pathways at the same receptor, via multiple active conformations (Kenakin, 2011, Onaran, Rajagopal, Costa, 2014, Rankovic, Brust, Bohn, 2016, Urban, Clarke, Von Zastrow, Nichols, Kobilka, Weinstein, Javitch, Roth, Christopoulos, Sexton, et al., 2007). Other terms for this phenomenon include functional selectivity and pluri-dimensional efficacy, while receptor promiscuity refers to the ability of a receptor to couple different G proteins with different affinities, via different active states. The possibility of multi-pathway activation may lead to a breakdown in the common classifications of ligands based on single active state theory (Kenakin, 2011), or errors in the interpretation of data using simple models (Tuček et al., 2002). Therefore, development of biased agonism theory has become an important field of pharmacological research.

Biased signalling has implications for drug discovery, including the prospect of clinical selectivity and the potential of reduced side effects (Kenakin, 2015, Kenakin, Christopoulos, 2013, Stott, Hall, Holliday, 2015). A schematic of biased agonism is shown in Fig. 1, indicating possibility for a ligand to activate two (or more) G protein pathways at the same receptor, one of which may be a “target” therapeutic pathway, while the other may be an unwanted “side-effect” pathway (panel (b)). To understand, quantify and exploit the potential for biased agonism, theoretical models for such schematics are required.

Fig. 1

Pluri-dimensional efficacy and biased agonism at a GPCR. (a) A classical view of signalling–two different receptors, each bound and activated (to a single active conformation) by a specific ligand, and bound by a specific G protein. The activated G protein subunit α signals to a downstream pathway specific to the G protein. (b) A two-active-state, two-G protein biased signalling schematic. The receptor has two active states, and the proportion of receptors in either active state, and the inactive state, may be affected (biased) by a single ligand. Two different G proteins, specific to the active conformations, couple to the receptors and signal to two pathways. (c) Pluri-dimensional efficacy - multi-active receptor with multiple G proteins, not necessarily each specific to a single receptor conformation. Here we have active states (represented by yellow, green, red and blue in the receptor block) and G proteins. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

A two-active-state model of ligand binding and receptor activation at equilibrium was presented in Leff et al. (1997). This equilibrium model addressed the limitations of single-active-state theory which could not recapitulate different pathway potency and efficacy patterns at the same receptor. It was found that theoretically, an agonist may enrich one active receptor state at the expense of another, and pathway-dependent efficacy was observed in simulations. For an intact system, however, pathway-dependent potency (with active receptor as the pathway readout) was not possible. G protein coupling and activation were not explicitly modelled, but their importance for future modelling was acknowledged. Later equilibrium models included the binding of G proteins (Ehlert, 2008, Scaramellini, Leff, 2002), which give further scope for pathway-dependent pharmacology. An alternative model for biased agonism is given in Roche et al. (2013), where downstream effects are modelled not explicitly via G protein binding, but by coupling the operational model of agonism (Black and Leff, 1983) to active receptor stimuli. This model does not include constitutive activity of the receptors. Further equilibrium modelling for promiscuous coupling of receptors to multiple G proteins has been presented in Kukkonen et al. (2001) and Tuček et al. (2002).

The direction and magnitude of a ligand’s bias towards one pathway over another has largely been quantified using equilibrium assumptions and empirical models such as the operational model (Gundry, Glenn, Alagesan, Rajagopal, 2017, Kenakin, 2014, Kenakin, Watson, Muniz-Medina, Christopoulos, Novick, 2012, Rajagopal, Ahn, Rominger, Gowen MacDonald, Lam, DeWire, Violin, Lefkowitz, 2011). A recent study (Herenbrink et al., 2016) has highlighted the role of “kinetic context” in approaching such calculations, whereby the apparent bias of a ligand towards any given pathway may vary over time. Interpretation of experimental readouts in terms of bias must therefore take into account the signalling dynamics and associated timescales of the measured pathway. Thus, dynamic models of GPCR biased signalling are proposed here to give new theoretical insights into the effects of biased agonists.

In Chen et al. (2003), an ordinary differential equation (ODE) model for the dynamics of biased signalling at GPCRs is presented. The steady-state behaviour of the model is analysed, with particular attention paid to the effect of G protein concentration, where the model output is active G protein. The dynamics in Chen et al. (2003) are not examined in detail, but extensive analysis of GPCR signalling dynamics has been presented elsewhere (Bridge, King, Hill, Owen, 2010, Woodroffe, Bridge, King, Chen, Hill, 2010, Woodroffe, Bridge, King, Hill, 2009), for mathematical models which also allow G proteins binding to inactive receptors, and constitutive receptor activity. In these models, the active G protein α subunit bound to guanosine triphosphate (α) is taken as a model readout which is representative of downstream signalling pathway activity.

In this paper we develop a new mathematical model for the dynamics of biased agonism at GPCRs. The model allows an in-depth theoretical analysis of time-dependent biased agonism at a GPCR for the first time, and is novel in its generality and detail; any number of active receptor states and G proteins may be considered, receptor states need not be specific to particular G proteins, and the response is at the level of α, downstream of active receptor and towards a dynamic functional response. In Section 2, we formulate a general ODE model for the dynamics of a receptor which can activate multiple G protein-mediated pathways. The general model has receptor with N* active conformations and N G proteins available for coupling, but our focus computationally (driven by Leff et al., 1997) throughout is the case . In Section 3, we present time course and concentration-response simulation results for our model, focusing on α dynamics. In particular, we highlight that our model has the propensity for agonist-inverse agonist interconversion both with respect to time and constitutive activity. A numerical analysis of the effects of multiple cooperativity factors is performed. In Section 4, we propose an heuristic method for quantifying dynamic bias, by way of bias factors, and show how these bias factors relate to our model parameters. It is shown that the bias rank order for a bank of ligands may change dynamically. In Section 5, we show that our model simulations fit well to new experimental data where biased agonism at the adenosine A1 receptor is suspected. We conclude in Section 6 with a discussion of our main results, underlining our contribution to the biased signalling literature.

Model formulation

Here we formulate an ODE model for the dynamics of signalling for multi-active state GPCRs capable of binding multiple G proteins, in response to a single ligand binding. The model allows for a receptor which may have an inactive conformation R, or one of N* active receptor conformations R*, for . Also, a receptor may couple to one of N G proteins G, for . The model encompasses ligand binding, receptor activation, G protein binding and the G protein cycle, whereby the model output is activated G protein α, which signals to second messengers, and is therefore taken as an indicator of pathway response, as in Woodroffe, Bridge, King, Chen, Hill, 2010, Woodroffe, Bridge, King, Hill, 2009 and Bridge et al. (2010).

A three-state (two active states) model

While the model is formulated for general N and N*, we largely focus throughout on the case . A schematic for the transitions between 18 receptor states for this particular case is shown in Fig. 2. R denotes inactive receptor, while R* () denotes the jth active state. Any species label including L represents a complex including ligand-bound receptor, while any species including G () is a complex including receptor coupled to the θth G protein. Double arrows represent the reversible binding and activation reactions between receptor states. As described in previous GPCR signalling studies (eg. Bridge, King, Hill, Owen, 2010, Woodroffe, Bridge, King, Chen, Hill, 2010, Woodroffe, Bridge, King, Hill, 2009), a R*G or LR*G complex may dissociate and exchange GDP for GTP on the α subunit of the G protein, leading to the signalling response and the G protein cycle.

Fig. 2

A multi-cubic ternary complex model schematic for biased signalling with two active receptor states (R*1, R*2) and two G proteins (G1, G2), giving 18 receptor species. Double arrows represent reversible binding and activation reactions between the receptor states. The four complexes R*1G1, LR*1G1, R*2G2 and LR*2G2 may dissociate, leading to the G protein cycle and increased active G protein signalling units and .

The (j, θ) receptor/G protein block

In order to formulate the ODE model for the schematic shown in Fig. 2 (or, indeed, the general N*, N-case), we consider the (j, θ) receptor/G protein block (where and for Fig. 2). Each such block is seen to be a cubic ternary complex schema for activation of receptor from inactive state R to active state R*, with coupling to G protein G (Woodroffe et al., 2009). In Fig. 3, the equilibrium rate constants K• and cooperativity factors μ, ν, ζ are labelled on each reversible reaction. For the individual kinetic rate constants and factors, we use lower case k, and subscripts + and to denote the forward and backward reactions respectively. The descriptions of the rate constants and cooperativity factors are given in Table 1. The G protein cycle and α responses follow from dissociation of R*G and LR*G according to the following reactions (see Woodroffe et al., 2009):

Fig. 3

The (j, θ) receptor/G protein block of the multi-cubic ternary complex schematic, for ligand binding to, and activation of, receptor j, with coupling to G protein θ. Equilibrium rate constants K• and cooperativity factors μ, ν, ζ are labelled on each reversible reaction.

Table 1

Equilibrium rate constants and cooperativity factors for the (j, θ) block of the biased signalling schematic.

Label	Description of equilibrium constant
KL	Association of ligand L and receptor R.
KGθ	Binding of G protein Gθ to receptor R.
Kactj	Activation of receptor R to give active state R*j.
μj, θ	Preference of Gθ for R*j over R. Equally, the factor increase in propensity for R → R*j activation when R is Gθ-bound.
νθ	Preference of L for RGθ over R. Equally, the preference of Gθ for LR over R.
ζj	Preference of L for R*j over R. Equally, the factor increase in propensity for R → R*j activation when R is L-bound.

Governing equations

Suppose in general that a receptor has N* distinct active states, and that each receptor may couple one of N distinct G proteins. Then applying mass action kinetics to our schematic and G protein cycle reactions gives a system of n nonlinear ODEs for the species concentrations, where The first term here is given by species L, R and LR. The second term is given by active non-coupled receptor states R*, LR* and the third term corresponds to G protein not coupled to active receptor (G, RG, LRG, βγ). Finally, the number of active receptor/G protein complexes, R*G and LR*G, is 2N*N, since we consider and . If ligand concentration is considered constant, then we will not have an ODE for [L] (so omitting Eq. (3b) below), and instead .

For the model “outputs”, or downstream responses of the system to an input ligand concentration, we take the concentrations for as we consider these as indicators of downstream activity in signalling pathways as in Bridge (2009) and Bridge et al. (2010). For our computational results, we will consider the case with two G proteins and two active receptor states, such that and our model has 18 receptor states and 8 non-receptor-bound G protein species (). Taking ligand concentration constant (as in previous studies), the system (3) in this case consists of 26 ODEs.

Initial conditions for our simulations have (the total receptor concentration), (the total concentration for each G protein), and all other species zero at .

Simulation results

Here we present numerical results (for concentrations) which illustrate the variety of dynamic behaviour which is possible for a system of two active states and two-G proteins. These results are intended to demonstrate potential dynamics rather than provide exhaustive or accurate predictions for any particular receptors or ligands. For all simulations, we first compute the system with for a long time (108 s) to allow the system to come to a steady-state equilibrium before the addition of ligand, and all parameters except those explicitly stated are maintained at the values in Table A.1.

Table A1

Parameter values for 2 G protein, 2 active receptor state model.

Label	Meaning	Cell or ligand specific	Value	Units	Source
k_L+	Ligand binding rate	Ligand	9.40E+04	M^−1 s^−1	Bridge et al. (2010)
k_L−	Ligand unbinding rate	Ligand	3.10E-01	s^−1	”
kact+1	Receptor activation rate to R*1	Cell	1.00E+00	s^−1	”
kact−1	Receptor deactivation rate from R*1	Cell	1.00E+03	s^−1	”
kact+2	Receptor activation rate to R*2	Cell	1.00E+00	s^−1	”
kact−2	Receptor deactivation rate from R*2	Cell	1.00E+03	s^−1	”
kG+1	G protein 1 binding rate	Cell	1.00E+08	M^−1 s^−1	”
kG−1	G protein 1 unbinding rate	Cell	1.00E-01	s^−1	”
kG+2	G protein 2 binding rate	Cell	1.00E+08	M^−1 s^−1	”
kG−2	G protein 2 unbinding rate	Cell	1.00E-01	s^−1	”
kGRA+1	G protein 1 reassociation rate	Cell	7.00E+05	M^−1 s^−1	”
kGRA−1	G protein 1 dissociation rate	Cell	1.30E-03	s^−1	”
kGRA+2	G protein 2 reassociation rate	Cell	7.00E+05	M^−1 s^−1	”
kGRA−2	G protein 2 dissociation rate	Cell	1.30E-03	s^−1	”
khyd+1	Hydrolysis rate of GαGTP1	Cell	1.00E-02	s^−1	”
khyd−1	Exchange rate of GTP for GDP at Gα1	Cell	1.00E-04	s^−1	”
khyd+2	Hydrolysis rate of GαGTP2	Cell	1.00E-02	s^−1	”
khyd−2	Exchange rate of GTP for GDP at Gα2	Cell	1.00E-04	s^−1	”
kGTP+1,1	R*1G1 dissociation rate	Cell	1.00E+00	s^−1	”
kGTP+1,2	R*1G2 dissociation rate	Cell	1.00E+00	s^−1	”
kGTP+2,1	R*2G1 dissociation rate	Cell	1.00E+00	s^−1	”
kGTP+2,2	R*2G2 dissociation rate	Cell	1.00E+00	s^−1	”
ν+θ	Forward cooperativity factor for ligand binding a Gθ bound receptor	Ligand	1.00E+00		”
ν−θ	Backward cooperativity factor for ligand binding a Gθ bound receptor	Ligand	1.00E+00		”
ζ+j	Forward cooperativity factor for ligand-bound Rj activation	Ligand	1.00E+03		”
ζ−j	Backward cooperativity factor for ligand-bound Rj activation	Ligand	1.00E+00		”
μ+j,θ	Forward cooperativity factor for Gθ-bound Rj activation	Cell	1.00E+00 (j=θ), 0 (j ≠ θ)		”
μ−j,θ	Backward cooperativity factor for Gθ-bound Rj activation	Cell	1.00E+00		”
Rtot	Total receptor concentration	Cell	4.15E-10	M	”
Gtot1	Total G1 concentration	Cell	4.15E-10	M	”
Gtot2	Total G2 concentration	Cell	4.15E-10	M	”
Ltot	Total Ligand concentration	Ligand	1.00E-07	M	”

Time courses

We consider time courses for the two α responses of interest. The signifcance of these computations is that two responses are generated from a single ligand at a single receptor, whereas previous GPCR dynamic models (eg. Bridge et al., 2010) have considered a single α response from each receptor. Further, the dynamic responses for the new model do not necessarily follow the previously reported behaviour of ligands classed as agonists, antagonists or inverse agonists for a single active state. An agonist is a ligand which encourages receptor activation, an antagonist is neutral in its action, and an inverse agonist discourages receptor activation. Within our model, therefore, an agonist for active state R* has ζ > 1, an antagonist has and an inverse agonist has ζ < 1.

Ligand is an agonist for both pathways

By varying the values of ζ1 and ζ2 , the preference of the ligand for a receptor in the active states 1 and 2 over the inactive receptor state, we vary the efficacy with respect to the G protein pathways 1 and 2 respectively. In Fig. 4, we show time courses of the responses to a ligand which is an equilibrium agonist for both pathways, for a range of concentrations. Three different ligand concentrations are used, and the responses for are shown. The higher efficacy with respect R*1 gives an increased response, and we note the peak-plateau dynamics. With increased ligand concentration, we see a higher response for both pathways, both at peak and plateau (end-point). Further, the peak timing is reduced with increased ligand concentration, in keeping with previous single active state studies (Bridge, 2009, Woodroffe, Bridge, King, Hill, 2009).

Fig. 4

The response (M) against time (in seconds) of two competing pathways with and after the addition of  M.

Ligand is agonist for one pathway and antagonist for the other

Neutral antagonists may be used as competitive ligands to endogenous agonists. Mathematical modelling of agonist-antagonist competition at a single active state GPCR has been considered in Bridge et al. (2010). Within our two-active state model, we may simulate the dynamics of a system for which a given ligand is an agonist for one pathway but an antagonist for the other. In Fig. 5, we show α and receptor time courses for this scenario, for a ligand which is an (equilibrium) agonist for pathway 1 ( ) and an (equilibrium) antagonist for pathway 2 ( ), over a range of ligand concentrations. We note the peak-plateau dynamics, and the nearly neutral effect on dynamics. However, closer inspection of reveals that the ligand in fact has an inverse agonist effect on pathway 2. Since the ligand is an agonist for pathway 1, its effect on overall receptor activation is an increase in pathway 1 active states, given by and a corresponding decrease in pathway 2 active states and free inactive receptor states, given, respectively, by and Upon ligand addition, there are, therefore, fewer receptors available to activate pathway 2, giving a decrease in and the inverse agonist effect of the “antagonist”. For  M, we also see the undershoot α response previously reported for inverse agonists (Bridge, 2009). We note that a true neutral antagonist effect with constant would be seen if we considered pathway 2 as an “isolated pathway” (see Leff et al., 1997) by setting .

Fig. 5

The response (M) against time (in seconds) for two pathways with and after the addition of  M so that the ligand is an agonist for pathway 1 and an antagonist for pathway 2. Receptor concentrations are given in the bottom panel. Here, .

Ligand is agonist for one pathway and inverse agonist for the other

Having seen apparent inverse agonist activity in the biased system for a ligand which would neutrally antagonise an isolated pathway, we now turn attention to a ligand which is a true inverse agonist for one pathway in the biased system, and an agonist for the other. In Fig. 6, we show α time courses for this scenario, for a ligand which is an (equilibrium) agonist for pathway 1 ( ) and an (equilibrium) inverse agonist for pathway 2 ( ), over a range of ligand concentrations. These simulations are for a system with increased R*2 constitutive activity, to represent conditions under which inverse agonism may be detectable. We note the peak-plateau dynamics, and drop-off in level. Further, we observe undershoot dynamics in the inversely agonised pathway, which may be seen in a single-active state system (Bridge, 2009). An interesting feature here is that while increasing ligand concentration decreases peak time as before, this is accompanied by an increase in trough time.

Fig. 6

The response (M) against time (s) after the addition of a range of ligand concentrations, where the ligand is an agonist for pathway 1 and an inverse agonist for pathway 2. Here, .

Time course surfaces

In order to summarise the effect that efficacy parameter ζ has on a system, in Fig. 7 we show time course surfaces for  M, where we vary over a spectrum of efficacy ranging from strong inverse agonist to strong agonist, while keeping all other parameters fixed. We clearly see that the stronger L is an agonist for pathway 1, the lesser its effect on pathway 2. When L is an agonist for both pathways, the peak-plateau dynamic response is clear for both and but increasing agonist strength for pathway 1, the pathway 2 response drops off. Similarly, for a pathway 2 antagonist, the dynamic response varies from apparent antagonism to inverse agonism with increasing . Also, for a pathway 2 inverse agonist, the magnitude of inverse agonism increases with .

Fig. 7

Time course surfaces for response (M) dynamically changing for system with varying agonist efficacy parameter acting under a ligand concentration of  M. Column (a): ; column (b): ; column (c): .

Observed agonist effect is system-dependent - constitutive activity and inter-conversion

A feature of the equilibrium three-state model in Leff et al. (1997) is that a ligand’s effect on a pathway can change qualitatively from agonist to inverse agonist, depending on the system-specific level of constitutive activity in that pathway. This so-called “inter-conversion” of ligand effect is demonstrated at steady-state in Leff et al. (1997), with respect to active receptor states. In Fig. 8, we show the effects of increasing the constitutive activity in pathway 2 by decreasing . With low constitutive activity (), the time courses for are indistinguishable. As pathway 2 constitutive activity is increased, it is clear that pathway 2 basal α increases at the expense of pathway 1 basal α, similarly to the active receptor trend in Leff et al. (1997) and Scaramellini and Leff (1998). The agonist effect on becomes less pronounced with decreased as the G protein response is largely effected via the basal activity, but the ligand remains a pathway 2 agonist. In contrast, with high activation of pathway 2 (), the long-time response decreases with “agonist” concentration, so that the ligand is now having an apparent inverse agonist effect, despite its isolated pathway classification as an agonist. It is also worth noting that for (), we observe non-monotonicity in the peak and plateau as functions of [L]. Thus non-monotonic concentration-response curves may result from multi-active state receptors with varying constitutive activity levels.

Fig. 8

Time courses for response (M) dynamically changing for systems with differing constitutive receptor activation level in pathway 2, varying . Here, .

Observed agonist effect may be time-dependent: G protein cycle and dynamic inter-conversion

With our new model, we are able to examine the dynamics under variation of constitutive receptor activation. In the final plot of Fig. 8, we see the phenomenon of dynamic inter-conversion between agonist and inverse agonist action. The ligand is an agonist for both pathways under equilibrium classification, but after initially displaying a typical agonist response, eventually drops below basal levels in an apparent inverse agonist response. The dynamic peak response to agonism occurs as in previous simulations (Woodroffe et al., 2009). The below basal long-time level is a result of G protein cycle dynamics on active receptor equilibration. As is inactivated and G1 reassociates, any new free receptors resulting from LRG complex dissociation are pulled towards a pathway 2 dominant equilibrium, and the receptor pool for G1 activation decreases below basal level.

Concentration-response relationships

Peak and plateau responses with varying ligand activation efficacy and constitutive receptor activity

The ligand concentration-dependent features which summarise the α equilibrium and dynamic behaviour may be summarised using conventional concentration-response curves. In Fig. 9, we show concentration response curves for both pathways, where the measured responses are the peak and plateau α levels. The non-monotonicity first noted in Section 3.1.5 is clearly a possibility. For a ligand which agonises both pathways, with high constitutive activity in one pathway, the plateau response in the other pathway is non-monotonic. The peak concentration-response curve is yet more complex; it is also non-monotonic, with a biphasic structure. We remark that biased agonism together with constitutive activity in our new model for α response is a mechanism by which non-monotonic concentration-response relationships can occur. Such behaviour cannot be observed for three-state models (Leff, Scaramellini, Law, McKechnie, 1997, Scaramellini, Leff, 1998) where the “readout” is a particular active receptor fraction.

Fig. 9

Concentration-response curves ( concentration against ligand concentration) where the ligand is an agonist for both pathways ( top row) and agonist for pathway 1 but an inverse agonist for pathway 2 ( bottom row). Constitutive activity for pathway 2 is low ( left column), medium ( middle column), and high ( left column). Here, .

Further demonstration of the dynamic and concentration-dependent features of the system is given in Fig. 10, where we observe decreasing peak timing for both pathways where the ligand is an agonist for both, but an increasing trough time at the pathway for which the ligand is an inverse agonist.

Fig. 10

Concentration-response curves ( maximum and minimum timing against ligand concentration) where the ligand is an agonist for both pathways ( top row) and agonist for pathway 1 but an inverse agonist for pathway 2 ( bottom row). Constitutive activity for pathway 2 is low ( left column), medium ( middle column), and high ( left column). Here, .

Effect of total receptor number on concentration-response

The total concentration of receptor can considerably affect the appearance of bias in a system (Rajagopal et al., 2011). In Fig. 11, we investigate the effect of differing receptor expression by examining concentration-response curves for two pathways being agonised by a ligand (L1) with different efficacies ( and at a range of receptor concentrations R (from 4.15  M to  M). As R is decreased, we observe both a rightward shift of the curves (increased EC50), together with a drop in the maximal responses, for both the peak and plateau values of α.

Fig. 11

Concentration-response curves for peak and plateau α for two pathways being agonised by L1, a ligand with different efficacies for the two pathways ( and ), under varying receptor concentrations.

Whilst overall efficacy depends partly on the preference of the ligand for an active rather than inactive receptor (controlled through variation of the ζ parameters), it is important to note that it can also depend on the preference of a ligand-bound receptor for each of the G proteins (mediated by the ν parameters). In the case of a system in which the ligand-dependent parameters affecting efficacy are chosen so as to counteract each other ( ), we see that (Fig. 12) not only the magnitude of the preference for one pathway, but even the direction (in terms of which pathway experiences the higher response) can be affected by a changing receptor concentration. In this case, peak α exhibits a change in direction, while plateau level does not.

Fig. 12

Concentration-response curves for two pathways being agonised by L2, a ligand with parameters vary under changing receptor concentrations.

Response surfaces

Parameter sensitivity and concentration-response relations may be conveniently summarised using response surfaces which show the effects of varying two system parameters (Bornheimer, Maurya, Farquhar, Subramaniam, 2004, Bridge, King, Hill, Owen, 2010, Woodroffe, Bridge, King, Hill, 2009). We now use this method to show the sensitivity of simulated response ( peak and plateau) to variations in system parameters, in particular the microaffinity coefficients ζ, ν and μ.

Effect of ζ - the possibility of biphasic relationships

In Fig. 13, we see the effect of varying and for a fixed ligand concentration. The reciprocal effects on the two G protein pathways mediated by the competing receptor states are clear. As is increased, both peak and plateau increase, accompanied by decreases in . Furthermore, it is clear that biphasic relationships are possible.

Fig. 13

Response surfaces for varying ligand efficacy, for fixed ligand concentration  M. Parameters and are varied through the spectrum of efficacy for each receptor active state.

Effect of G protein non-specificity and receptor “cross-states”

Thus far in our computations, we have focussed on systems in which the two G proteins are each specific to a particular active receptor conformation. By setting we have simulated systems whereby G protein 1 can neither activate a pre-coupled receptor towards R*2, nor bind to R*2, and vice-versa. It is a novel aspect that our model allows receptor “cross-states”, where there is not exclusive specificity of each G protein for one particular active receptor state. We see in Fig. 14 the effects of non-exclusive specificity and accessibility of these cross states on the α responses. The general trend is that with all cross states signalling (with ), increasing μ21 gives a decreased peak and slight increase in plateau .

Fig. 14

Response surfaces for varying ligand accessibility of receptor cross states for fixed ligand concentration  M.

Our model allows for variation in specificity of not only the G proteins for each receptor conformation, but also the propensity for G protein cycling with respect to these active states, controlled by . With cross states which do not signal (with ), increasing μ12 now gives a decreased peak and plateau as the general trend, with non-monotonicity, which may be explained by considering the effects on basal conditions. Further explanation and discussion of these effects is given in Appendix B.

Effect of ν, the preference of ligand for specific G protein-coupled receptor

The microaffinity constant ν controls the preference of ligand for RG over R. The effect of varying ν is as expected, in that increasing ν increases both peak and plateau (see Fig. C.1 in Appendix C).

Fig. C1

Response surfaces for varying with expected monotonic relationships. Here,  M.

Detecting and quantifying bias

A balanced agonist is one which signals with equal efficacy to available downstream pathways, whereas a biased agonist has different efficacies for signalling to different pathways (Rajagopal et al., 2011). There is a need to detect and quantify the level of bias towards one pathway over another, considering that physiologically and clinically, certain pathways represent therapeutic targets while others are “side effect” pathways (Gundry et al., 2017). Here, we employ current quantification methods for the level of ligand bias within our two-pathway system.

Bias factors and the operational model of agonism

The operational model of agonism (Black and Leff, 1983) provides a standardised and widely adopted method for estimating ligand affinity and “operational efficacy” parameters from functional response data in the form of hyperbolic concentration-response curves. Briefly, for a single downstream readout E resulting from ligand concentration [A] at a receptor, where E is the maximum response of the system, K is the ligand’s equilibrium dissociation constant, and τ is a measure of ligand efficacy, in particular measuring the propensity of the ligand and the system to yield a response. A modified form of the model is sometimes used to account for nonzero basal responses (Shonberg et al., 2014), namely Further generalisation of this model is possible by introducing a Hill coefficient to the signal transduction sub-model (Black, Leff, 1983, Shonberg, Lopez, Scammells, Christopoulos, Capuano, Lane, 2014). Recently, the operational model has been used to quantify the level of bias in systems exhibiting multi-dimensional efficacy (ie. the activation of multiple pathways at a single receptor). Bias is typically defined with respect to a reference ligand, and bias factors are computed using fitted values of τ (Rajagopal et al., 2011) or both τ and K (Gundry, Glenn, Alagesan, Rajagopal, 2017, Kenakin, 2014, Kenakin, Christopoulos, 2013). Here, we follow the transduction coefficients method (Kenakin and Christopoulos, 2013) by defining a transduction coefficient for a ligand A at a given pathway as The difference in transduction coefficients for two ligands A and B is usually written in Δlog  notation, with The relative bias factor for a ligand A, relative to ligand B, for pathway 1 over pathway 2, is usually defined by first calculating its logarithm, written in ΔΔlog  notation as so that

Bias factor’s dependence on ζ and ν

The bias factor, is a standard measure of a ligand’s bias for effecting a response in pathway 1 over pathway 2. While this is defined in terms of the parameters τ and K which are fitted to the semi-mechanistic operational model, rather than explicitly in terms of the parameters in our new α model, we expect correlations between the bias factor and ligand-specific parameters in our model. In particular, when α is measured at equilibrium and taken as the response E, we expect, on the whole, bias factor should increase with decreased ζ2 and ν2, which control a ligand’s effect on R*2 activation and G2 coupling to the receptor. In Fig. 15, we show the bias factor for a bank of ligands generated by varying and while keeping all other parameters fixed. The correlation is clear. The overall trend is the expected increase in with decreasing and while the relationship is approximately a power law over much of the parameter space shown.

Fig. 15

Bias factor surface for fixed ligand 1 parameters, varying ligand 2 parameters ζ2 and ν2. Bias factor represents the bias for pathway 1 over pathway 2 signalling. Here, and .

Kinetic context and dynamic bias factors

It has recently been demonstrated that binding, activation and signalling dynamics may significantly affect bias measurements, and hence the classification of biased ligands, and that “kinetic context” is an important consideration in the quantification of bias (Herenbrink et al., 2016). Although bias calculations based on the operational model implicitly assume equilibrium conditions, this method is shown to be an effective and simple heuristic approach to investigating and quantifying dynamic bias in Herenbrink et al. (2016). In Fig. 16, we show bias factor time courses for a bank of ligands, generated by constructing a concentration-response curve for each ligand at each time point, then fitting each of these curves to the operational model (using optimisation routines in MATLAB). Here, the bias factor is calculated with respect to the reference ligand (ligand 7), and we see that the long-time bias factor indeed increases with decreased and/or . We also plot, in the right hand panel, an alternative bias factor based on our model parameters, specifically which should also indicate ligand bias. We note the excellent agreement between the dynamic bias factors from operational model fitting and our alternative bias factor. Dynamically, there is the possibility of a change of order of bias factors, and this phenomenon is even more marked for the bank of ligands shown in Fig. 17. Clearly, the order of bias factors may change dynamically, so the classification of ligands requires consideration of kinetic context, as described in Herenbrink et al. (2016).

Fig. 16

Bias factor dynamics for a bank of nine ligands. The operational model bias factor represents the bias for pathway 1 over pathway 2 signalling is shown for each time point, and the alternative bias factor is shown in the right hand panel. Reference ligand is ligand 7, a strong agonist for pathway 1. Here, and .

Fig. 17

Bias factor dynamics for a bank of four ligands. Bias factor represents the bias for pathway 1 over pathway 2 signalling. Here, and .

Bias dynamics beyond the operational model of agonism

The operational model of agonism is a semi-empirical model used to provide summary measures of signalling efficacy, by implicitly assuming equilibrium conditions; in itself, it cannot be used to simulate biased signalling dynamics. Analysis of biased signalling using a dynamic model which does not rely on the operational model is therefore desirable (for example for fitting to time course data). Our new mechanistic model clearly fulfills this requirement, and will be further applied (beyond our scope here) to define new dynamic bias factors which do not use the operational model assumptions at all.

Fitting to a model of downstream functional antagonism via biased signalling

Our model outputs thus far have been the α levels of the two G proteins in the system, which represent responses downstream of ligand binding, and may correspond to a downstream functional experimental readout. We now consider whether our model can be used to explain, and fit to, experimental end point data in a system where biased agonism is suspected.

Experimental method

The adenosine A1 receptor (A1R) is well known for mediating the protective effects of adenosine in the heart (Donato, Gelpi, 2003, Minamino, 2012). How these effects are brought about is not fully understood, as the A1R is able to couple to multiple signalling pathways (Baltos et al., 2016). This makes interpretation of physiological effects difficult to attribute to an individual, signalling pathway. While the A1R is a predominantly G-coupled receptor, which inhibits the accumulation of the second messenger, cyclic adenosine monophosphate (cAMP), it has been observed that at higher agonist concentrations, the levels of cAMP begin to rise again producing a non-monotonic response profile. This accumulation of cAMP arises through the ability of the A1R to switch its G protein coupling and now promote activation of G (Baker, Hill, 2007, Cordeaux, IJzerman, Hill, 2004). The extent to which an individual agonist either inhibits or stimulates cAMP production at the A1R may vary.

To obtain data to enable fitting of our models, experiments were performed using Chinese hamster ovary-K1 (CHO-K1) cells stably expressing the A1 receptor (these cells do not endogenously express any of the four adenosine receptor subtypes, and therefore provide a null background (Knight et al., 2016)), treated with a range of concentrations of a single agonist each time, and the effect on intracellular concentration of cAMP determined (see Appendix D and Knight, Hemmings, Winfield, Leuenberger, Frattini, Frenguelli, Dowell, Lochner, Ladds, 2016, Weston, Lu, Li, Barkan, Richards, Roberts, Skerry, Poyner, Pardamwar, Reynolds, et al., 2015, Weston, Winfield, Harris, Hodgson, Shah, Dowell, Mobarec, Woodlock, Reynolds, Poyner, et al., 2016 for details). In particular, the experiments were carried out for three different agonists individually, namely5’(N-ethyl carboxamido) adenosine (NECA), and two test compounds which bind the A1R, denoted here as Compound 6 (Cmpd6) and Compound 20 (Cmpd20) (Knight et al., 2016). The measured response was the accumulated cAMP concentration in the presence of the phosphodiesterase (PDE) inhibitor rolipram, which blocks cAMP degradation.

For each concentration of all agonists, two experiments were performed. Firstly, intact (wild-type) cells were used, which allow for coupling and activation of both G and G proteins to the A1R, thereby allowing activation of both the G pathway which inhibits cAMP production via an increased α signal, and the G pathway which stimulates cAMP production via an increased α signal. For these cells, the recorded response is the percentage inhibition of cAMP when compared with cells treated with forskolin (which promotes maximal stimulation of cAMP production Barritt, 1992). The second experimental condition is for cells that have been treated with pertussis toxin (PTX), which both inhibits binding of G to its receptor and blocks its signal transduction, thereby locking α in its inactive, GDP-bound state (Mangmool and Kurose, 2011). For these cells, the recorded response is the percentage stimulation of cAMP when compared with forskolin-stimulated cells. Time courses of cAMP were not recorded, and the signalling readout in each case is taken at the endpoint of the experiment ( s).

Further details of the experiments are given in Appendix D.

Experimental results

In Fig. 18, we show the experimental cAMP endpoint signals in response to three ligands individually in turn (NECA, Cmpd6 and Cmpd20) for the two different experimental conditions. For wild-type cells the log concentration response curves for the inhibition of cAMP show non-monotonic behaviour with a downturn at higher concentrations, whereas the log concentration response curves for the production of cAMP in PTX-treated cells show, with the exception of one data point for the NECA experiment, monotonic behaviour. By blocking the inhibitory pathway, we largely see “standard” monotonic behaviour, which suggests that the non-monotonic wild-type response results from crosstalk between the inhibitory and stimulatory pathways. Since in each case a single ligand has been introduced and A1R is the only receptor in the cells, we hypothesise that this target receptor may exhibit biased agonist effects, via two active conformations, one of which is specific to the G protein and the other to the G protein.

Fig. 18

Using the model of biased agonism with functional antagonism at level of α to fit cAMP readouts (signalstim(1800) and signalinhib(1800)) for three different ligands (NECA, Cmpd6 and Cmpd20). Experimental data points are given (red squares with dashed lines for percentage cAMP inhibition for wild-type cells, blue circles with dashed lines for percentage cAMP production in PTX-treated cells, each compared with forskolin-treated cells which give maximal cAMP response) for end point readouts over a range of ligand concentrations. Solid curves are the fitted log concentration response curves for our model. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Modelling considerations

Since the data shown in Fig. 18 are hypothesised to result from biased agonism with competition between two activated G protein pathways, we now seek to fit our model to the data, in order to add support to this hypothesis and understand the possible underlying mechanisms.

Within our modelling framework, we let G1 and G2 represent the G and G proteins respectively. We simulate the PTX effect of blocking G binding and activation by setting . Since cAMP is produced in response to G activation (Barritt, 1992, Leander, Friedman, 2014), for a simple, minimal model of cAMP levels in PTX-treated cells, with blocked cAMP degradation, we take the cAMP production rate proportional to α levels, so that and hence the stimulation signal is given by where C is a constant.

For the wild-type cells in which both stimulatory and inhibitory cAMP pathways are intact, we require a model for crosstalk between G and G pathways. Here we use a simple “functional antagonism” model for the competing effects of these pathways. Functional antagonism refers to the response of a cell in which signalling via one pathway is antagonised by signalling via another pathway, and simple theoretical models have been presented which are based on differences between pathway signals (Leff, Martin, Morse, 1985, Mackay, 1981, Szabadi, 1977). Here, we use such a model where α and α are stimuli to the cAMP stimulatory and inhibitory pathways respectively, and the cAMP production rate is simply a scaled difference between the two α levels. The inhibition signal is then given by where C is a constant, and C is the constant as in (14). Since functionally opposite signalling can result in non-monotonic concentration-response curves with downturns (PliSka, 1994, Szabadi, 1977) such as those seen in the cAMP inhibition curves in Fig. 18, our biased agonism model augmented by the simple functional readout models (14) and (15) may be able to recapitulate the experimental data, at least qualitatively. We proceed to employ parameter estimation methods to pursue a fit to the concentration-response curves for each ligand.

Parameter estimation

We fit the experimental data to the model given by (3a), (3b), (3c), (3d), (3e), (3f), (3g), (3h), (3i), (3j), (3k), (3l), (3m) with together with (14) and (15), where simulations are first run to a time of 108 s with to pre-equilibrate the system before ligand addition. For each ligand, the experimental data for the intact and PTX cells were fitted simultaneously, using optimisation algorithms to minimise the squared error between simulation and data points. The methods used were the trust region algorithm implemented in PottersWheel (Maiwald and Timmer, 2008), followed by a genetic algorithm routine implemented in MATLAB (Matlab, 2017). A subset of the kinetic parameters were varied; for each reversible reaction, we fixed one rate constant (typically for the reverse reaction), and allowed one rate constant to float. Further, we consider systems where the active receptor cross states are inaccessible, such that are fixed for j ≠ θ, since these have been shown to largely have little effect. Fitted parameters for the NECA data set were used as initial parameter guesses for Cmpd6 and Cmpd20, to speed up the overall fitting for these compounds. For each ligand, the model can clearly fit the data very well qualitatively. In Fig. 18, we see that the model fit for the stimulation curve is monotonic, with maximal and basal signals, and EC50 values in good agreement with the data. Further, the fitted inhibition curve in each case is non-monotonic, with the concentration which gives the peak value in good agreement with the data. The basal, peak and plateau levels are in good agreement with the data, and the model recapitulates the differences in peak “spread” between the three ligands. Values for the fitted parameters are given in Table E.1.

Table E1

Parameter estimates for fitting biased agonism with functional antagonism model to cAMP data (see Fig. 18).

Parameter	Cmpd6	Cmpd20	NECA
k_L+	1.20E+04	7.49E+03	1.52E+04
ζ+1	4.47E+03	1.52E+02	4.55E+03
ζ+2	2.70E+04	7.00E+03	2.54E+04
ν+1	5.82E-01	1.22E+01	5.82E-01
ν+2	6.00E+00	2.54E+02	8.36E+00
kact+1	8.80E-01	1.44E+00	9.32E-01
kact+2	6.47E-01	1.40E+00	5.94E-01
kG+1	6.34E+08	5.69E+08	6.55E+08
kG+2	2.89E+09	3.72E+09	2.36E+09
kGRA+1	7.55E+05	7.99E+05	7.63E+05
kGRA+2	4.10E+05	3.50E+05	3.28E+05
khyd+1	2.42E-03	4.40E-03	2.42E-03
khyd+2	3.20E-03	3.00E-03	3.30E-03
kGTP+1,1	1.90E-01	2.01E-01	1.75E-01
kGTP+2,2	1.87E-01	2.27E-01	1.87E-01
μ+1,1	1.28E+00	2.16E+00	1.28E+00
μ+2,2	8.60E-01	8.50E-01	8.60E-01
Rtot	4.15E-10	4.15E-10	4.15E-10
Gtot1	3.79E-10	4.80E-10	4.49E-10
Gtot2	8.22E-10	7.20E-10	8.28E-10
Cs	4.30E+08	4.36E+08	3.24E+08
Ci	4.40E+08	5.30E+08	4.99E+08

In Fig. 19, simulations from the fitted parameter sets for each ligand show the ∫α contributions to the overall measured signal for both cell types. In each case, the stimulatory responses for the intact cells and the PTX cells are almost indistinguishable, and while the individual ∫α curves are monotonic, the difference between them for the intact cells is not.

Fig. 19

Log concentration response curves for and levels using fitted parameter values.

Having estimated parameters which fit the experimental data (taken at a single time point  s), we may now simulate the underlying α dynamics up to this time point. In Fig. 20, we show time courses for α and α levels, using the NECA-fitted parameters. With the ligand being an agonist for both G protein pathways, the peak-plateau α dynamics are clear, and consistent with the temporal characteristics observed in our earlier numerical simulations. Peak values are monotonic with [L], with peaking later than . We conclude that our model recapitulates, and fits to, experimental data well in the cases shown, therefore adding support to the biased agonism conjecture for the experiments discussed, and validating our model. Our functional model for cAMP production is very simple, comprising a linear combination of α and α. It is reasonable to expect that a more detailed model of cAMP signalling with more degrees of freedom would result in an even better fit to the data.

Fig. 20

Underlying α and α dynamics for NECA-fitted parameters, for a range of agonist concentrations. Time is on a logarithmic scale to clearly show the peak-plateau time scales.

Discussion

Biased agonism is now a widely accepted phenomenon for signalling via GPCRs (Kenakin, Christopoulos, 2013, Urban, Clarke, Von Zastrow, Nichols, Kobilka, Weinstein, Javitch, Roth, Christopoulos, Sexton, et al., 2007), and exploiting this is a potential route to developing novel therapeutics (Kenakin, 2015, Kenakin, Christopoulos, 2013). Theoretical (mathematical) models are key tools towards understanding biased signalling, and have previously been presented for equilibrium conditions (Leff, Scaramellini, Law, McKechnie, 1997, Scaramellini, Leff, 2002). These models have enabled a foundation biased agonism theory to be established, largely at the level of receptor activation. For functional readouts downstream of the receptor, further detail has been added at the level of G protein binding  (Ehlert, 2008), and simple empirical models for pathway signalling (Roche et al., 2013). In this paper, we have developed a new model for biased agonism which includes the detail of G protein activation via a cubic ternary complex/G protein cycle model, with α as a readout, and serving as an indicator/proxy of pathway activity. This model is general in terms of the number of active receptor conformations and G proteins, and also that it is not specific to any particular pathway; it can thus be used to model biased signalling at any GPCR, and detailed further signalling components can be added downstream of α to model particular pathways as desired. Potentially, our model could provide a foundation for simulating single-ligand multi-pathway dynamics, such as recent experimental work revealing dynamic biased signalling behaviour at the dopamine D2 receptor (Herenbrink et al., 2016).

An important advance in the present study is the analysis of signalling dynamics as predicted by our model. The role of kinetic context in the investigation of biased agonism has recently been highlighted (Herenbrink et al., 2016) and, as such, a model and method for analysing dynamics represents a timely contribution to the literature. A number of dynamic features have been observed here, including the apparent inverse agonist effect of an “antagonist”, dynamic inter-conversion of agonist effect, and the time-dependence of bias factor order. Non-monotonic concentration-response relationships for endpoint signals are possible from our model, for both a single α readout within a two-pathway system, and downstream crosstalk between two α signals.

The current standard method for quantifying bias from experimental data uses parameters fitted to the equilibrium operational model of agonism (Black, Leff, 1983, Kenakin, Watson, Muniz-Medina, Christopoulos, Novick, 2012). Calculating bias factors using this empirical model applied to timecourse data shows the dynamic nature of bias (Herenbrink et al., 2016), and our model and computations have reproduced this phenomenon. We propose that such analysis may provide important new insights into, and quantitative characterisation of, experimental timecourse results. The use of the operational model appears to be the current state-of-the-art in bias quantification, but it has a number of limitations: It is empirical rather than mechanistic, it does not consider dynamics, and it does not account for constitutive activity (Stott et al., 2015). An alternative model which includes constitutive activity is given in Slack and Hall (2012), but this equilibrium model has yet to be fully explored with respect to biased signalling. While beyond the scope of our current work, a valuable future investigation will focus on further formulation and definition of dynamic bias factors, including constitutive activity.

We have shown our model to be capable of reproducing endpoint trends in experimental data for cAMP levels in response to ligands at the A1R receptor, through multi-pathway (α and α) signalling with functionally opposite downstream signals. This endpoint analysis has resulted in parameterisations of the model which then predict the underlying α dynamics, qualitatively consistent with our earlier agonist-induced simulations. This validation of our model allows us to propose its use for further study of downstream signalling, and fitting to time-course data when it becomes available. For example, for any future dynamic cAMP experimental readouts, our simple functional models (14) and (15) can be used to fit to time-courses, with better fits expected by letting a greater number of parameters float, or using a more detailed cAMP model (eg. Leander and Friedman, 2014). The simulation and fitting in the current work also clearly shows that single-ligand multi-pathway activation at a single receptor provides a mechanism for non-monotonic concentration-response relations either for α itself or for downstream signals, by way of functional antagonism. While functional signalling experiments often results in monotonic concentration-response curves, relationships with downturns at high concentrations are not uncommon (Calabrese, Baldwin, 2001, PliSka, 1994, Zhu, Liu, Qin, Chen, Liu, 2013), and the current work provides a plausible mechanistic model for understanding such results in systems where multi-pathway signalling via a single receptor is possible.

The mathematical work here represents a theoretical framework for further study of the potential benefits of developing biased agonists as therapeutics. The multidimensionality of GPCR signalling now constitutes a new paradigm in drug discovery, and the potential benefits of new understanding of multi-pathway signalling lie in the development of “functionally selective” drugs which preserve efficacy in target pathways, while minimising activation of unwanted side-effect pathways at the same receptor (Rankovic et al., 2016). Further mechanistic modelling encompassing G protein binding and activation, downstream signalling, dynamics and complexity of the level we have studied here is acknowledged as a potentially very valuable advance towards such drug discovery goals (Stott, Hall, Holliday, 2015, Urban, Clarke, Von Zastrow, Nichols, Kobilka, Weinstein, Javitch, Roth, Christopoulos, Sexton, et al., 2007).

Author contributions

LJB developed models, performed computations and analysis, and wrote the manuscript. JM performed computations and analysis, and wrote the manuscript. EF performed wet-lab experiments. IW performed wet-lab experiments. GL developed models, analysed data and wrote the manuscript.

Table F1

Parameter values for 2 G protein, 2 active receptor state model, now measuring concentrations in nM.

Label	Meaning	Value (typical simulation)	Value (fitted for NECA)	Units
k_L+	Ligand binding rate	9.40E-05	1.52E-05	nM^−1 s^−1
kG+1	G protein 1 binding rate	1.00E-01	6.55E-01	nM^−1 s^−1
kG+2	G protein 2 binding rate	1.00E-01	2.36E+00	nM^−1 s^−1
kGRA+1	G protein 1 reassociation rate	7.00E-04	7.63E-04	nM^−1 s^−1
kGRA+2	G protein 2 reassociation rate	7.00E-04	3.28E-04	nM^−1 s^−1
Rtot	Total receptor concentration	4.15E-1	4.15E-01	nM
Gtot1	Total G1 concentration	4.15E-1	4.49E-01	nM
Gtot2	Total G2 concentration	4.15E-1	8.28E-01	nM
Ltot	Total ligand concentration	1.00E+02	–	nM