Kinetic operational models of agonism for G-protein-coupled receptors

Kinetic operational models of agonism for G-protein-coupled receptors

Accepted Manuscript Kinetic Operational Models of Agonism for G-Protein-Coupled Receptors Samuel R.J. Hoare , Nicolas Pierre , Arturo Gonzalez Moya ,...

2MB Sizes 0 Downloads 54 Views

Accepted Manuscript

Kinetic Operational Models of Agonism for G-Protein-Coupled Receptors Samuel R.J. Hoare , Nicolas Pierre , Arturo Gonzalez Moya , Brad Larson PII: DOI: Reference:

S0022-5193(18)30066-3 10.1016/j.jtbi.2018.02.014 YJTBI 9357

To appear in:

Journal of Theoretical Biology

Received date: Revised date: Accepted date:

24 July 2017 7 February 2018 13 February 2018

Please cite this article as: Samuel R.J. Hoare , Nicolas Pierre , Arturo Gonzalez Moya , Brad Larson , Kinetic Operational Models of Agonism for G-Protein-Coupled Receptors, Journal of Theoretical Biology (2018), doi: 10.1016/j.jtbi.2018.02.014

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT Highlights For GPCRs, kinetics is of increasing importance for developing new therapeutics



Here a new model is described for quantifying kinetic pharmacology



Efficacy is defined by a single parameter that is simple to measure



The model enables binding kinetics to be measured using functional data



Templates are provided for application of the equations to experimental data

AC

CE

PT

ED

M

AN US

CR IP T



1

ACCEPTED MANUSCRIPT Kinetic Operational Models of Agonism for G-Protein-Coupled Receptors Samuel R.J. Hoare a, Nicolas Pierre b, Arturo Gonzalez Moya b, and Brad Larson c

Pharmechanics, LLC, 14 Sunnyside Drive South, Owego NY 13827, USA

b

Cisbio Bioassays, 135 South Road, Bedford, MA 01730, USA

c

BioTek Instruments, Inc, 100 Tigan Street, Winooski, VT 05404, USA

CR IP T

a

Corresponding author: Samuel R.J. Hoare

AN US

Pharmechanics, LLC 14 Sunnyside Drive South

M

Owego

NY 13827

ED

USA

AC

CE

Tel: US (619) 203-2886

PT

Email: [email protected]

2

ACCEPTED MANUSCRIPT Abstract The application of kinetics to research and therapeutic development of G-protein-coupled receptors has become increasingly valuable. Pharmacological models provide the foundation of pharmacology, providing concepts and measurable parameters such as efficacy and potency that have underlain decades of successful drug discovery. Currently there are few pharmacological models that incorporate kinetic activity in such a way as to yield experimentally-accessible drug parameters. In this study, a kinetic model of pharmacological response

CR IP T

was developed that provides a kinetic descriptor of efficacy (the transduction rate constant, kτ) and allows measurement of the kinetics of receptor-ligand binding from functional data. The model assumes: 1) receptor interacts with a precursor of the response (“Transduction potential”) and converts it to the response. 2) The response can decay. Familiar response vs time plots emerge, depending on whether transduction potential is

AN US

depleted and/or response decays. These are the straight line, the “association” exponential curve, and the riseand-fall curve. Convenient, familiar methods are described for measuring the model parameters and files are provided for the curve-fitting program Prism 7.0 (GraphPad Software) that can be used as a guide. The efficacy

M

parameter kτ is straightforward to measure and accounts for receptor reserve; all that is required is measurement of response over time at a maximally-stimulating concentration of agonist. The modular nature of the model

ED

framework allows it to be extended. Here this is done to incorporate slow agonist-receptor equilibration and antagonist-receptor binding. In principle, the modular framework can incorporate other cellular processes, such

PT

as receptor desensitization. The kinetic response model described here can be applied to measure kinetic

CE

pharmacological parameters than can be used to advance the understanding of GPCR pharmacology and optimize new and improved therapeutics.

AC

Keywords: drug discovery; drug kinetics; ligand binding; receptor binding; kinetics; receptor theory 1. Introduction

The superfamily of G-protein-coupled receptors is the largest class of receptors in the human genome

(Fredriksson et al., 2003). The superfamily comprises receptors for a large diversity of activating agents, including photons, small organic compounds, peptides and glycoproteins. GPCRs are the targets of a substantial fraction of drugs in clinical use (Rask-Andersen et al., 2014; Sriram and Insel, 2018). Pharmacological theory 3

ACCEPTED MANUSCRIPT has been a framework underlying successful drug development for decades. Classical receptor theory provides the pharmacological parameters, such as potency and efficacy, that are used in the drug development process to optimize efficacy and safety margin. For agonist ligands, the pre-eminent pharmacological model is the operational model of agonism. In this model, agonist efficacy is quantified using a single term, which accounts for receptor reserve (the transducer ratio τ). GPCR-specific models include the ternary complex model, and its extensions that incorporate constitutive activity (De Lean et al., 1980; Samama et al., 1993; Weiss et al., 1996).

CR IP T

GPCR models are now based on the structural concept of conformational ensembles, which can be selected and/or induced to yield ligand-specific and system-specific receptor signaling and receptor regulation events (Kenakin, 2002; Kenakin, 2015; Manglik and Kobilka, 2014).

Recently, the temporal dimension of GPCR pharmacology has attracted considerable attention. Kinetic

AN US

pharmacology has been advocated as an approach to develop more effective therapeutics (Copeland, 2016; Copeland et al., 2006; Guo et al., 2014; Hoffmann et al., 2015). Knowledge of the duration of receptor occupancy by ligand (its residence time) can be necessary to account for the magnitude of drug effect

M

(Copeland, 2016; Copeland et al., 2006; Guo et al., 2014; Hoffmann et al., 2015; Hothersall et al., 2016). A long residence time offers the potential to increase the duration of therapeutic efficacy beyond that predicted by the

ED

pharmacokinetic profile of the drug (Lacourcière and Asmar, 1999; Vauquelin and Van Liefde, 2006). In principle, this property could allow once-a-day dosing of compounds that are rapidly cleared from the body, a

PT

scenario that could increase the therapeutic window (Cusack et al., 2015). The timing of agonist response is a

CE

determinant of therapeutic efficacy. For long-acting beta agonists, a long duration of exposure is necessary to produce sustained bronchodilation (Zafar et al., 2014). For parathyroid hormone receptor agonists used to treat

AC

osteoporosis, intermittent exposure promotes bone formation whereas continuous administration can result in bone loss (Frolik et al., 2003). The temporal profile of receptor signaling influences the measurement of biased agonism (Klein Herenbrink et al., 2016). Biased agonism holds the promise of activating the therapeuticallyrelevant response while minimizing potentially deleterious signaling through the same receptor (Kenakin, 2015; Whalen et al., 2011). In cell systems, the differences in time course of the responses being compared results in

4

ACCEPTED MANUSCRIPT estimates of signaling bias being dependent on the time at which the response is measured (Klein Herenbrink et al., 2016). At present, pharmacological models for GPCRs do not incorporate kinetics. Typically, time is not a factor in the model or the equations describing it. By contrast, for enzymes and ion channels, time is central component of the models, yielding kinetic drug parameters such as the Vmax for enzymes, and the ligand residence time and the desensitization rate of ion channels (Hodgkin and Huxley, 1952; Michaelis et al., 2011).

CR IP T

For GPCRs, a kinetic model of functional response would be valuable for quantifying the kinetic pharmacology of GPCRs. In some instances, there is a pressing need for models that can be used to analyze kinetic functional data. For example, biosensors are being used routinely to measure the time course of GPCR signaling. This technology has been developed from research into the kinetics of GPCR signaling in live cells (Lohse et al.,

AN US

2008; Sungkaworn et al., 2017; Vilardaga et al., 2009). Presently, no pharmacological model is available to analyze these data. Consequently, the maximum value of this technology might not be being realized. At present, models of GPCR response are typically used to simulate system behavior, rather than to provide

M

estimates of pharmacological activity by fitting curves to data (Riccobene et al., 1999; Shea and Linderman, 1997; Waelbroeck, 2001; Woodroffe et al., 2009). These models provide considerable insight into mechanisms

ED

of signaling but are not designed to provide estimates of macroscopic pharmacological parameters like efficacy and potency.

PT

The goal of this study was to develop a kinetic model of GPCR response that meets the following

CE

criteria:

AC

1. Can be applied to typical kinetic data sets, specifically those in which time and ligand concentration are the only independent variables.

2. Minimal number of fitted pharmacological parameters, for the benefit of structure-activity analysis and to optimize curve-fitting. 3. Equations formatted for use in familiar curve-fitting programs, such as Prism (GraphPad Software, Inc. La Jolla, CA). Prism templates are supplied in the Supplemental information that can be immediately 5

ACCEPTED MANUSCRIPT applied to experimental data. The files contain simulations that are ready to use; all that is required is to modify the parameters of interest. 4. Modular framework, to allow straightforward extension of the model to incorporate additional signaling components and dimensions of pharmacological activity.

In this article, a model is described meeting these criteria. It is a minimal model designed to be extended

CR IP T

in subsequent studies, for example to incorporate receptor desensitization. Practical application of the model to experimental data is emphasized in the body of the article. The mathematical development and formulations are

AC

CE

PT

ED

M

AN US

detailed in the Appendices.

6

ACCEPTED MANUSCRIPT 2. The Model 2.1. Framework In developing a generalized pharmacological model, it is necessary to base the model on principles common to the receptors of interest. A common principle in GPCR signal transduction is intermolecular interaction, between components in the signal transduction cascade. In the first step, the receptor interacts with proteins – G-proteins, arrestins and kinase enzymes. The subsequent cascade is typically one of interaction

CR IP T

between pairs of signaling molecules. For example, Gs interacts with adenylyl cyclase, which generates cAMP, which activates protein kinase A, which leads, through a series of steps, to interaction between cAMP responseelement binding protein and a target sequence on DNA, leading to expression of a protein. Collectively these components are described as a signal transduction pathway. In the operational model, the signal transduction

AN US

pathway is quantified using a single parameter, Em, defined as the “Maximum effect in a system” (Black and Leff, 1983). In the kinetic models formulated here, the transduction system is divided into two temporallyseparated, quantitative components. The first component is the quantity of the transduction system accessible to

M

the receptor before addition of the agonist. The second component is the quantity of the response generated by the agonist-occupied receptor. The amount of response is the amount of the transduction system that is

ED

converted to the response by the agonist-occupied receptor. The first component is the amount of the transduction system material available for generating the response. It is termed here the transduction potential,

PT

by analogy with the concept of potential energy in thermodynamics, and is denoted EP. Explicit examples of

CE

transduction potential and the corresponding response are: The amount of inactive G-protein and active Gprotein; the concentration of Ca2+ in intracellular stores and in the cytoplasm; the level of non-phosphorylated

AC

and phosphorylated proteins; the amount of a transcription factor and the amount of newly-synthesized protein. Agonist binding to the receptor results in conversion of the transduction potential to the response. The

rate of this response-generating process defines the efficacy of the agonist. It is modeled operationally, i.e. with an empirical parameter. This parameter is a rate constant, termed here the transduction rate constant, and denoted kτ. (The script τ is used because the transduction rate constant is a kinetic corollary of τ, the transducer ratio, in the operational model – see Appendix B.5.1.). kτ is the maximum rate of response generation by the 7

ACCEPTED MANUSCRIPT agonist in the system. It is not necessary to know the rate-limiting molecular interaction in the signaling pathway that determines this rate constant. Although the molecular interactions are well established for most GPCR signaling cascades, for example the cAMP pathway, identifying the rate limiting step is not straightforward (Frace et al., 1993). The model is not intended to measure the rate constant of a known molecular interaction in the signaling pathway. In a kinetic model of response, it is necessary to consider the kinetics of agonist binding to the receptor.

CR IP T

In this study, the model framework is divided into two sections, depending on the relative rates of agonistreceptor interaction and response generation. In the first section, the agonist is assumed to rapidly equilibrate with the receptor, i.e. is effectively at equilibrium throughout the time course of the response measurement

AN US

(Section 3). In the second section, receptor-ligand equilibration is considered sufficiently slow that it slows generation of the response (Section 7). An experimental method for distinguishing these models is presented in Section 3.5. The time required to approach equilibrium is a function of the receptor residence time; 90% of the equilibrium occupancy value is reached within a time interval equal to 2.3 times the residence time. Rapid

M

equilibration is typically observed with lower potency agonists (in the M potency range), for which the residence time is typically less than one minute. For example, dopamine and ropinirole at the D 2 dopamine

ED

receptor [22 sec and 19 sec, respectively, (Klein Herenbrink et al., 2016)], and less than 10 sec for a range of

PT

lower potency agonists at the M3 muscarinic acetylcholine receptor (Sykes et al., 2009). Slow equilibration, resulting from a residence time of the order of minutes to hours, is commonly-encountered when using high-

CE

affinity synthetic agonists (Guo et al., 2012; Hothersall et al., 2017; Klein Herenbrink et al., 2016; Rosethorne

AC

et al., 2016) and endogenous peptide and protein ligands (Ferrandon et al., 2009; Hothersall et al., 2016). A second common principle of GPCR-mediated responses is that responses decay over time. For

example, GTP-bound G-protein is converted back to the inactive form, GDP-bound G-protein, by hydrolysis of GTP by the intrinsic GTPase activity of the G-protein. cAMP is metabolized by phosphodiesterase enzymes. Calcium released from intracellular stores is re-sequestered. Synthesized protein is degraded. This principle is incorporated into the model as an exponential decay process (Section 2.2.). A third commonly-encountered 8

ACCEPTED MANUSCRIPT principal is receptor desensitization. Receptor desensitization is not included in the present model. Hence the model in the present study is applicable only to response systems without desensitization. In principal, the model can be extended to incorporate receptor desensitization, as a reduction of the receptor concentration. This is currently being evaluated.

2.2. Formulating the model

CR IP T

In kinetic functional assays, response is measured at multiple time points at multiple concentrations of agonist. In this study, equations were derived that can be used to analyze these data sets to provide estimates of agonist efficacy and potency. The equation format is of a type that can be loaded into commonly-used curvefitting software, such as Prism (GraphPad Software, Inc., La Jolla, CA). Specifically, equations are of the form

AN US

y = f(t), where y is response, and f(t) is a function of time and pharmacological parameters. Throughout this

CE

PT

ED

M

study, this form is described as an E vs t equation. The model is described by the following general scheme:

AC

The response (E) is generated by conversion of the transduction potential (EP) by the agonist-occupied receptor (RA). The forward process, the generation of the response, is described by the law of mass action (Hill, 1910; Langmuir, 1918). Specifically, the rate of response generation is a product of the interacting components and a rate constant. The interacting components are RA and EP. The rate constant is termed here the response generation rate constant, is denoted kE, and is of units of receptor units.t–1. kE is the agonist-dependent parameter in the model, being the rate constant of response generation by the agonist-occupied receptor. The response 9

ACCEPTED MANUSCRIPT decay process is assumed to follow the first order exponential decay function and so is defined by the product of E and a rate constant. The rate constant is termed the response decay rate constant, is denoted kD, and is of units t–1.

2.3. Efficacy in the kinetic models

CR IP T

In the kinetic models, a simple efficacy parameter emerges, termed the transduction rate constant (kτ). It can be described as the maximum rate of response generation by the agonist within the system. kτ is straightforward to measure, as described in Section 3. All that is required is measurement of response over time at a maximally-stimulating agonist concentration. Agonists with different efficacies will manifest different

AN US

values of kτ. The parameter accounts for receptor reserve, like τ in the operational model, because it incorporates the receptor concentration (Appendix C.2.). However, it is not necessary to manipulate the

ED

M

signaling system to estimate kτ, as is often the case for measuring τ (Appendix C.2.). kτ is defined as,

kτ comprises agonist-dependent and system-dependent parameters. kE is the agonist-dependent parameter. kE is

PT

also system-dependent because it is dependent on the particular transduction pathway to which the receptor

CE

couples. The other parameters are system-dependent. The receptor concentration, [R]TOT, can differ between different systems (for example, different cell or tissue types), as can the amount of transduction potential,

AC

EP(TOT)

3. Models with rapid agonist-receptor binding kinetics Fitting the E vs t data to the model equations enables estimation of the pharmacological parameters. In the following sections, curve fitting procedures are shown for models in which agonist rapidly equilibrates with 10

ACCEPTED MANUSCRIPT the receptor. (The slow agonist equilibration is described in Section 7.). The selection of the model and corresponding equation is dependent of the shape of the E vs t profile. The models give rise to E vs t profiles that are commonly observed in pharmacological systems. These include the straight line, the horizontal exponential curve, and the rise-and-fall curve. (The horizontal exponential curve is defined by the equation y = constant.(1 – e–k.t). A familiar example is the ligand binding association curve.) The profile is dependent on the model variant (Appendix A). The profile is dependent on the fate of E, specifically whether or not it decays, and

CR IP T

if it does whether or not it recycles to EP. The profile is also dependent on whether or not EP is depleted by conversion to E. In the following sections, the model mechanism that gives rise to each of these E vs t curves is presented. The curve fitting method is then described. The equations are derived in Appendix B.

AN US

3.1. Model 1: Linear E vs t profile

The simplest E vs t plot observed in pharmacological assays is the straight line. In the kinetic model, this profile emerges when EP is not depleted and E does not decay (Model 1, Appendix B.1.). The model is represented by

AC

CE

PT

ED

M

Scheme 1:

Scheme 1

The time course profile predicted by this model is shown in the simulation in Fig. 1. In the E vs t plot (Fig. 1B) it is evident that the time course is linear at all agonist concentrations and the gradient of the line increases as the agonist concentration increases. Kinetic data are often graphed as an agonist concentration-response curve 11

ACCEPTED MANUSCRIPT over multiple time points (Fig. 1C). When graphed in this way, the model predicts an increase of maximal response as time increases (Fig. 1C). Maximal response increases without limit as time increases. There is no shift of the dose response over time (Fig. 1C). This is evident when plotting EC50 vs time (Fig. 1D); there is no change of EC50 as time increases. Finally, the model predicts there is no receptor reserve in the system. This is discussed in Appendix C.1. The E vs t data can be fit to equations to obtain estimates of the pharmacological parameters. These

CR IP T

parameters are agonist efficacy and affinity, kτ and KA respectively. The parameters can be estimated by globally fitting E vs t data at multiple agonist concentrations to the model equation, Eq. 1 (derived in Appendix

AN US

B.1.):

Eq. 1

M

The fit is shown in Fig. 2A; simulated data were fit to Eq. 1 using GraphPad Prism 7. The fit details are shown

ED

in the legend to Fig. 2 and the Prism file is available in the Supplementary material accompanying this article (“Fig 2 Model 1 global fit”). The Prism file can be used as a template by investigators wishing to employ the

PT

model.

A more familiar method for analyzing the data is also available. The mathematical basis of this method

CE

is given in Appendix B.1. The E vs t data for each concentration of agonist are fit to a straight line equation E =

AC

slope  t. The simulated data in Fig. 2B were fit to the equation using the “Line through origin” Equation in Prism 7.0. The fitted value of the slope is then plotted against the agonist concentration (Fig. 2C). In the next step, these data are fitted to the sigmoid curve equation,

op

og

og

12

i

ope

ACCEPTED MANUSCRIPT This can be done in Prism using the “log(agonist) vs. response -- Variable slope (four parameters)” equation, with the parameter “Bottom” set to zero. The y value at maximally-effective agonist concentrations (“Top” in the equation above) is kτ, and the A50 value is KA. This data analysis method is shown in the Prism file in Supplementary information, “Fig 2 Model 1 individual straight line fits.” With respect to the Hill Slope parameter, in theory the value is unity for the model as formularized in Appendix B. Allowing the Hill Slope to vary accommodates experimental variability in the precision of serial dilution of the agonist, or biological

CR IP T

processes that can deviate the concentration range in the vicinity of the receptor from the concentration range added to the assay. However, a slope value deviating from unity can also result from biological cooperativity mechanisms within the system. Whether cooperativity mechanisms compromise the model assumptions will need to be considered, as is the case for any pharmacological model in which a cooperative mechanism is not

AN US

explicitly incorporated.

An especially convenient method is available for measuring kτ. The time course of response is measured for a maximally-stimulating concentration of agonist. The time course data are fit to a modified form of Eq. 1,

ED

M

derived as the limit of Eq. 1 as [A] approaches infinity:

PT

This fit provides an estimate of kτ. Alternatively, the time course data can be fit to the general straight line

AC

CE

equation, in which kτ is the slope of the line.

3.2. Model 2: Horizontal exponential E vs t profile, t1/2 independent of [A] In this profile, the E vs t curve is a horizontal exponential curve. Initially the response rises steeply over

time and then shallows, finally approaching a plateau at which response does not change over time (Fig. 3B). Model 2 gives rise to a profile in which the t1/2 is the same for all agonist concentrations. In this model, EP is not depleted and E decays (Appendix B.2.). This model is represented by Scheme 2: 13

Scheme 2

CR IP T

ACCEPTED MANUSCRIPT

AN US

The level of response at the plateau, formally the horizontal asymptote, is the response at steady-state. At steady state, the rate of response generation is equal to the rate of response decay. There is no receptor reserve in this model, as discussed in Appendix C.1. The full E vs t profile predicted by Model 2 is shown in Fig. 3 and

ED

M

was simulated using Eq. 2 (derived in Appendix B.2.):

CE

PT

Eq. 2

In the E vs t plot, it is evident that the time course is a horizontal exponential curve (Fig. 3B). As the

AC

concentration of agonist increases, the t1/2 does not change (Fig. 3B). Fig. 2C shows the data graphed as an agonist concentration-response curve over multiple time points. In this plot, Model 2 predicts an increase of the maximal response to agonist over time (Fig. 3C). The maximal response approaches a limit as time increases (Fig. 3C). The EC50 of the agonist does not change over time. This is evident from the lack of a leftward or rightward shift of the agonist concentration response curve (Fig. 3C) and from plotting the EC50 vs time (Fig. 3D). 14

ACCEPTED MANUSCRIPT kτ and KA can be estimated by fitting the E vs t data to Eq. 2. E vs t data for all agonist concentrations are fit globally. An example is shown in Fig. 4A, in which simulated data were fit to Eq. 2 using Prism 7.0. The fit details are shown in the legend to Fig. 4 and in the Prism file in Supplementary information, “Fig 4 Model 2 global fit.” This fitting procedure also provides an estimate of the response decay rate constant, kD. A more familiar method for analyzing the data is also available. The mathematical basis of this method is given in

horizontal exponential equation:

te u

CR IP T

Appendix B.2. The E vs t plot for each concentration of agonist is fit individually to the following general

AN US

This can be done in Prism using the “One-phase association” equation, with the parameter “Y0” set to zero. The resulting fit parameters are the rate constant, K, and the response at the plateau, (i.e. the y value at the horizontal asymptote). kτ and KA can be determined from these parameters as follows. First, the rate constant is multiplied

M

by the plateau. Next, the resulting value is plotted against the agonist concentration. This plot is shown in Fig. 4B. These data are then fit to the sigmoid curve equation, as described for Model 1 in Section 3.1. The y value

ED

at maximally-effective agonist concentrations is kτ, and the midpoint A50 value is KA (Fig. 4B). The response decay parameter kD can also be determined. It is the value of K from the fit to the general horizontal exponential

PT

equation. This curve fitting procedure is shown in the Prism file in Supplementary information, “Fig 4 Model 2

CE

individual horizontal exponential fits.”

kτ can be measured from the time course of response to a maximally-stimulating concentration of

AC

agonist using only a single agonist concentration. The time course data are fit to a modified form of Eq. 2, derived as the limit of Eq. 2 as [A] approaches infinity:

15

ACCEPTED MANUSCRIPT Alternatively, the time course data can be fit to the general horizontal exponential equation, and kτ calculated as the rate constant multiplied by the plateau.

3.3. Model 3: Horizontal exponential E vs t profile, t1/2 dependent on [A] Model 3 applies to responses in which the E vs t profile is a horizontal exponential curve and the t1/2 of the curve is dependent on the agonist concentration. Specifically, the t1/2 decreases as the agonist concentration

CR IP T

increases. This profile results from the kinetic model when EP is depleted by its conversion to E, and E recycles back to EP (Appendix B.3.). This mechanism is encountered in G-protein activation. The response, GTP-bound G-protein, is converted back to the transduction potential, GDP-bound G-protein, by the response-degrading mechanism, hydrolysis of bound GTP to GDP by the intrinsic GTPase activity of the G-protein (Gilman, 1987).

AN US

Receptor reserve is present in Model 3. It is a consequence of depletion of transduction potential, as described

AC

CE

PT

ED

M

in Appendix C.1. The mechanism is defined in Scheme 3 below:

Scheme 3

The full E vs t profile predicted by Model 3 is shown in Fig. 5 and was simulated using Eq. 3 (derived in Appendix B.3.):

16

ACCEPTED MANUSCRIPT ⁄ ⁄

(

) Eq. 3

CR IP T

where,

At each concentration of agonist, the E vs t profile is a horizontal exponential curve (Fig. 5B). As the

AN US

agonist concentration increases, the plateau (response as t  ∞) approaches a limit. In addition, the t1/2 is reduced. This is evident as a left-ward shift of the midpoint of the E vs t curve (Fig. 5B). The reduction of t1/2 also approaches a limit as the agonist concentration increases. In Fig. 5C, the data are plotted as agonist

M

concentration response curves at multiple time points. The maximal response to agonist increases as time progresses, approaching a limit at later time points. The potency of the agonist increases over time, evident as a

ED

leftward shift of the curve (Fig 5C). This effect is clearly evident when the EC50 is plotted against the agonist concentration (Fig. 5D). In this respect, Model 3 differs from Model 2, in which there is no change of agonist

PT

EC50 over time (Fig. 3D).

CE

The pharmacological parameters for the agonist, kτ and KA, can be estimated from a global fit of the E vs t data for all agonist concentrations to Eq. 3. An example is shown in Fig. 6A, in which simulated data were fit

AC

to Eq. 3 using Prism 7.0. The fit details are shown in the legend to Fig. 6 and the Prism file is available in the Supplementary material (“Fig 6 Model 3 global fit”). This fitting procedure also provides an estimate of kD and of EP(TOT). A more familiar method for analyzing the data is also available. The mathematical basis of this method is given in Appendix B.3. This method is the same as that for Model 2 (Section 3.2.). The E vs t plot for each concentration of agonist is fit individually to the general horizontal exponential equation:

17

ACCEPTED MANUSCRIPT te u

The resulting parameters are the rate constant K and the plateau. The rate constant is multiplied by the plateau. The resulting value is then plotted against the agonist concentration (Fig. 6B). These data are then fit to the sigmoid curve equation in Section 3.1. The y value at maximally-effective agonist concentrations is kτ, and the midpoint A50 value is KA (Fig. 6B). This procedure is detailed in the Prism file in Supplementary information,

CR IP T

“Fig 6 Model 3 individual horizontal exponential fits.” kD can also can be determined from the fit to the general equation: The rate constant is plotted against the agonist concentration. kD is the value of the rate constant at the lower asymptote, i.e. the value as [A] approaches zero.

kτ can be estimated from the time course of response to a maximally-stimulating concentration of

AN US

agonist. The time course data are fit to a modified form of Eq. 3, derived as the limit of Eq. 3 as [A] approaches infinity:

(

)

)

ED

M

(

This fit provides an estimate of kτ. Alternatively, time course data for a maximally-stimulating agonist

PT

concentration can be fit to the general horizontal exponential equation, and kτ calculated as the rate constant

AC

CE

multiplied by the plateau. (Multiplying the plateau by the rate term in the above equation gives kτ.)

3.4. Model 4: Rise-and-fall exponential profile In this E vs t profile, the response initially rises rapidly over time, and then shallows and reaches a peak.

Subsequently the response decreases over time, declining to zero at infinite t (Fig. 7A). This profile emerges from the kinetic model when EP is depleted and the response decays (Appendix B.4.). A familiar example of this mechanism is release and subsequent sequestration of cytosolic Ca2+ into the endosomal compartment. 18

ACCEPTED MANUSCRIPT Receptor reserve is present in Model 4, being a consequence of depletion of transduction potential, as described

AN US

CR IP T

in Appendix C.1. The mechanism is defined in Scheme 4 below:

Scheme 4

The full E vs t profile predicted by Model 4 is shown in Fig. 7 and was simulated using Eq. 4 (derived in

ED

M

Appendix B.4.):

PT

Eq. 4

AC

CE

where,



19

ACCEPTED MANUSCRIPT At each concentration of agonist, the E vs t profile is a rise-and-fall exponential curve that approaches zero at later time points (Fig. 7B). The rate of rise increases as the agonist concentration increases, as does the peak response. Both of these effects approach a limit at maximally-effective agonist concentrations (Fig. 7B). At later time points, a surprising phenomenon is apparent when comparing the response at different agonist concentrations. For example, at the 10 min time point, the response to the higher concentrations (46 – 4,600 nM) is less than the response to an intermediate concentration (10 nM). In other words, the dose response at 10

CR IP T

min is bell-shaped. This is clearly evident when the data are plotted as agonist concentration-response curves at various times (Fig 7C, see the curve for the 10 min time point). This phenomenon has been observed experimentally in Ca2+ mobilization assays (Princen et al., 2003). The change of potency over time was investigated at the earlier time points, for which the data are described by a sigmoid curve (Fig. 7D). Potency

AN US

increases over time between 0.3 and 4 min; the curve is left-shifted as the measurement time progresses (Fig 7C) and the EC50 decreases (Fig. 7D).

The pharmacological parameters for the agonist, kτ and KA, can be estimated from a global fit of the E vs

M

t data for all agonist concentrations to Eq. 4. An example is shown in Fig. 8A, in which simulated data were fit to Eq. 4 using Prism 7.0. The fit details are shown in the legend to Fig. 8 and the Prism file is available in the

ED

Supplementary material (Fig 8 Model 4 global fit.”). This fitting procedure also provides an estimate of kD and of EP(TOT). An alternative method for analyzing the data is also available. The mathematical basis of this method

PT

is given in Appendix B.4. The E vs t plot for each concentration of agonist is fit individually to a general rise-

CE

and-fall exponential equation. This equation is less familiar to pharmacologists than other exponential equations such as the association equation and the exponential decay equation. It is commonly-encountered in

AC

pharmacodynamics, being the equation that describes systemic drug concentration in a single-compartment drug absorption and elimination model (Mayersohn and Gibaldi, 1970). It can be written as,

20

ACCEPTED MANUSCRIPT where C is a fitted parameter and K1 and K2 are rate constants. Obtaining kτ and KA from fitting the data to this equation is surprisingly easy. The fitted value of C is plotted against the agonist concentration (Fig. 8B). These data are then fit to the sigmoid curve equation in Section 3.1. The C value at maximally-effective agonist concentrations is kτ, and the midpoint A50 value is KA (Fig. 8B). This procedure is provided in the Prism file, “Fig 8 Model 4 individual rise-and-fall exponential fits” in Supplementary information. The value of kD also can be determined, from the fit to the general rise-and-fall exponential equation – it is the value of the rate constant

CR IP T

that is the same at all agonist concentrations.

It is possible to measure kτ using only a maximally-stimulating concentration of agonist. The time

AN US

course data are fit to a modified form of Eq. 4, derived as the limit of Eq. 4 as [A] approaches infinity:

where,

ED

M



This fit provides an estimate of kτ, EP(TOT) and kD. Alternatively, the time course data can be fit to the general

PT

rise-and-fall exponential equation given above. kτ is the fitted value of C. In rise-and-fall responses, investigators typically use the response value at the peak of the E vs t curve in

CE

quantifying drug effect. This approach greatly simplifies data analysis because there is no need to fit a curve to the E vs t data. Unfortunately, it is not possible to obtain estimates of kτ and KA by plotting the peak response

AC

against the agonist concentration, as described in Appendix B.4.

4. Experimental examples The validity of the models was assessed by fitting the model equations to experimental data. For Model 1, a linear E vs t profile has been observed for glucagon-stimulated cAMP accumulation in hepatic membranes 21

ACCEPTED MANUSCRIPT (Rodbell et al., 1974). The experimental data are reproduced in Fig. 9. The linear profile is consistent with the molecular properties of the system. The response, cAMP generation by adenylyl cyclase, did not decline, likely because the membrane preparation contained little cyclic nucleotide phosphodiesterase activity (Pohl et al., 1969). The substrate did not deplete because an ATP-regenerating system was present (Rodbell et al., 1974). Depletion of GDP-bound G-protein by generation of active G-protein was likely minimal owing to the endogenous level of receptor expression and the short incubation period. These data were fit to Model 1, using

CR IP T

Eq. 1 (red dotted lines in Fig. 9). The fitted values were: affinity (KA), 0.76 nM; and efficacy (kτ) 22 pmol cAMP / mg protein / min. A second example of Model 1 is β2 adrenoceptor-mediated cAMP response element gene transcription, in which the response to agonists increases in a linear fashion up to five hours, with no change of the agonist EC50 (Baker et al., 2004).

AN US

For Model 2, the example is provided by new experiments performed measuring the time course of cAMP generation in whole cells expressing the β2 adrenoceptor (Fig. 10, Table 2). The response to moderatepotency agonists ritodrine and salbutamol was measured as described in the legend to Fig. 10. The E vs t profile

M

was a horizontal exponential curve, implicating a response decay process (Section 3.2.), for example metabolism of cAMP by phosphodiesterase enzymes. When fit to a general horizontal exponential equation

ED

(Section 3.2.), the observed rate was similar for the two effective concentrations, consistent with Model 2 (0.21  0.09 and 0.19  0.11 min-1 for 200 nM and 20 M ritodrine, respectively; 0.14  0.03 and 0.17  0.01 min-1

PT

for 200 nM and 20 M salbutamol, respectively, mean  range, n = 2). The data were then fit globally to Eq. 2,

CE

with a baseline constant added, using Prism 7.0 (Fig. 10). In all cases, the fit was not improved with the use of Model 3 (Eq. 3), tested using an extra sum-of-squares F-test using Prism 7.0, (Motulsky, 2017), indicating

AC

Model 2 was the simplest of the horizontal exponential curve models to fit the data. The fitted parameter values (mean and range, n = 2) are shown in Table 2. The kτ values for ritodrine and salbutamol were similar (7.0 and 5.1 nM cAMP / min respectively, Table 2). This finding is consistent with the similar efficacies of these agonists for stimulating cAMP at single time point, 30 min (Baker, 2010). The agonist affinity value from the fits to Model 2 (KA) was 170 nM for ritodrine and 82 nM for salbutamol (Table 2). These values are approximately 10-fold lower than the Ki value measured using a whole cell radioligand binding assay [1,500 22

ACCEPTED MANUSCRIPT nM and 980 nM for ritodrine and salbutamol, respectively (Baker, 2010)]. This difference could be explained by the assays sampling different distributions of the GPCR conformational ensemble (Kenakin, 2002; Motulsky et al., 1985). Regarding kD, the range of kD values for the two agonists overlapped (Table 1), with a mean value of 0.30 min-1 for ritodrine and 0.17 min-1 for salbutamol. The similarity of the kD value for the two agonists is consistent with the model because the response decay mechanism in the model is downstream of the receptor (Scheme 2, Section 3.2.).

CR IP T

An example for Model 3 is the results of DAMGO-stimulated [35S]GTPS binding to membranes from cells expressing the -opioid receptor (Traynor et al., 2002). The GTPS binding assay provides a convenient measure of upstream agonist response (Strange, 2010). The steady-state concentration was dependent on DAMGO concentration (Traynor et al., 2002). In addition, the observed rate was dependent on DAMGO

AN US

concentration, and the observed rate at maximally-stimulating concentrations of agonist was dependent on the efficacy of the agonist. These observations are consistent with the E vs t equation of Model 3 (Eq. 3, Section 3.3.). A data simulation with Eq. 3 was performed to obtain values of the observed rate and steady-state

M

response. The results are compared with the experimental results in Table 3. Both the simulated rate and steadystate response were reasonably consistent with the experimental data. In the simulation, the affinity of the

ED

agonist DAMGO was 130 nM, and kτ was 0.016 pmol/mg/min. It is possible to assess the degree of receptor reserve in the response, by calculating the value of τ, the transducer ratio of the operational model. The fitted

AC

CE

PT

values of kτ, EP(TOT) and kD are entered into the following equation, derived in Appendix B.5.1:

The calculated τ value was 1.2. This value indicates a small degree of receptor reserve; the EC50 is 2.2-fold lower than the KA value. Only slight receptor reserve is anticipated in [35S]GTPS binding assays because the response is only one step downstream of the receptor.

23

ACCEPTED MANUSCRIPT An example for Model 4 is bombesin-stimulated Ca2+ mobilization via the BB1 bombesin receptor (Fig. 10) expressed in COS-7 cells (Liu et al., 2003). Data are shown for a maximally-stimulating concentration of bombesin (1 M). The kτ value was estimated by fitting data to the modified form of Eq. 4 that assumes a maximally-stimulating concentration of agonist (Section 3.4.). The estimated value of kτ was 1,500 FIU.sec-1. The degree of receptor reserve can be determined by calculating τ, as described for Model 3 in the preceding paragraph. The calculated value was 6.7. This value indicates appreciable receptor reserve. Receptor reserve is

CR IP T

commonly encountered in intracellular calcium responses because the response is three steps downstream of the receptor (receptor > G-protein activation > IP3 generation > calcium release from intracellular stores).

AN US

5. Measuring antagonist-receptor binding kinetics in functional assays

One of the benefits of developing the kinetic agonism model is it provides a framework for quantifying the kinetics of antagonist binding. Currently, ligand binding kinetics for GPCRs is routinely measured using

M

radioligand binding assay. Binding of a tracer ligand is measured over time in the absence and presence of test compound and the data fit to an equation (the Motulksy and Mahan equation) to obtain an estimate of the

ED

compound association and dissociation rate constant (Motulsky and Mahan, 1984). There is currently no equivalent method in functional assays. In this study, equations were derived describing inhibition of agonist

PT

response by antagonist over time. The equations derived here are for responses that conform to two criteria –

CE

the agonist rapidly equilibrates with the receptor, and the response is without receptor reserve. The models without receptor reserve are Models 1 and 2. In the section below, methods are described for Model 2, the

AC

horizontal exponential E vs t profile with t1/2 being independent of agonist concentration. (Equations for Model 1 are given in the Appendix, but equations for Models 3 and 4 were not obtained – see Appendix A.3.). Antagonist is assumed to competitively inhibit agonist binding in a one-step receptor interaction. The scheme for the receptor interactions and generation of response is,

24

ACCEPTED MANUSCRIPT

CR IP T

where B is antagonist, and k3 and k4 the antagonist association and dissociation rate constant, respectively. Table 1 lists the model parameters used in this scheme and in the equations shown in the following sections. Equations were derived for four schedules of ligand addition to the receptor system (Appendix C). These are

AN US

illustrated schematically in Fig. 12. The schedules are:

1. Schedule A: Simultaneous addition of agonist and antagonist. 2. Schedule B: Antagonist pre-incubation

PT

ED

4. Schedule D: Agonist pre-incubation

M

3. Schedule C: Antagonist pre-incubation and washout

5.1. Schedule A: Simultaneous addition of agonist and antagonist

CE

Typically, in a kinetic competition experiment the receptor is exposed to both ligands simultaneously. In the context of a functional assay, the receptor response system, e.g. whole cells, is added to an assay vessel, e.g.

AC

the well of a microtiter plate, that already contains agonist and antagonist. The time course of the response is then measured. This is Schedule A (Fig. 12A). In Fig. 13, simulated data for such an experiment are shown, simulated with random scatter to mimic experimental variability. Response to agonist is measured in the presence of multiple concentrations of antagonist, and in the absence of antagonist. The data are then fit to an equation that provides estimates of the antagonist association and dissociation rate constant. This is Eq. 5, derived in Appendix D.1: 25

ACCEPTED MANUSCRIPT

*

+

Eq. 5

)

AN US

(

CR IP T

where,

E[B]=0,t∞ is the response at the plateau in the absence of antagonist. This equation resembles the familiar

M

Motulsky and Mahan equation for kinetics of competitive binding (Motulsky and Mahan, 1984). For ease of

ED

use, it has been entered as a custom equation in a Prism 7.0 file provided in the Supplemental information accompanying this article (“Fig 12 Antagonist kinetics, simultaneous addition”). The fitting method is described

PT

in the legend to Fig. 13 and in the Prism file. The shape of the curve is dependent on the antagonist dissociation rate constant (Supplementary information file, “Antagonist binding kinetics in functional assays”, Supplemental

CE

Fig. S1). For example, in Fig. 13 an “overshoot” is observed in the time course profile because the antagonist

AC

dissociation rate constant is less than the decay rate constant. Importantly, a limit of detection is approached as the dissociation rate constant (k3) approaches and then exceeds the response decay rate constant (k4) (Supplementary Fig. S1D). Under these conditions, data might be better fitted by an equation that assumes rapid equilibration between antagonist and receptor. This equation is given in Appendix D.1.4.

5.2. Schedule B: Antagonist pre-incubation 26

ACCEPTED MANUSCRIPT Often in functional assays of antagonist activity it is necessary or desirable to add the antagonist prior to the agonist. For example, when using adherent cells in a microtiter plate. Antagonist is pre-incubated with cells for a specific period of time, then agonist is added and the time course of response measured. This is Schedule B (Fig. 12B). An example, using simulated data, is shown in Fig. 14. The response to agonist is measured in the presence of multiple concentrations of antagonist, and in the absence of antagonist. The data are fit to Eq. 7,

CR IP T

derived in Appendix D.2:

+

Eq. 6

AN US

*

where,

(

)

)

(

)

CE

PT

ED

M

(

A Prism 7.0 file is provided in Supplemental information to facilitate application of this equation to new

AC

data (“Fig 13 Antagonist kinetics, antagonist pre-incubation”). The fitting method is described in the legend to Fig. 14 and in the Prism file. The shape of the curve is dependent on the antagonist dissociation rate constant, k4 (Supplementary information file, “Antagonist binding kinetics in functional assays”, Supplemental Fig. S2). A limit of detection is reached as k4 approaches or exceeds kD (Supplementary Fig. S2D). Under these conditions, fitting to an equation that assumes rapid equilibration might be more appropriate than fitting to Eq. 6. This equation is provided in Appendix D.2.4. 27

ACCEPTED MANUSCRIPT

5.3. Schedule C: Antagonist pre-incubation and washout In this schedule, antagonist is pre-incubated with the receptor. Subsequently, unbound antagonist is removed by washing, under conditions in which the receptor-antagonist complex is preserved and the antagonist concentration is reduced to zero. Agonist is then added and the time course of the response measured (Fig. 12C). This paradigm can be applied to ex vivo receptor occupancy assays, used to evaluate target engagement

CR IP T

of therapeutic candidates. It is the functional analogue of a ligand binding method described previously (Malany et al., 2009). The antagonist is administered to the animal for a specified period of time during which it occupies the receptor. The relevant tissue is then dissected, the preparation washed to remove unbound antagonist, and the response to agonist measured over time. These data can be fit to Eq. 7. The curve fitting procedure provides

AN US

an estimate of the degree of receptor occupancy by antagonist. The analysis also provides an estimate of the

M

antagonist dissociation rate constant. Eq. 7, derived in Appendix D.3, is:

PT

ED

where,

Eq. 7

CE

ρB(t=0) is fractional receptor occupancy by antagonist at the moment of agonist addition.

AC

A simulated example of the data and the curve fitting procedure is shown in Fig. 15. A test compound is applied to an animal at the doses shown in the figure. A vehicle control is included, i.e. without antagonist. Data for the vehicle control are analyzed to obtain an estimate of kD, using the following equation (Appendix D.3.3):

28

ACCEPTED MANUSCRIPT Data for the antagonist-treated samples are then analyzed using Eq. 7, as described in the legend to Fig. 15 and in the Prism 7.0 file in Supplementary information, “Fig 14 Antagonist kinetics, antagonist pre-incubation and wash.” The resulting receptor occupancy values can then be plotted against the dose (see Fig. 15, inset). From this occupancy vs dose curve the ED50 for receptor occupancy can be determined. The data analysis method reaches a limit of detection as the antagonist dissociation rate constant approaches or exceeds the response decay rate constant (Supplementary information file, “Antagonist binding kinetics in functional assays”,

CR IP T

Supplemental Fig. S3D). This is an important consideration. If the antagonist dissociates too rapidly, it will not affect the response to agonist and it would be erroneously concluded that antagonist does not engage the target in vivo.

AN US

5.4. Schedule D: Agonist pre-incubation

In Schedule D, receptor is first exposed to the agonist, generating a response. Antagonist is then added, without removal of agonist, and the response measured (Fig. 12D). This paradigm is designed to be used when

M

the onset of the response to agonist is rapid, such that response at early time points is difficult to measure accurately or cannot be measured at all. The analysis assumes the response to agonist has reached steady-state

ED

before application of the antagonist. Response to agonist is measured in the presence of multiple concentrations

PT

of antagonist, and in the absence of antagonist. The data are then fit to Eq. 8, derived in Appendix D.4:

]

AC

CE

[

Eq. 8

The fitting method is described in the legend to Fig. 14 and in a Prism 7.0 file provided in Supplemental information (“Fig 15 Antagonist kinetics, agonist pre-incubation”). This file can be used as a template to analyze new data. The shape of the curve is dependent on the antagonist dissociation rate constant, k4 29

ACCEPTED MANUSCRIPT (Supplementary information file, “Antagonist binding kinetics in functional assays”, Supplemental Fig. S4). For rapidly-dissociating compounds (as k4 approaches or exceeds kD), fitting to an equation that assumes rapid equilibration might be more appropriate than fitting to Eq. 8 (Supplementary Fig. S4D). This equation is provided in Appendix D.4.4.

6. Testing for rapid vs slow agonist-receptor binding kinetics

CR IP T

In the sections above, agonist is assumed to rapidly equilibrate with the receptor. In the sections below, the alternative scenario is evaluated, that agonist equilibrates sufficiently slowly to affect the kinetics of the response. A method is required to discriminate between rapid and slow agonist-receptor equilibration. Here a straightforward test is described. The simulated data are shown in Fig. 17.

AN US

In this experiment, receptor is exposed to agonist for a period of time during which response is generated (the stimulation phase). At a specified time point, a large excess of a competitive antagonist is added, sufficiently large that the only RA complexes in the system are those that formed during the stimulation phase.

M

After the addition of antagonist, the level of response continues to be measured over time (the antagonist block phase). The theoretical basis of the test is as follows. If the dissociation of agonist from receptor is sufficiently

ED

rapid, [RA] is reduced to zero rapidly after addition of antagonist, so no further response generation occurs. Under this condition, the only process occurring during the antagonist block phase is decay of the response. By

PT

contrast, if dissociation of agonist-receptor complex is sufficiently slow, RA will be present in sufficient

CE

abundance to generate additional response during the antagonist block phase. Under this condition, two processes will occur in the dissociation phase, response generation and response decline.

AC

In Fig. 17, the manifestation of these scenarios on the time course of response is presented, for Models 1 to 4. The time of antagonist addition is indicated by the arrow. The solid line represents the response for rapid dissociation of RA, and the dashed line the response for slow dissociation of RA (residence time of 12 min). The fast vs slow kinetic scenarios are clearly differentiated. For Model 1 (Fig. 17A), when RA dissociates rapidly the gradient of the E vs t plot immediately falls to zero after antagonist addition. This makes sense because no more response is generated, and the response that has been generated during the stimulation phase does not decline 30

ACCEPTED MANUSCRIPT because there is no response decay in Model 1. When RA dissociates slowly in Model 1, response continues to increase after antagonist addition (Fig. 17A). This can be explained by the presence of persistent RA complexes that continue to generate the response. At later times, the curve approaches a horizontal plateau because the RA concentration approaches zero, due to complete dissociation of the agonist. In Models 2 and 3, for the rapid equilibration model there is an immediate decline of the response level after the antagonist is added (Fig. 17B and C). The response declines exponentially, defined by the exponential

CR IP T

decay equation (E = E(t=0)e–k.t). The half-time of response decline is equal to the t1/2 of the response decay rate constant, kD (Fig. 17B and C). As a result, the rapid equilibrium assumption can be quantitatively verified, as follows: The exponential decay equation is fitted to the antagonist block phase data. The fitted value of the rate constant k is then compared with the value of kD. (kD measurement is described in Sections 3.1. and 3.2. for

AN US

Models 2 and 3, respectively.) The shape of the E vs t profile in the antagonist block phase make sense. Additional response cannot be generated after antagonist addition, owing to the loss of RA. The response generated during the stimulation phase declines owing to the response decay mechanism in the models. In the

M

case of slow agonist dissociation, there is a delay in the decline of the response (Fig. 17B and C) and as a result the half-time of decline is greater than the kD t1/2 (Fig. 17B and C). This delay results from the presence of

ED

persistent RA complexes during the antagonist block phase. For Model 4, the rise-and-fall E vs t profile, the experiment most clearly discriminates fast vs slow

PT

equilibration when antagonist is added before the peak of the response is reached (Fig. 17D). When RA

CE

dissociates rapidly, the result is the same as that for Models 2 and 3 (see above): There is an immediate, exponential decline, the half-time of which is equal to the t1/2 of the response decay rate constant, kD. (The rate

AC

of response decline is determined from fitting the data in Fig. 17D to the exponential decay equation, and kD is determined as described in Section 3.4.) When RA dissociates slowly, the shape of the E vs t profile is different. The response continues to rise, peaks, then declines (Fig. 17D). If this profile is observed, visual inspection is sufficient to conclude that agonist dissociates slowly from the receptor. The half-time of response decline is greater than the kD t1/2 (Fig. 17D). This profile results from persistent RA complexes generating additional response during the antagonist-block phase. 31

ACCEPTED MANUSCRIPT

7. Models with slow agonist-receptor binding kinetics In this section, the effect of slow agonist-receptor equilibration on the kinetic models is determined. In these models, the time course of the response is determined by two rate processes – the rate of agonist recognition, and the rate of response generation by the agonist occupied receptor. [A kinetic operational model has been presented for the case of rate-limiting agonist association, in which receptor and transduction system

CR IP T

are in equilibrium (Slack and Hall, 2012).] Four aspects are considered. 1) The shape and equations of the E vs t profile. 2) Methods for measuring agonist efficacy. 3) The capacity of the models for measuring agonist binding parameters [affinity, association rate constant (k1) and dissociation rate constant (k2)]. 4) Effect of slow equilibration on agonist EC50. Note in the sections below, Models 5,6,7 and 8 are the slow equilibration

AN US

analogues of Models 1,2,3 and 4, respectively. The E vs t equations are derived in Appendix E.

7.1. Model 5: No EP depletion, no E decay, slow agonist equilibration

M

In this model, the response mechanism is the same as that in Model 1. Transduction potential is not depleted by generation of the response and the response does not decay (Section 3.1.). The E vs t data for the

ED

model with slow equilibration is defined by Eq. 9, derived in Appendix E.1:

PT

(

CE

Eq. 9

AC

where,

)

In Fig. 18, the data profile is simulated using this equation. When agonist equilibrates slowly, the E vs t profile is a concave curve at low concentrations (Fig. 18A) and a straight line at maximally stimulating concentrations (dashed line in Fig. 18A). This pattern differs from that for rapid agonist equilibration, in which the E vs t 32

ACCEPTED MANUSCRIPT profile is a straight line at all concentrations of agonist (Fig. 1B). The agonist binding parameters can be estimated by globally fitting Eq. 9 to time course data obtained at multiple concentrations of agonist, i.e. analysis with two independent variables, t and [A]. (In Supplementary information, a Prism 7.0 file is provided showing the fitting procedure using simulated data, “Model 5 slow agonist equilibration curve fitting”) This global analysis yields estimates of the agonist binding rate constants k1 and k2. A simple method is available for measuring agonist efficacy (kτ). The time course of response is

CR IP T

measured at a maximally-effective concentration of agonist. Under these conditions, the Eq. 9 reduces to,

This is the same equation as that for Model 1 (Section 3.1.). This is explained by the law of mass action; as the

AN US

agonist concentration approaches the maximally-effective concentration, the time taken for equilibrium occupancy to be obtained approaches zero, which is the assumption of the rapid equilibration model (Model 1). This is the equation for a straight line. As a result, kτ can be measured simply, as the gradient of the E vs t plot

M

at a maximally-stimulating concentration of agonist. This plot is the dashed line in Fig. 18A and is shown in the Prism file “Model 5 slow agonist equilibration curve fitting” in Supplementary information.

ED

The time-dependence of EC50 is now examined, by transposing the E vs t data at multiple agonist concentrations (Fig. 18A) to E vs [A] data at multiple time points (shown in Fig. 18B). As time progresses, there

PT

is a leftward shift of the dose-response curve (Fig. 18B), visualized easily in the EC50 vs time plot (Fig. 18C).

CE

This pattern is explained by the slow equilibration of agonist with receptor. At early time points, equilibrium has not been reached so there is a deficit of agonist occupied receptors. Consequently, there is a deficit of

AC

response, which translates to a reduction of agonist potency. As time progresses, equilibrium is more closely approached so the concentration of agonist-occupied receptors increases, resulting in increased agonist potency. Extrapolating to infinite time, the true agonist potency is reached (EC50 of 100 nM, Fig 17C). This profile is in marked contrast to that for a rapidly-equilibrating agonist (Fig. 1D, Section 3.1.), in which the true agonist potency is achieved immediately. As a result, there is no change of agonist EC50 (Fig. 1D). This difference is important when considering receptor reserve (Appendix C.1.). Working under the assumption of rapid 33

ACCEPTED MANUSCRIPT equilibration, a change of EC50 is used to conclude the presence of receptor reserve in the response system. However, a change of EC50 over time could indicate slow equilibration of agonist with receptor. If the latter is true, we could falsely conclude there is receptor reserve in the system. This is one reason why a test to discriminate between slow and rapid equilibration is necessary (Section 6).

7.2. Model 6: No EP depletion, E decay, slow agonist equilibration

CR IP T

In this model, the response mechanism is the same as that for Model 2. Transduction potential is not depleted by generation of the response and the response decays (Section 3.2.). The E vs t data for the model

AN US

with slow equilibration is defined by Eq. 10, derived in Appendix E.2:

(

Eq. 10

ED

M

where,

)

PT

In Fig. 18D, the data profile is simulated using this equation. When agonist equilibrates slowly, the E vs t profile is a bi-exponential curve at low concentrations (Fig. 18D) and a horizontal exponential curve at

CE

maximally-stimulating concentrations (dashed line in Fig. 18D). This pattern differs from that for rapid agonist

AC

equilibration, in which the E vs t profile is a horizontal exponential curve at all concentrations of agonist (Fig. 3B). The agonist binding parameters can be estimated by globally fitting Eq. 6 to time course data obtained at multiple concentrations of agonist, i.e. analysis with two independent variables, t and [A]. (In Supplementary information, a Prism 7.0 file is provided showing the fitting procedure using simulated data.) This global analysis yields estimates of the agonist binding rate constants k1 and k2.

34

ACCEPTED MANUSCRIPT A simple method is available for measuring agonist efficacy (kτ). The time course of response is measured at a maximally-effective concentration of agonist. Under these conditions, the E vs t equation reduces

CR IP T

to,

Note this is the same equation as that for Model 2 in the case of a maximally-stimulating concentration of agonist (Section 3.2.). This is explained by the law of mass action, as described for Model 5 (Section 7.1.) This

AN US

is the equation for a horizontal exponential curve (Section 3.2.):

te u

where kτ / kD is the plateau and kD the rate constant, K. As a result, kτ can be measured as follows: The E vs t

M

data at maximally stimulating concentration of agonist (the dashed line in Fig 8D) are fit to the general

ED

horizontal exponential equation. kτ is then calculated, by multiplying the plateau by the rate constant. This is shown in the Prism file, “Model 6 slow agonist equilibration curve fitting” in Supplemental information.

PT

Agonist potency increases over time, manifest as a leftward shift of the dose-response curve (Fig. 18E) and a decrease of the EC50 (Fig. 18F). This results from a deficit of agonist-occupied receptors before

CE

equilibrium is reached (see Section 7.1.). This profile contrasts with that for rapid equilibration rapidly (Model

AC

2), in which the EC50 is constant over time (Fig. 3D). As described in Section 7.1., these considerations complicate the detection of receptor reserve in the system, necessitating the experiment to discriminate fast vs slow agonist-receptor equilibration (Section 6).

7.3. Model 7: EP depletion, E recycles to EP, slow agonist equilibration

35

ACCEPTED MANUSCRIPT In this model, the response mechanism is the same as that in Model 4. Transduction potential is depleted by generation of the response and the response recycles back to the transduction potential (Section 3.3.). When agonist equilibrates slowly, the E vs t profile resembles a biexponential curve at low concentrations (Fig. 18G) and is a horizontal exponential curve at maximally-stimulating concentrations (dashed line in Fig. 18G). This pattern differs from that for rapid agonist equilibration, in which the E vs t profile is a horizontal exponential curve at all concentrations of agonist (Fig. 5B). An E vs t equation was not obtained for this model (except in

CR IP T

the case of a maximally-stimulating agonist concentration, see below) because the response generation mechanism is second order with respect to time (see Appendix A.3.). In order to model the response behavior, E vs t data were simulated by numerical solution of the differential equations, using the Euler method [(Vauquelin et al., 2001), see Appendix E.3. for details, simulated curves in Fig. 18G].

AN US

Efficacy (kτ) can be measured from the response to maximally-stimulating concentrations of agonist (dashed line in Fig 17G). An E vs t equation can be readily derived for this response (Eq. 11, derived in

(

(

)

)

ED

M

Appendix E.3.). The equation is,

PT

Eq. 11

CE

[Note this is the same as that for rapid agonist equilibration (Section 3.3.), a result of the law of mass action (see

AC

Section 7.1.] This equation takes the form of the general horizontal exponential equation (Section 3.2.):

te u

kτ is determined by multiplying the plateau value by the rate constant, K. (In Eq. 11, multiplying the plateau by the rate term gives kτ.)

his method is shown in the rism fi e “Model 7 curve fitting at maximally-

stimu ting A” in upp ement ry inform tion. 36

ACCEPTED MANUSCRIPT In the simulation of the E vs [A] curves, agonist potency increases over time, manifest as a leftward shift of the dose-response curve (Fig. 18H) and a decrease of the EC50 (Fig. 18I). The EC50 reduction over time is greater than that for the rapid agonist equilibration model (Model 3, compare Fig. 18I with Fig. 5D). This results from a deficit of agonist-occupied receptors before equilibrium is reached (see Section 7.1.).

7.4. Model 8: EP depletion, E decay, slow agonist equilibration

CR IP T

In this model, the response mechanism is the same as that in Model 3. Transduction potential is depleted by generation of the response and the response decays (Section 3.4.). When agonist equilibrates slowly, the E vs t profile resembles a bell-shaped curve at low concentrations (Fig. 18J) and is a rise-and-fall exponential curve at maximally-stimulating concentrations (dashed line in Fig. 18J). Mathematically, this pattern differs from that

AN US

for rapid agonist equilibration, in which the E vs t profile is a is a rise-and-fall exponential curve at all concentrations of agonist (Fig. 7B). However, the distinction is not obvious from visual inspection of the data (comparing Fig. 18J with 7B). An E vs t equation was not obtained for this model (except in the case of a

M

maximally-stimulating agonist concentration, see below) because the response generation mechanism is second order with respect to time (see Appendix A.3.). In order to model the response behavior, E vs t data were

ED

simulated by numerical solution of the differential equations, using the Euler method [(Vauquelin et al., 2001), see Appendix E.4. for details, simulated curves in Fig. 18J].

PT

Efficacy (kτ) can be measured from the response to maximally-stimulating concentrations of agonist

CE

(dashed line in Fig 17J). An E vs t equation can be readily derived for this response (Eq. 12, derived in

AC

Appendix E.4.). The equation is,

(

)

Eq. 12

This equation takes the form of the general rise-and-fall exponential equation (Section 3.4.): 37

ACCEPTED MANUSCRIPT

kτ is determined simply, as the value of C from the curve fit. (By inspection, it is evident that C is equal to kτ in Eq. 12). his method is shown in the rism fi e “Mode 8 curve fitting t m xim

y-stimu ting A” in

CR IP T

Supplementary information. The agonist dose-response curve at multiple time points is shown in Fig. 18K. While at early time points the data conform to a sigmoid curve, the dose response is bell-shaped at higher concentrations. In Fig. 18L, the EC50 from the sigmoid curve is plotted against time (0.3, 0.6, 2 and 4 min time points). The EC50 decreases over

AN US

time (Fig. 18L), as observed in the rapid agonist equilibration model (Fig. 7D).

8. Measuring the agonist residence time in functional assays

Knowing the residence time of an agonist is potentially valuable in interpreting and predicting

M

pharmacological response characteristics. This applies to in vitro response systems, for example biased agonism

ED

(Klein Herenbrink et al., 2016), and potentially to clinical efficacy, for example the S1P1 agonist FTY720-P (a metabolite of fingolimod) for the treatment of multiple sclerosis (Sykes et al., 2014b). In recent studies, the

PT

agonist residence time has been measured by recording the decline of functional response to the agonist by stopping further agonist association with the receptor, by washout or by addition of a sufficiently high

CE

concentration of antagonist [reviewed in (Hothersall et al., 2016), a specific example being the A2A receptor

AC

(Hothersall et al., 2017)]. This assay paradigm can be run in a screening mode, particularly when using a continuous-read functional assay. However, there are currently no mechanistic equations available that can be used to analyze the data to determine the agonist residence time. In this study, such equations were obtained for Models 5 and 6, in which transduction potential is not depleted. Equations were not obtained for Models 7 and 8 because the response generation mechanism is second order with respect to time (see Appendix A.3.).

38

ACCEPTED MANUSCRIPT Here a method is presented for measuring agonist residence time for Model 5. (Model 6 is rarely encountered in the literature. However, the relevant equations are derived in Appendix E.6.4.). The experimental paradigm follows that in (Hothersall et al., 2017). Agonist is applied for a period long enough for the response to approach steady-state (the plateau as t approaches infinity). A high concentration of competitive antagonist is then added, and the response measured over time (see scheme in Fig. 19). The resulting data are

CR IP T

then fit to Eq. 13 (Appendix E.6.3.):

(

)

AN US

where Et=0 is the response to agonist at the time of antagonist addition (on the plateau of the E vs t curve), k2 the agonist dissociation rate constant, and kD the response decay rate constant. In Supplementary information, a Prism file showing the fitting procedure is given “Agonist residence time assay curve fitting”). The fits yield estimates of both k2 and kD. Alternatively, kD can be entered as a constant from a separate assay, for example

M

from measuring the time course of agonist response generation at a maximally-stimulating concentration of

ED

agonist.

The shape of the curve is shown in Fig 18 for varying values of k2, i.e. varying values of the agonist

PT

residence time. For a rapidly-dissociating agonist (“Very short” in Fig. 19), the response decline conforms to a mono-exponential decay function, with the t1/2 being equal to the kD t1/2. In this example, dissociation is assumed

CE

to be so rapid that [RA] declines to zero immediately upon antagonist addition, as described in Section 6. An A2A receptor agonist that dissociates rapidly (assessed in a washout radioligand binding experiment) displays

AC

this profile (compound 1aa, Fig. 3A in (Hothersall et al., 2017)). As the residence time increases up to 30 min (Fig. 19), it becomes increasingly evident that the curve is bi-exponential; there is a lag phase and then the response declines more rapidly. This curve shape has been observed experimentally (for example compound 3cd in (Hothersall et al., 2017)). At relatively long residence times, little decline of the response is evident over

39

ACCEPTED MANUSCRIPT the duration of the time course (120 and 420 min in Fig. 19), as observed experimentally (compound 3ch in

AC

CE

PT

ED

M

AN US

CR IP T

(Hothersall et al., 2017)).

40

ACCEPTED MANUSCRIPT Discussion In this study a new kinetic pharmacological model is presented for quantifying kinetic drug parameters from functional data. The model quantifies efficacy in kinetic terms (the transduction rate constant) and allows receptor binding constants for antagonist and agonist ligands to be measured in functional assays. These parameters are straightforward to measure in conventional kinetic functional assays. The model is based on the known chemical characteristics of GPCR signaling, intermolecular interaction between signaling molecules,

CR IP T

and decline of signal resulting from response decay processes (Section 2.1.). For example, activation of Gprotein resulting from receptor-G-protein interaction, and deactivation of G-protein resulting from hydrolysis of bound GTP to GDP. Introducing a parameter representing the response precursor, termed the transduction potential, permitted straightforward mathematical derivation of the model equations (Appendix B). The model

AN US

is modular, designed to be extended or modified. The model was validated by comparison with experimental data: It yields known general response vs time profiles of GPCR responses (Section 3) and fits specific experimental data (Section 4). The model was also validated theoretically for the particular forms in which

M

transduction potential is depleted (Models 3 and 4.). In these cases, the underlying assumption of the operational model, that response is a hyperbolic function of receptor concentration, is apparent in the kinetic

ED

equations (Appendix B.5.). While developed for GPCRs, the model is potentially applicable to other metabotropic receptor classes that mediate their effects via intermolecular interaction between signaling

PT

pathway components, for example nuclear hormone receptors (Alexander et al., 2017b) and catalytic receptors

CE

(Alexander et al., 2017a).

The modular nature of the model allows it to be extended to incorporate additional components of the

AC

signaling system and other dimensions of pharmacological activity. An example is the inclusion of receptor binding kinetics (Sections 5 – 8), done simply by making [RA] a time-dependent variable in the formulation. Mathematically, the Laplace transform method is well-suited for addition of kinetic components, provided that no second order processes are present in the mechanism. This approach has been used in pharmacokinetic modeling to add compartments in multi-compartment models (Benet and Turi, 1971; Mayersohn and Gibaldi, 1971). Receptor desensitization can be incorporated by allowing the active receptor concentration to change 41

ACCEPTED MANUSCRIPT over time. This is presently being investigated. The model can also be extended to include two (or more) pools of transduction potential. This concept allows for pre-coupling of receptor with response, i.e. before addition of the agonist, for example pre-coupling of receptor with G-protein. The pre-coupled pool would be rapidly transformed to response (quantified as a high value of kE) and the non-pre-coupled pool would be transformed more slowly because it first needs to access the receptor (resulting in a lower value of kE). Multiple transduction potential pools are known to be involved in Ca2+ signaling, for example Ca2+ stored in the endoplasmic

CR IP T

reticulum and present in the extracellular milieu. Two pools of transduction potential can yield a rise-and-fall to plateau E vs t profile, i.e. one in which the response plateaus at a certain level after prolonged exposure to the agonist exposure (data not shown). This profile is frequently observed experimentally (Klein Herenbrink et al., 2016).

AN US

The model is designed to provide estimates of a small number of pharmacological descriptors, using familiar curve fitting techniques, that can be used by medicinal chemists in drug discovery and development projects. These descriptors are the efficacy parameter kτ, and the receptor binding constant(s) (KA for the rapid

M

equilibration model, k1 and k2 for the slow equilibration model). The nature of these parameters is discussed in depth below. The curve fitting is straightforward and familiar to pharmacologists. It can be performed in the

ED

widely-used data analysis software GraphPad Prism (GraphPad Software, Inc., La Jolla, CA). Prism templates are provided in Supplemental information to facilitate rapid application of the equations. (Note specific

PT

recommendations on curve-fitting practices, such as which parameters to share among data sets or hold

CE

constant, are avoided in anticipation that these will be specific to the experiment being conducted.) Two types of curve fitting protocol are described, global fitting to the entire data set or individual fitting to the data for

AC

each agonist concentration. A benefit of the former is that analysis is performed in a single step, but it requires fitting to an unfamiliar, rather complex equation. The benefit of the latter is that data are fit to a familiar, simple equation (for example, the straight line) but an additional step is required, plotting the fitted parameters versus the agonist concentration, to obtain the model parameter estimates. The individual curve fit method also incorporates the Hill slope parameter, as described in Section 3.1. Efficacy is especially straightforward to measure in the kinetic models. Efficacy is defined by a single parameter, kτ. All that is required is measurement 42

ACCEPTED MANUSCRIPT of response over time at a maximally-stimulating agonist concentration. This method can be applied whether or not agonist equilibrates rapidly with the receptor; the E vs t equation is the same for both at a maximallystimulating agonist concentration (Section 7). Finally, the minimal number of parameters should aid application of kinetic functional assays to medicinal chemistry, especially when the testing large numbers of compounds typically encountered in the hit-to-lead and lead-optimization stages of drug discovery. Presently, response values at multiple time points are used as the descriptors of efficacy. Depending on the circumstances, these

CR IP T

values could be replaced with a single efficacy parameter with the application of the kinetic model (kτ), from which the response value at any time point of interest could be back-calculated.

With regard to receptor theory, the transduction rate constant can be viewed as a kinetic extension of classical measurements of ligand efficacy. These measurements are equal to or a function of the maximal

AN US

response generated by the agonist. The transduction rate constant is the maximal rate of response generation by the agonist in the system. In classical measurements, efficacy is a function of the quantity of the response system. This is implicit in systems without receptor reserve, the maximum system potential being the Emax of a

M

full agonist, and explicit in the operational model, the quantity described as the “Maximum effect in a system” and denoted Em (Black and Leff, 1983). Likewise, kτ is a function of the quantity of response system, the

ED

transduction potential (EP). Specifically, the total transduction potential, EP(TOT), is one of the terms in the massaction equation defining kτ (kτ = [R]TOTEP(TOT)kE). In the operational model, relative agonist efficacy within a

PT

system is determined by a parameter quantifying the interaction of the receptor with the response system (Black

CE

and Leff, 1983). The parameter is the amount of agonist-occupied receptor required to elicit a response half that of Em, denoted KE. The lower the value of KE, the less RA required to generate a given level of response, and so

AC

the higher the efficacy of the agonist. In the kinetic model also, efficacy is determined by receptor-system interaction, governed by the response generation rate constant, kE. In this case, the higher the value of kE, the faster the generation of response by RA and so the higher the efficacy of the agonist. In both the operational model and the kinetic model, the receptor concentration is a determinant of agonist efficacy. Decreasing the receptor concentration decreases the level of response in the operational model and decreases the rate of response in the kinetic model. Finally, in both models an efficacy parameter emerges from the equations that 43

ACCEPTED MANUSCRIPT can be used to determine relative efficacy of a series of agonists within the same assay system (τ in the operational model, kτ in the kinetic model). In both cases the parameter can be used to distinguish efficacy of full agonists, i.e. those that produce the same maximal response at steady-state. It emerges that kτ is directly proportional to τ. Specifically, kτ is an extension of τ, in which τ is multiplied by two of the kinetic model constants: kτ = τkDEP(TOT) (Appendix B.5.1.). Affinity is defined by a single parameter, KA. (In the models with slow agonist-receptor equilibration, KA

CR IP T

is obtained by dividing k2 by k1.) It is important to note that, for GPCRs, agonist affinity in functional assays is likely system-dependent (Colquhoun, 1998). Multiple GPCR conformations are present in response systems, for example inactive and active conformations (Kenakin, 2002). The measured affinity value is a function of the microscopic affinity for particular conformations, and the relative amounts of the conformational states. Since

AN US

the mass distribution of conformational states can be dependent on system parameters, for example the concentration of G-protein, agonist affinity is a system-dependent parameter for GPCRs. Consequently, the value of KA from the kinetic model can be used to determine relative affinity values of a collection of agonists

M

when tested in the same assay system but cannot be used to compare agonist affinity values from different systems.

ED

The models can be used to measure the kinetics of receptor-ligand interaction. Binding assays have been the assay of choice for this purpose, in part because an equation is available for analyzing binding data [the

PT

Motulksy and Mahan kinetics of competitive binding equation (Motulsky and Mahan, 1984)]. Here, analogous

CE

equations were derived that can be applied to functional data to obtain estimates of the association and dissociation rate constant of test ligands. For antagonist ligands, equations were derived for estimating the

AC

binding rate constants when the agonist equilibrates rapidly with the receptor and the agonist response conforms to models without receptor reserve (Models 1 or 2). Consequently, application of the method requires selection and/or design of the assay to satisfy these criteria, for example selection of a rapidly-equilibrating agonist that conforms to Model 2. A surprisingly efficient method for quantifying antagonist binding emerges when using a continuous-read assay format, such as a biosensor. All that is required is a dilution series of antagonist and a single concentration of agonist, read at multiple time points. Notably, antagonist affinity can be determined 44

ACCEPTED MANUSCRIPT from these data sets (as the ratio of the dissociation to the association rate constant values). When using a continuous-read assay, this method for measuring antagonist binding is less resource intensive than the current standard method, Schild analysis (Arunlakshana and Schild, 1959). The kinetic method requires only a dilution series of antagonist and a single concentration of agonist. Schild analysis requires a matrix of a dilution series of both agonist and of antagonist. Agonist binding kinetics can also be assessed using the model. Binding kinetics of the agonist is

CR IP T

specifically incorporated into Models 5,6,7 and 8 (which are extensions of the response kinetic mechanisms in Models 1,2,3 and 4, respectively). From a qualitative perspective, whether or not agonist equilibrates rapidly with the receptor, relative to the time course of the response, can be determined using the method in Section 6. From a quantitative perspective, equations are derived for measuring the receptor residence time of the agonist

AN US

(Section 8). For certain applications it might be desirable to minimize the effect of slow agonist equilibration. For example, when measuring antagonist binding kinetics. This can be done by selecting and/or designing the assay appropriately. The slower the response, the less sensitive the assay will be to slow agonist equilibration. Responses with rapid onset can be substituted with slower responses in the same pathway. Examples include

M

substituting Ca2+ mobilization with inositol phosphate accumulation, or G-protein activation with gene

ED

transcription. Alternatively, the response can be triggered following a period of agonist equilibration. In a labeled GTPγS binding assay, agonist and receptor can be incubated together before addition of the labeled

PT

GTPγS. In a cAMP assay, following a period of agonist equilibration an inhibitor of cAMP metabolism can be

CE

added (e.g. 3-isobutyl-1-methylxanthine). Finally, slow equilibration can be avoided by simply increasing the concentration of agonist. This method is convenient for measuring agonist efficacy, since this measurement is

AC

routinely made using maximally-stimulating concentrations of agonist. These methods for measuring receptor binding kinetics in functional assays provide alternatives or

complements to established methods, such as radioligand binding (Aranyi, 1980; Aranyi, 1982; Bot et al., 2017; Dowling and Charlton, 2006; Fleck et al., 2012; Guo et al., 2012; Klein Herenbrink et al., 2016; Louvel et al., 2014; Motulsky and Mahan, 1984; Ramsey et al., 2011; Riddy et al., 2015; Rosethorne et al., 2016; Sykes et al., 2009; Sykes et al., 2014a; Sykes et al., 2012; Waelbroeck et al., 1991; Wittmann et al., 2011). An existing 45

ACCEPTED MANUSCRIPT functional assay method is a variant of Schild analysis, called the hemi-equilibrium method, in which antagonist dissociation rate is determined from the reduction of agonist maximal response (Kenakin et al., 2006; Mould et al., 2014; Paton and Rang, 1966; Riddy et al., 2015). In summary, the kinetic models developed in this study should aid the quantification of drug activity

AC

CE

PT

ED

M

AN US

CR IP T

necessary to maximize the potential of kinetics for GPCR research and therapeutic development.

46

ACCEPTED MANUSCRIPT Acknowledgments The authors would like to thank Scott Struthers, Philip Strange and Simeon Ramsey for invaluable insights and advise on experimental application and historical perspective. The authors are grateful for the guidance of the

AC

CE

PT

ED

M

AN US

CR IP T

anonymous reviewers.

47

ACCEPTED MANUSCRIPT Table 1. Parameters used in the kinetic models and equations Term

Description

Units

Model parameters

Receptor

Receptor unit

[R]TOT

Total concentration of receptor

Receptor unit

E

Response level

Response unit

EP

Transduction potential

Response unit

EP(TOT)

Total transduction potential

Response unit

kE

Response generation rate constant

kD

Response decay rate constant

A

Agonist

KA

Agonist equilibrium dissociation constant

k1

Agonist association rate constant

k2

Agonist dissociation rate constant

CR IP T

R

Receptor unit-1min-1 Min-1

Macroscopic pharmacological parameters

AN US

Molar Molar

Molar-1min-1 Min-1

Transduction rate constant, EP(TOT)[R]TOTkE

Response unit.min-1

τ

Transducer ratio, kτ / kDEP(TOT)

Unitless

ρA

Fractional receptor occupancy by agonist, [A] / ([A] + KA)

M



Unitless

ED

Parameters for antagonist binding kinetic equations B

Antagonist

k3

Antagonist association rate constant

k4

Antagonist dissociation rate constant

Min-1

KB

Antagonist equilibrium dissociation constant

Molar

Apparent antagonist association rate constant, [B]k3(1 – ρA) + k4

Min-1

PT

Fractional receptor occupancy by antagonist at response initiation Unitless

AC

ρB(t=0)

Molar-1min-1

CE

kB(app)

Molar

E[B]=0, t∞

Response to agonist in absence of antagonist at steady-state

Response unit

G[B]=0

Gradient of response in absence of antagonist (Model 1)

Response unit.min-1

M

Macroscopic fitting parameter, kB(app)kD

Min-2

N

Macroscopic fitting parameter, k4(kB(app) – kD)

Min-2

P

Macroscopic fitting parameter, kB(app)kD(1 – ρB(t=0))

Min-2

Q

Macroscopic fitting parameter, ρB(t=0)kD / (k4 – kD)

Unitless

48

ACCEPTED MANUSCRIPT

Table 2. Model 2 fit values for agonist-stimulated cAMP accumulation via the β2 adrenoceptor expressed in HEK293 cells. cAMP was measured as described in the legend to Fig. 5 and the data fit to Eq. 2 (Section 3.2.). Data are mean  range (n=2).

kτ (nM cAMP.min-1)

-log KA

KA

kD

(nM)

(min-1)

CR IP T

Ligand

7.0  2.7

6.78  0.05

170

0.30  0.12

Salbutamol

5.1  1.0

7.10  0.08

82

0.17  0.01

AC

CE

PT

ED

M

AN US

Ritodrine

49

ACCEPTED MANUSCRIPT

Table 3. Agonist-stimulated [35S]-GTPS binding via the -opioid receptor in C6 membranes (Traynor et al., 2002), an example of a horizontal exponential E vs t profile in which the t1/2 is dependent on agonist concentration (Model 3, EP depletion, E recycled to EP). Comparison of experimental data, and simulated data from Model 3. Experimental data are from (Traynor et al.,

CR IP T

2002). Data were simulated using Eq. 3 (Section 3.3.) with the following parameter values: KA, 130 nM; kτ, 0.012 response units.min-1; kD, 0.028 min-1; EP(TOT), 0.35 pmol/mg. The calculated value of τ, calculated using the equation at the end of Appendix B.5.1., was 1.2. 1Termed Bmax in

AN US

(Traynor et al., 2002).

[35S]GTPS binding at steady-state1

DAMGO

kobs (min-1)

concentration (nM) 0.063

50

0.039

20

0.031

Experiment

Simulation

0.057

0.16

0.18

0.037

0.095

0.089

0.032

0.061

0.049

AC

CE

PT

ED

1000

Simulation

M

Experiment

(pmol/mg)

50

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

AN US

CR IP T

Figure legends

51

AC

CE

PT

ED

M

AN US

CR IP T

ACCEPTED MANUSCRIPT

52

ACCEPTED MANUSCRIPT

Fig. 1. Time course data profile of Model 1. In this model, the E vs t profile is a straight line. The profile arises from the kinetic model when transduction potential (EP) is not depleted and response (E) does not decay. This is represented by the mechanism scheme in (a). This

CR IP T

mechanism results in a linear E vs t profile shown in (b). As the agonist concentration increases, the slope of the line increases, approaching a limit at maximally-effective agonist concentrations. In (c) the data are replotted as agonist concentration-response curves at multiple time points. When viewed in this manner, Model 1 predicts sigmoidal concentration-response curves at all

AN US

agonist concentrations. The maximal response to agonist (the upper plateau of the sigmoid curve) increases without limit as the agonist concentration increases. The potency of the agonist does not change over time, evident when plotting the agonist EC50 against time (d). Data were

AC

CE

PT

ED

M

simulated using Eq. 1, with the following parameters: kτ, 21 response units.min-1; KA, 100 nM.

53

AC

CE

PT

ED

M

AN US

CR IP T

ACCEPTED MANUSCRIPT

Fig. 2. Curve-fitting method for Model 1, linear E vs t profile. Two methods are shown for estimating kτ and KA. In the first, shown in (a) the E vs t data are fit to Eq. 1. A global fit is used, fitting E vs t data for multiple concentrations of agonist simultaneously. A Prism 7.0 file is 54

ACCEPTED MANUSCRIPT

provided in Supplementary information that can be used as a template to analyze new data (“Fig 2 Model 1 global fit”). The data shown were simulated according to Eq. 1, with a kτ value of 21 response units.min-1 and KA value of 100 nM. The fitted values are: kτ, 21.0 response units.min-1; KA, 98.0 nM. Data points are mean  standard error of duplicate determinations, generated with

CR IP T

Gaussian random error with a standard deviation of 10% (Motulsky, 2018a; Motulsky, 2018b). In the second method, the E vs t data for each agonist concentration are fit individually to a straight line equation (Response = slope  t). The fitted value of the slope is then plotted against the agonist concentration as shown in (b). These data are then fit to the sigmoid curve equation in

AN US

Section 3.1. kτ is the slope at maximally-effective agonist concentrations (“Top”) and KA the agonist concentration at the midpoint of the curve (“A50”). The fitted values are kτ, 20.9 response units.min-1; KA, 97.0 nM. The procedure is shown in the Supplementary information file, “Fig 2

M

Model 1 individual straight line fits.” (Note, for the sake of clarity, the fits to the straight line equation are not shown in the E vs t plot in (a); the fitted lines almost wholly overlap those

AC

CE

PT

ED

shown for the global fit to Eq. 1.)

55

AC

CE

PT

ED

M

AN US

CR IP T

ACCEPTED MANUSCRIPT

56

ACCEPTED MANUSCRIPT

Fig. 3. Time course data profile for Model 2. In this model, the E vs t profile is a horizontal exponential curve for which the t1/2 is independent of agonist concentration. The profile arises from the kinetic model when transduction potential (EP) is not depleted and response (E) decays. This is represented by the mechanism scheme in (a). This mechanism results in the horizontal

CR IP T

exponential E vs t profile shown in (b). This plot shows the t1/2 is independent of the agonist concentration. As the agonist concentration increases, the response at the plateau (as t  ∞) increases, approaching a limit at maximally-effective agonist concentrations. In (c) the data are replotted as agonist concentration-response curves at multiple time points. When viewed in this

AN US

manner, Model 2 predicts sigmoidal concentration-response curves at all agonist concentrations. The maximal response to agonist (the upper plateau of the sigmoid curve) increases and then approaches a limit as the agonist concentration increases. The potency of the agonist does not

M

change with time, evident when plotting the agonist EC50 against time (d). Data were simulated using Eq. 2, with the following parameters: kτ, 120 response units.min-1; KA, 100 nM; kD, 0.3

AC

CE

PT

ED

min-1.

57

AC

CE

PT

ED

M

AN US

CR IP T

ACCEPTED MANUSCRIPT

58

ACCEPTED MANUSCRIPT

Fig. 4. Curve-fitting method for Model 2, horizontal exponential E vs t profile, t1/2 independent of agonist concentration. Two methods are shown for estimating kτ and KA. In the first, shown in (a), the E vs t data are fit to Eq. 2. A global fit is used, fitting E vs t data for multiple concentrations of agonist simultaneously. A Prism 7.0 file is provided in Supplementary

CR IP T

information that can be used as a template to analyze new data (“Fig 4 Model 2 global fit”). The data shown were simulated according to Eq. 2, with a kτ value of 120 response units.min-1, a KA value of 100 nM and kD value of 0.3 min-1. The fitted values are: kτ, 121 response units.min-1;

AN US

KA, 99.9 nM; kD, 0.303 min-1. Data points are mean  standard error of duplicate determinations, generated with Gaussian random error with a standard deviation of 10% (Motulsky, 2018a; Motulsky, 2018b). In the second method, the E vs t data for each agonist concentration are fit individually to a horizontal exponential equation, Response = Plateau  (1 – e-K.t). The fitted

M

value of the plateau is then multiplied by the fitted value of K, the rate constant. The resulting value is then plotted against the agonist concentration, shown in (b). These data are then fit to the

ED

sigmoid curve equation in Section 3.1. kτ is the Plateau  K value at maximally-effective agonist concentrations (“Top”) and KA the agonist concentration at the midpoint of the curve (“A50”).

PT

The fitted values are kτ, 120 response units.min-1; KA, 95.7 nM. The procedure is shown in the

CE

Supplementary information file, “Fig 4 Model 2 individual straight line fits.” (Note, for the sake of clarity the fits to the horizontal exponential equation are not shown in the E vs t plot (a); the

AC

fitted lines almost wholly overlap those shown for the global fit to Eq. 2.)

59

AC

CE

PT

ED

M

AN US

CR IP T

ACCEPTED MANUSCRIPT

60

ACCEPTED MANUSCRIPT

Fig. 5. Time course data profile for Model 3. In this model, the E vs t profile is a horizontal exponential curve for which the t1/2 is dependent on the agonist concentration. The profile arises from the kinetic model when transduction potential (EP) is depleted and response (E) decays and recycles back to EP. This is represented by the mechanism scheme in (a). This mechanism results

CR IP T

in the horizontal exponential E vs t profile shown in (b). This plot shows the t1/2 is dependent on the agonist concentration. The t1/2 decreases as the agonist concentration increases, approaching a limit at maximally-effective agonist concentrations. In addition, as the agonist concentration increases, the response at the plateau (as t  ∞) increases than approaches a limit at maximally-

AN US

effective agonist concentrations. In (c) the data are replotted as agonist concentration-response curves at multiple time points. When viewed in this manner, Model 3 predicts sigmoidal concentration-response curves at all agonist concentrations. The maximal response to agonist

M

(the upper plateau of the sigmoid curve) increases and then approaches a limit as the agonist concentration increases. The potency of the agonist increases over time. This increase is evident

ED

when plotting the agonist EC50 against time (d). Data were simulated using Eq. 3, with the

AC

CE

response units.

PT

following parameters: kτ, 120 response units.min-1; KA, 100 nM; kD, 0.3 min-1; EP(TOT), 500

61

AC

CE

PT

ED

M

AN US

CR IP T

ACCEPTED MANUSCRIPT

62

ACCEPTED MANUSCRIPT

Fig. 6. Curve-fitting method for Model 3, horizontal exponential E vs t profile, t1/2 dependent on agonist concentration. Two methods are shown for estimating kτ and KA. In the first, shown in (a) the E vs t data are fit to Eq. 3. A global fit is used, fitting E vs t data for multiple concentrations of agonist simultaneously. A Prism 7.0 file is provided in Supplementary information that can be

CR IP T

used as a template to analyze new data (“Fig 6 Model 3 global fit”). The data shown were simulated according to Eq. 3, with the following parameter values: kτ, 750 response units.min-1; KA, 100 nM; kD, 0.3 min-1; EP(TOT), 500 response units. The fitted values are: kτ, 750 response

standard error

AN US

units.min-1; KA, 95.2 nM; kD, 0.315 min-1; EP(TOT), 504 response units. Data points are mean  of duplicate determinations, generated with Gaussian random error with a

standard deviation of 10% (Motulsky, 2018a; Motulsky, 2018b). In the second method, the E vs t data for each agonist concentration are fit individually to a horizontal exponential equation,

M

Response = Plateau  (1 – e-K.t). The fitted value of the plateau is then multiplied by the fitted value of K, the rate constant. The resulting value is then plotted against the agonist concentration

ED

as shown in (b). These data are then fit to the sigmoid curve equation in Section 3.1. kτ is the

PT

Plateau  Rate constant value at maximally-effective agonist concentrations (“Top”) and KA the agonist concentration at the midpoint of the curve (“A50”). The fitted values are kτ, 767 response

CE

units.min-1; KA, 102 nM. The procedure is shown in the Supplementary information file, “Fig 6 Model 3 individual straight line fits.” (Note, for the sake of clarity, the fits to the horizontal

AC

exponential equation are not shown in the E vs t plot; the fitted lines almost wholly overlap those shown for the global fit to Eq. 3.)

63

AC

CE

PT

ED

M

AN US

CR IP T

ACCEPTED MANUSCRIPT

64

ACCEPTED MANUSCRIPT

Fig. 7. Time course data profile for Model 4. In this model, the E vs t profile is a rise-and-fall exponential curve. The profile arises from the kinetic model when transduction potential (EP) is depleted and response (E) decays. This is represented by the mechanism scheme in (a). This mechanism results in the rise-and-fall exponential E vs t profile shown in (b). The peak response

CR IP T

increases as agonist concentration increases, approaching a limit at maximally-effective agonist concentrations. The steepness of the rise phase also increases with increasing agonist concentration. The response declines to zero for all concentrations of agonist. In (c) the data are replotted as agonist concentration-response curves at multiple time points. When viewed in this

AN US

manner, the curve is sigmoid at early time points (0.3 to 4 min) but bell-shaped at later time points (10 and 20 min). In the early time point range, the potency of the agonist increases over time. This increase is evident when plotting the agonist EC50 against time (d). Data were

M

simulated using Eq. 4, with the following parameters: kτ, 1050 response units.min-1; KA, 100 nM;

AC

CE

PT

ED

kD, 0.3 min-1; EP(TOT), 700 response units.

65

AC

CE

PT

ED

M

AN US

CR IP T

ACCEPTED MANUSCRIPT

66

ACCEPTED MANUSCRIPT

Fig. 8. Curve-fitting method for Model 4, rise-and-fall E vs t profile. Two methods are shown for estimating kτ and KA. In the first, shown in (a), the E vs t data are fit to Eq. 4. A global fit is used, fitting E vs t data for multiple concentrations of agonist simultaneously. A Prism 7.0 file is provided in Supplementary information that can be used as a template to analyze new data (“Fig

CR IP T

8 Model 4 global fit”). The data shown were simulated according to Eq. 4, with the following parameter values: kτ, 1050 response units.min-1; KA, 100 nM; kD, 0.3 min-1; EP(TOT), 700 response units. The fitted values are: kτ, 1054 response units.min-1; KA, 102 nM; kD, 0.300 min-1; EP(TOT),

AN US

699 response units. Data points are mean  standard error of duplicate determinations, generated with Gaussian random error with a standard deviation of 10% (Motulsky, 2018a; Motulsky, 2018b). In the second method, the E vs t data for each agonist concentration are individually fit to a rise-and-fall exponential equation, Response = C / (K1 – K2)  (e-K2.t – e-K1.t). The fitted value

M

of C is then plotted against the agonist concentration as shown in (b). These data are then fit to the sigmoid curve equation in Section 3.1. kτ is the value of C at maximally-effective agonist

ED

concentrations (“Top”) and KA the agonist concentration at the midpoint of the curve (“A50”). The fitted values are kτ, 1050 response units.min-1; KA, 98.1 nM. The procedure is shown in the

PT

Supplementary information file, “Fig 8 Model 4 individual straight line fits.” (Note, for the sake

CE

of clarity, the fits to the rise-and-fall equation are not shown in the E vs t plot (a); the fitted lines

AC

almost wholly overlap those shown for the global fit to Eq. 4.)

67

M

AN US

CR IP T

ACCEPTED MANUSCRIPT

Fig. 9. Glucagon-stimulated cAMP accumulation in rat hepatic membranes, an example of

ED

Model 1, no EP depletion, no E decay, a linear E vs t profile (Section 3.1.). Partially purified plasma membranes, largely devoid of phosphodiesterase activity, were incubated with a range of

PT

concentrations of glucagon and cAMP measured. Depletion of substrate, ATP, was minimized

CE

with the use of a cAMP-regenerating system. Data are reproduced from (Rodbell et al., 1974) with permission. The dashed red line is the fit of the data to Model 1 (Eq. 1 in Section 3.1.). The

AC

fitted parameter values are: kτ, 22 response units.min-1; KA, 0.76 nM. The fit was performed using Prism 7.0 (GraphPad Software, Inc.). Data were fit to the black straight lines of the original figure, which were reconstructed using a plot digitizer program [WebPlotDigitizer (Rohatgi, 2018)].

68

AC

CE

PT

ED

M

AN US

CR IP T

ACCEPTED MANUSCRIPT

69

ACCEPTED MANUSCRIPT

Fig. 10. Agonist-stimulated cAMP accumulation in HEK293 cells expressing the β2 adrenoceptor, an example of Model 2 (no EP depletion, E decay, Section 3.2.) (a) Ritodrine. (b) salbutamol. cAMP was measured using an HTRF-based competitive immunoassay (Cisbio, Bedford, MA). 5 l of assay buffer containing ritodrine or salbutamol, together with 3-isobutyl-

CR IP T

1-methylxanthine (final assay concentration of 1mM), was added to a 384-well plate. 5 l Taglite HEK293 cells stably expressing the β2 adrenoceptor (catalogue number C1SU1BETA2) were then added using the injector of a Synergy Neo2 Hybrid Multi-Mode Reader (BioTek, Winooski,

AN US

VT). Following the indicated incubation times, 6 l lysis/detection buffer were injected into the wells to stop the production of cAMP. Following the final time point, 5 l cAMP detection reagent was added (d2-labeled cAMP and Eu3+-labeled anti-cAMP antibody), the plates incubated at room temperature for 1 hr, then the HTRF signal measured in the Synergy Neo2

M

reader. Data were converted to nM cAMP from a cAMP standard curve. Data were fit to an extended version of Eq. 2 (Table 2); a constant representing the baseline was added to the

ED

equation. Data were fit globally, with time and agonist concentration as the independent

PT

variables, using Prism 7.0 (GraphPad Software, Inc.). Data points are the mean  standard error of duplicate measurements. Data are from representative experiments performed twice. The

AC

CE

averages of the fitted parameter values are given in Table 2.

70

PT

ED

M

AN US

CR IP T

ACCEPTED MANUSCRIPT

CE

Fig. 11. Bombesin-stimulated Ca2+ mobilization via the BB1 bombesin receptor expressed in COS-7 cells. Data are from (Liu et al., 2003), reconstructed using a plot digitizer program

AC

[WebPlotDigitizer (Rohatgi, 2018)]. The line is the fit to the form Eq. 4 that assumes a maximally-stimulating concentration of agonist (Section 3.4.). The fitted values were, kτ, 1,500 FIU.sec-1; kD, 0.042 sec-1; EP(TOT), 5,200 FIU. The calculated value of of τ, calculated using the equation at the end of Appendix B.5.1., was 6.7.

71

AC

CE

PT

ED

M

AN US

CR IP T

ACCEPTED MANUSCRIPT

Fig. 12. Ligand addition schedules for measuring antagonist binding kinetics in a functional assay. The order of antagonist and agonist addition to the system can be assay-specific and the antagonist binding kinetic equations are formulated to accommodate four schedules. (a)

72

ACCEPTED MANUSCRIPT

Simultaneous addition of agonist and antagonist (Eq. 5). (b) Antagonist pre-incubation (Eq. 6).

AC

CE

PT

ED

M

AN US

CR IP T

(c) Antagonist pre-incubation and washout (Eq. 7). (d) Agonist pre-incubation (Eq. 8).

Fig. 13. Measuring antagonist binding kinetics in a functional assay in which receptor is exposed simultaneously to antagonist and agonist (Section 5.1.). Simulated E vs t data are shown for inhibition of a functional response to agonist by multiple concentrations of antagonist. The response mechanism is Model 2 (Section 3.2.). The time course data for all conditions (agonist alone, and agonist with three concentrations of antagonist) were analyzed simultaneously by 73

ACCEPTED MANUSCRIPT

global fitting to Eq. 5 (Section 5.1.) using Prism 7.0. The fitted values are; k3, 1.01  107 M-1min1

; k4 0.0234 min-1; kD, 0.101 min-1; E[B]=0,t∞, 399 response units. A Prism 7.0 file detailing the

analysis is provided in Supplementary information (“Fig 12 Antagonist kinetics, simultaneous addition”). In this file, Eq. 5 has been entered as a “User-defined equation” titled “Antagonist vs

CR IP T

agonist response kinetics, simultaneous addition.”. Data points are mean  standard error of duplicate determinations, generated with Gaussian random error with a standard deviation of 10% (Motulsky, 2018a; Motulsky, 2018b). The parameter values used to simulate the data were:

AN US

k3, 1.0  107 M-1min-1; k4 0.0231 min-1; kD, 0.10 min-1; E[B]=0,t∞, 400 response units; [A], 300

AC

CE

PT

ED

M

nM; KA, 100 nM

74

AC

CE

PT

ED

M

AN US

CR IP T

ACCEPTED MANUSCRIPT

Fig. 14. Measuring antagonist binding kinetics in a functional assay in which receptor is preincubated with antagonist (Section 5.2.). Simulated E vs t data are shown for inhibition of a functional response to agonist by multiple concentrations of antagonist. The response mechanism 75

ACCEPTED MANUSCRIPT

is Model 2 (Section 3.2.). Antagonist is pre-incubated with receptor for 30 minutes prior to the addition of agonist and measurement of the response. The time course data for all conditions (agonist alone, and agonist with three concentrations of antagonist) were analyzed simultaneously by global fitting to Eq. 6 (Section 5.2.) using Prism 7.0. The fitted values are; k3,

CR IP T

0.971  107 M-1min-1; k4 0.0226 min-1; kD, 0.101 min-1; E[B]=0,t∞, 401 response units. A Prism 7.0 file detailing the analysis is provided in Supplementary information (“Fig 13 Antagonist kinetics, antagonist pre-incubation”). In this file, Eq. 6 has been entered as the “User-defined equation” titled “Antagonist vs agonist response kinetics, antagonist pre-incubation.” Data points

AN US

are mean  standard error of duplicate determinations, generated with Gaussian random error with a standard deviation of 10% (Motulsky, 2018a; Motulsky, 2018b). The parameter values used to simulate the data were: k3, 1.0  107 M-1min-1; k4 0.0231 min-1; kD, 0.10 min-1; E[B]=0,t∞,

AC

CE

PT

ED

M

400 response units; [A], 300 nM; KA, 100 nM; antagonist pre-incubation time, 30 min.

76

AC

CE

PT

ED

M

AN US

CR IP T

ACCEPTED MANUSCRIPT

Fig. 15. Measuring antagonist binding kinetics in a functional assay in which receptor is preincubated with antagonist then washed out. The response mechanism is Model 2 (Section 3.2.). 77

ACCEPTED MANUSCRIPT

This procedure, applicable to ex vivo receptor occupancy assays, is detailed in Section 5.3. It is the functional analogue of a ligand binding method described previously (Malany et al., 2009). Simulated E vs t data are shown for inhibition of a functional response to agonist by antagonist. It is assumed an antagonist is administered to animals at the doses shown, then the relevant tissue

CR IP T

removed and free antagonist washed out. The response to agonist in the washed preparation is then measured. The experiment includes a „no antagonist‟ control, i.e. animals dosed with vehicle only. The data were analyzed as follows. First, the data for the vehicle control were fit to the equation,

AN US

Et = E[B]=0,t∞(1 – e–kD.t) in Section 5.3. to determine the value of kD. The fitted values are, kD, 0.103 min-1; E[B]=0,t∞, 400 response units. Next, the time course data for all doses were analyzed simultaneously by global fitting to Eq. 7 (Section 5.3.) using Prism 7.0. The fitted values are; k4

M

0.0468 min-1; ρB(t=0), 13.4, 35.7, 67.2, 87.2 and 95.8% for 0.3, 1, 3, 10 and 30 mg/kg, respectively; E[B]=0,t∞, 399, 397, 397, 400 and 398 response units for 0.3, 1, 3, 10 and 30 mg/kg,

ED

respectively. kD was entered as a constant, the value from the control experiment (0.103 min-1). A Prism 7.0 file detailing the analysis is provided in Supplementary information (“Fig 14

PT

Antagonist kinetics, antagonist pre-incubation and wash”. In this file, Eq. 7 has been entered as

CE

the “User-defined equation” titled “Antagonist vs agonist response kinetics, antagonist pre-inc. & wash.” Data points are mean  standard error of duplicate determinations, generated with

AC

Gaussian random error with a standard deviation of 10% (Motulsky, 2018a; Motulsky, 2018b). The parameter values used to simulate the data were: kD, 0.10 min-1; k4 0.0462 min-1; ρB(t=0), 15, 35, 68, 87 and 96% for 0.3, 1, 3, 10 and 30 mg/kg, respectively; E[B]=0,t∞, 400 response units for all samples; [A], 300 nM; KA, 100 nM. Inset: occupancy vs dose data. The data were analyzed using a four parameter-logistic equation, giving an ED50 value of 1.57 mg/kg. 78

AC

CE

PT

ED

M

AN US

CR IP T

ACCEPTED MANUSCRIPT

Fig. 16. Measuring antagonist binding kinetics in a functional assay in which receptor is preincubated with agonist (Section 5.4.). The response mechanism is Model 2 (Section 3.2.). 79

ACCEPTED MANUSCRIPT

Simulated E vs t data are shown for inhibition of a functional response to agonist by multiple concentrations of antagonist. Agonist is pre-incubated with receptor until the plateau is reached. Antagonist is then added and response continues to be measured. The time course data for all conditions (agonist alone, and agonist with six concentrations of antagonist) were analyzed

CR IP T

simultaneously by global fitting to Eq. 8 (Section 5.4.) using Prism 7.0. The fitted values are; k3, 0.973  107 M-1min-1; k4 0.0223 min-1; kD, 0.0996 min-1; E[B]=0,t∞, 398 response units. A Prism 7.0 file detailing the analysis is provided in Supplementary information (“Fig 15 Antagonist kinetics, agonist pre-incubation”). In this file, Eq. 8 has been entered as the “User-defined

AN US

equation” titled “Antagonist vs agonist response kinetics, agonist pre-incubation.” Data points are mean  standard error of duplicate determinations, generated with Gaussian random error with a standard deviation of 10% (Motulsky, 2018a; Motulsky, 2018b). The parameter values

M

used to simulate the data were: k3, 1.0  107 M-1min-1; k4 0.0231 min-1; kD, 0.10 min-1; E[B]=0,t∞,

AC

CE

PT

ED

400 response units; [A], 300 nM; KA, 100 nM; agonist pre-incubation time, 30 min.

80

AC

CE

PT

ED

M

AN US

CR IP T

ACCEPTED MANUSCRIPT

81

ACCEPTED MANUSCRIPT

Fig. 17. Method for testing rapid vs slow agonist-receptor binding kinetics. Agonist is applied for a period of time during which response is generated (the stimulation phase). A large excess of a competitive antagonist is then added, at the time point indicated by the arrow. The response continues to be recorded (the antagonist block phase). The solid line following antagonist

CR IP T

addition is the response for an agonist with rapid receptor binding kinetics, and the dashed line the response for an agonist with slow binding kinetics. (a) Model 1, no EP depletion, no E decay; (b) Model 2, no EP depletion, E decay; (c) Model 3, EP depletion, E recycles to EP; (d) EP depletion, E decay. Data for the antagonist block phase were simulated using the equations in

AN US

Appendix E.5. The values of the model parameters kτ, kD, KA and EP(TOT) for Models 1, 2, 3 and 4 are the same as in Fig. 1, 3, 5 and 7, respectively. The concentration of agonist is 300 nM. For the slow binding kinetics curves, the agonist dissociation rate constant (k2) is 0.06 min-1 (t1/2 of

M

12 min). Antagonist was added at 10 min (a-c) and 0.4 min (d). Solid black symbol – response value that is 50% of the response at the time of antagonist addition. For the rapidly-equilibrating

AC

CE

PT

ED

models this level is reached at a time after antagonist addition equal to the kD t1/2 (2.3 min).

82

PT

ED

M

AN US

CR IP T

ACCEPTED MANUSCRIPT

CE

Fig. 18. E vs t profiles for kinetic models of agonist-stimulated response with rate-limiting

AC

agonist-receptor binding kinetics. (a-c) Model 5: No EP depletion, no E decay, slow agonist equilibration. (d-f) Model 6: No EP depletion, E decay, slow agonist equilibration. (g-i) Model 7: EP depletion, E recycles to EP, slow agonist equilibration. (j-l) Model 4: EP depletion, E decay, slow agonist equilibration. The upper graphs are the E vs t data, the middle graphs the agonist dose response data at multiple time points, and the lower graphs the agonist EC50 plotted against time. (The dashed line in the time course graphs (a,d,g, and j) is the curve for the corresponding 83

ACCEPTED MANUSCRIPT

model with rapid agonist equilibration (Models 1-4), at 22,000 nM, a maximally-stimulating concentration.) Data were simulated using Eqs. 9-12, with agonist association and dissociation rate constants of 3  105 M-1min-1 and 0.03 min-1, respectively. The resulting residence time is 23 min. Agonist affinity (KA) is 100 nM. The remaining parameters are those used to simulate the

CR IP T

data for Models 1-4 (Fig. 1,3,5 and 7), so that the profiles of rapid and slow agonist equilibration can be compared. kD is 0.3 min–1 for all models. The kτ values are, in units of response units.min1

: Model 5, 21; Model 6, 120; Model 7, 750; Model 8, 1,050. The EP(TOT) values are, in units of

response units: Model 7, 500; Model 8; 700. The dose-response data in panels b, e, h and k were

AN US

fitted to a four-parameter logistic equation using Prism 7.0. In the lower graphs, the fitted value of EC50 is plotted against the time point (c,f,i and l). (Note at the later time points for Model 8

AC

CE

PT

ED

M

the dose-response was bell shaped so no EC50 value was obtained.)

84

ACCEPTED MANUSCRIPT

Fig. 19. Method for measuring agonist residence time from reduction of response following antagonist addition. The response mechanism is Model 7 (No EP depletion, E decline, slow agonist equilibration). The lower panel shows the sequence and timing of ligand addition. Receptor is first exposed to agonist for a period long enough for a stable response (the plateau) to

CR IP T

be reached. A large excess of antagonist is then added, sufficient to prevent any further association of agonist with receptor (t = 0). After this point, response declines as agonist dissociates from the receptor. The data in the graph were simulated using Eq. 13 (Appendix E.6.3.), with a E(t=0) value of 200 response units and a kD value of 0.3 min-1 The agonist

AN US

dissociation rate constant values (k2) are 0.00149, 0.00578, 0.0231, 0.139, and 0.693 min-1, with corresponding residence times of 420, 120, 30, 5 and 1 min, respectively. The control (brown line) is the response when antagonist is omitted. The red line is for an agonist that dissociates

M

sufficiently rapidly that the agonist-occupied receptor concentration immediately falls to zero on agonist addition. Under these conditions, the rate of response decline is the same as the rate of

ED

response decay, i.e. kD (Section 6).

CE

PT

Appendix A: Mathematical methods

A.1. Formularization strategy

AC

In this study, equations were derived describing the level of response over time according to various mechanisms. The final equations take the form y = f(t), where y is response, and f(t) is a function of time and pharmacological parameters. Throughout this study, this form is described as an E vs t equation. In the first step of deriving the relevant E vs t equation, the model mechanism is formularized using differential equations. These equations are then solved

85

ACCEPTED MANUSCRIPT

analytically to obtain the E vs t equation. For simple systems, such as Models 1 and 2, the differential equations were solved directly, by integration, to obtain the E vs t solution. For more complex mechanisms, such as the models that incorporate kinetics of ligand binding, the Laplace transform method is used here to obtain the E vs t equation from the differential equations

CR IP T

(Hoare, 2017; Mayersohn and Gibaldi, 1970). This solution procedure, involving an intermediate step (the Lapace transform), is useful for mechanisms in which there is more than one timedependent term. For example, response level, E, and receptor concentration [RA]. In the Laplace transform domain, expressions for these terms can be substituted into one another algebraically.

AN US

This obviates the complexities associated with directly solving the differential equations, in particular the different algebraic relationships in the differential equations vs their analytic

M

solutions.

A.2. Verification method

ED

For the equations solved with the Laplace transform method, a procedure was used to verify the E vs t equation. (These are Eqs. 4-10 and 12-14). Data simulated with the E vs t

PT

equation were compared with data simulated by numerical solution of the differential equations.

CE

Numerical solution was performed using the Euler method, as described in (Vauquelin et al., 2001), using 1,200 data points. The dt value was the last time point in the E vs t equation

AC

simulation divided by 1,200. In all cases, the correlation coefficient (R2) comparing the analytic and numerical solutions was > 0.999. We advocate this validation procedure be used when solving more complex differential equations, because the analytical solution procedure requires many steps and the predicted E vs t profiles are often counter-intuitive. In the initial exploration of the models here it revealed mathematical errors on two occasions.

86

ACCEPTED MANUSCRIPT

A.3. Second order mechanisms In some models, E vs t equations were not obtained. Specifically, models incorporating both the kinetics of ligand binding and depletion of transduction potential (for example, Models

CR IP T

7 and 8). The ease of applying the methods above is dependent on the “reaction order” of the interactions with respect to time. The reaction order of the forward process, response generation by the receptor, is dependent on whether EP or [RA] vary with time. If only one of these terms is time-dependent, the E vs t equation can be readily derived from the differential equation. If both

AN US

[RA] and EP are time dependent, the reaction is second order with respect to time. Second order reaction processes are challenging to formularize analytically. The Laplace transform method cannot be applied to second order reaction processes (Mayersohn and Gibaldi, 1970). It might be

M

possible to obtain analytic E vs t equations for the second order mechanisms here using different mathematical techniques. In the case of second order mechanisms, data were analyzed by

PT

ED

numerical solution of the differential equations as described above (Appendix A.2.).

CE

Appendix B. Derivation of kinetic model equations, for rapid agonist-receptor equilibration

AC

In this section, the mathematical aspects of the kinetic agonism model are considered for the models with rapid agonist-receptor equilibration, and the equations are derived. The theoretical framework of the model is given in Section 2.1. The model parameters are defined in Section 2.2. and listed in Table 1. In the kinetic response model, the transduction potential, EP, is

87

ACCEPTED MANUSCRIPT

converted to the response, E, by the agonist occupied receptor. The response can decay. This

AN US

CR IP T

mechanism is illustrated by the following general scheme:

M

Equations defining the level of response over time were derived from first principles. Specifically: 1) The rate of response generation is defined by the concentration of the interacting

ED

components and a rate constant, kE. 2) The rate of response decay is defined by the exponential decay rate constant kD. 3) The rate of change of response with time is equal to the rate of

PT

response generation minus the rate of response decline. Based on these assumptions, we can

AC

CE

write the differential equation for E:

From this general model, four variants are considered that result in differently-shaped E vs t profiles. The profile is dependent on the fate of the response, specifically whether or not it decays and, if it does, whether or not it recycles back to EP. The profile is also dependent on 88

ACCEPTED MANUSCRIPT

whether or not EP is depleted by conversion to E. In order of increasing complexity, the E vs t profiles are linear; horizontal exponential, t1/2 independent of [A]; horizontal exponential, t1/2 dependent on [A]; and rise-and-fall. (The horizontal exponential profile is defined by the equation y = Plateau  (1 – e–K.t). A familiar example is the ligand binding association curve. t1/2

CR IP T

is the time at which half of the steady-state response has been reached.) These variants are considered in the proceeding sections.

B.1. Model 1: No EP depletion, no E decay

AN US

In this model, there is no depletion of EP and E does not decay. The resulting E vs t plot

M

is a straight line (Fig. 1B). The mechanism is represented by Scheme 1:

Since EP is not depleted by generation of E, EP is independent of time and is a constant in the

ED

differential equation for E. Since E does not decline, the rate constant kD is zero and the second term of the differential equation is zero. Since agonist binding is assumed to be at equilibrium,

AC

CE

PT

[RA] is constant. The resulting differential equation is,

where EP(TOT) is the total amount of transduction potential in the system. (Note, in this model and Model 2 EP(TOT) is equal to EP because EP is not depleted. In Models 3 and 4, in which EP is depleted, EP(TOT) is EP at the moment of agonist addition, which is equal to the sum of EP and E at

89

ACCEPTED MANUSCRIPT

any given time point. Using EP(TOT) allows a common terminology to be applied for all the

CR IP T

models.) Integration yields the E vs t equation,

The agonist concentration, [A], is incorporated by substituting [RA] with the equation for

AN US

equilibrium occupancy of receptor by agonist:

where KA is the agonist equilibrium dissociation constant. Substituting into the E vs t equation

ED

M

gives,

PT

Eq. 1

CE

where,

AC

kτ is the transduction rate constant, the agonist efficacy parameter in the model. The properties of this parameter are described in Section 2.3. It is evident that when t is the independent variable, Eq. 1 is an equation for a straight line (see Fig. 1B for simulated E vs t plot). Eq. 1 can be used directly to fit E vs t data to obtain estimates of the model parameters, using a global fit with two independent variables, time and agonist concentration (Section 3.1., 90

ACCEPTED MANUSCRIPT

Fig. 2A, Supplementary information file “Fig 2 Model 1 global fit”). Here an alternative method is developed, utilizing familiar empirical equations. The E vs t data for each agonist

s ope.

s ope

AN US

It is evident by comparison with Eq. 1 that the slope is equal to,

CR IP T

concentration is fit individually to a general linear equation, with y intercept of zero:

From this equation it is evident that plotting the slope against the agonist concentration gives a

M

hyperbola. It is conventional to use the logarithm of the ligand concentration in pharmacological

ED

data analysis. In addition, a slope factor is introduced to allow for physical or biological effects that deviate the range of the agonist serial dilution series in the vicinity of the receptor from that

PT

intended (see Section 3.1.). This is handled by incorporating the Hill Slope. The resulting

op

og

og

i

ope

AC

CE

sigmoid equation, used to analyze the slope versus log[A] data, is,

In this plot, y is the slope from the E vs t plot. kτ is the slope value at infinite [A], i.e. the asymptote of the sigmoid curve. KA is the antilog of the x value at the midpoint of the curve, log(A50). The application of this method to experimental data is presented in Section 3.1., Fig. 2 and the Supplementary information file, “Fig 2 Model 1 individual straight line fits.” 91

ACCEPTED MANUSCRIPT

B.2. Model 2: No EP depletion, E decay In this model, there is no depletion of EP and E decays. The resulting E vs t profile is a

CR IP T

horizontal exponential curve, in which the t1/2 of the response is independent of agonist concentration (Fig. 3B). The mechanism is represented by Scheme 2:

In this model, the value of EP is independent of time because the level of EP is unaffected by the

CE

PT

ED

Integration yields the E vs t equation,

M

AN US

generation of E (the time-dependent process). The resulting differential equation for E is,

AC

The full equation, incorporating the concentration of agonist, is:

Eq. 2

92

ACCEPTED MANUSCRIPT

It is evident that when t is the independent variable, Eq. 2 is a horizontal exponential equation (see Fig. 3B for simulated E vs t plot). The t1/2 of this E vs t equation is –ln(0.5)/kD. [A] is not present in this term, indicating t1/2 is independent of agonist concentration. Eq. 2 can be used directly to fit E vs t data to obtain estimates of the model parameters,

CR IP T

using a global fit with two independent variables, time and agonist concentration (Section 3.2., Fig. 4A, Supplementary information file, “Fig 4 Model 2 global fit”). It is also possible to analyze the data utilizing familiar empirical equations. First, the E vs t data for each agonist

AN US

concentration is fit individually to the general horizontal exponential equation,

te u

M

where K is a rate constant. It is evident by comparison with Eq. 2 that the plateau is equal to,

PT

ED

te u

This is because as t approaches infinity, the exponential term approaches zero. It is also evident

CE

that K is equal to kD. Next, the plateau value is multiplied by the rate constant value. The

AC

resulting equation is,

te u

This is the same equation as that for Model 1, except “Slope” in model 1 is replaced with Plateau  K. The data analysis now proceeds as described for Model 1 (Appendix B.1.): The Plateau  K

93

ACCEPTED MANUSCRIPT

versus log[A] data are fit to the sigmoid curve equation. kτ is the Plateau  K value at infinite [A]. KA is the antilog of the x value at the midpoint of the curve, log(A50). The application of this method to experimental data is presented in Section 3.2., Fig. 4 and the Supplementary

CR IP T

information file, “Fig 4 Model 2 individual horizontal exponential fits.”

B.3. Model 3: EP depletion, E recycles to EP

In this model, EP is depleted by generation of E, and E recycles back to EP. The resulting E vs t profile is a horizontal exponential curve, in which the t1/2 of the response is dependent on

PT

ED

The differential equation for E is,

M

AN US

the agonist concentration (Fig. 5B). The mechanism is represented by Scheme 3:

CE

In formulating the E vs t equation, we need to consider that EP no longer remains constant over

AC

time. Consequently, the differential equation for E contains two time-dependent parameters, E and EP. For the formulation to proceed, we need an equation in which the only time-dependent parameter is the one being measured, E. Mathematically, this can be handled by an equation summing the amounts of E and EP:

94

ACCEPTED MANUSCRIPT

where EP(TOT) is the total amount of transduction potential. (EP(TOT) is also equal to the transduction potential prior to the addition of agonist.) EP(TOT) is a constant in the differential

CR IP T

equation because it is independent of time. Next, we solve the sum equation, for EP:

In this equation, EP is defined by the time-dependent variable being measured, E, and the time-

AN US

independent parameter EP(TOT). Next, we substitute this expression into the differential equation

ED

M

for E, giving,

This is the desired equation – the only time-dependent parameter is E. Integrating we obtain the

CE

PT

E vs t equation:

AC

(

)

The half-time of the response is –ln(0.5)/([RA]kE + kD). Since this term contains [RA], which is dependent on [A], the half-time is dependent on the agonist concentration. On incorporating the concentration of agonist, the equation becomes:

95

ACCEPTED MANUSCRIPT

( (

(

)

)

)

CR IP T

kτ can be incorporated as follows. kτ is defined as in Model 1 (Appendix B.1.):

AN US

Solving for [R]TOTkE gives,

M

Substituting into the E vs t equation and re-arranging gives the final equation used for fitting the

⁄ ⁄

(

) Eq. 3

AC

CE

PT

ED

data:

where,

96

ACCEPTED MANUSCRIPT

Eq. 3 can be used directly to fit E vs t data to obtain estimates of the model parameters, using a global fit with two independent variables, time and agonist concentration (Section 3.3., Fig. 6A, Supplementary information file, “Fig 6 Model 3 global fit”). It is also possible to

CR IP T

analyze the data using familiar empirical equations. First, the E vs t data for each agonist concentration is fit individually to the general horizontal exponential equation,

AN US

te u

where K is a rate constant. It is evident by comparison with Eq. 3 that the plateau is equal to,

M



ED

te u



)

PT

(

AC

CE

And that the rate constant K is equal to,

Next, the plateau value is multiplied by the rate constant value. After some rearranging, the resulting equation is,

97

ACCEPTED MANUSCRIPT

te u

This is the same equation as that for Model 2. The data analysis now proceeds as described for

CR IP T

Model 2 (Appendix B.2.): The Plateau  K versus log[A] data are fit to the sigmoid curve equation. kτ is the Plateau  K value at infinite [A]. KA is the antilog of the x value at the midpoint of the curve, log(A50). The application of this method to experimental data is presented in Section 3.3. Fig. 6 and the Supplementary information file, “Fig 6 Model 3 individual

AN US

horizontal exponential fits.”

B.4. Model 4: EP depletion, E decay

In this model, there is both depletion of EP and decay of E. A rise-and-fall E vs t profile results

ED

M

(Fig. 7B). The mechanism is represented by Scheme 4:

AC

CE

PT

The differential equation for E is,

Under these conditions, EP is a time-dependent variable, as in Model 3. However, in contrast to Model 3, the equation for the sum of response terms incorporates the response that is lost owing to the counteracting process, here termed D:

98

ACCEPTED MANUSCRIPT

CR IP T

Substituting into the differential equation for E gives,

This equation contains two time-dependent terms, E and D. In order to proceed with the

AN US

formulation, we want an equation in which the only time-dependent parameter is the one being measured, E. Mathematically, this can be handled using the Laplace transform method (Hoare, 2017; Mayersohn and Gibaldi, 1970). First, the differential equation for D is written. D is

ED

M

generated from E at the rate kD:

PT

The expressions for E and D are combined using their Laplace transforms. The transforms for E

AC

CE

and D are, respectively,

̅

̅

̅

̅

̅

Solving the second equation for ̅ and substituting into the transform for E gives, 99

ACCEPTED MANUSCRIPT

̅

̅

̅

̅

AN US

Taking the inverse Laplace transform gives the E vs t equation:

CR IP T

We next solve for ̅ :

)

M

(

PT

ED

On incorporating the concentration of agonist, the equation becomes:

(



)

/

CE

.



AC

kτ can be incorporated as done in Model 3 (Appendix B.3.), giving the final equation used for fitting the data:

Eq. 4 100

ACCEPTED MANUSCRIPT

where,

CR IP T



AN US

Eq. 4 can be used directly to fit E vs t data to obtain estimates of the model parameters, using a global fit with two independent variables, time and agonist concentration (Section 3.4., Fig. 68A, Supplementary information file, “Fig 8 Model 4 global fit”). Here an alternative

M

method is developed, utilizing familiar empirical equations. First, the E vs t data for each agonist

PT

ED

concentration is fit individually to the general rise-and-fall exponential equation,

CE

where C is a fitted parameter and K1 and K2 are rate constants. It is evident by comparison with

AC

Eq. 4 that C is equal to,

101

ACCEPTED MANUSCRIPT

This is the same equation as that for Model 1, except “Slope” in Model 1 is replaced with C. The data analysis now proceeds as described for Model 1 (Appendix B.1.): The C versus log[A] data are fit to the sigmoid curve equation. kτ is the C value at infinite [A]. KA is the antilog of the x value at the midpoint of the curve, log(A50). The application of this method to experimental data

CR IP T

is presented in Section 3.4., Fig. 8 and the Supplementary information file, “Fig 8 Model 4 individual rise-and-fall exponential fits.”

In rise-and-fall responses, investigators typically use the response value at the peak of the E vs t curve in quantifying drug effect. Unfortunately, this peak response parameter has little

AN US

value in the kinetic model for quantifying efficacy or affinity. In the kinetic model, the peak response is defined by a complicated function of the response parameters. The time of peak

M

response is given by,

)

ED

n(

(

(

)

AC

CE

PT

and the response at peak time given by,

(

)

)

It isn‟t obvious how this expression can be used to quantify kτ and KA when [A] is the independent variable.

B.5. Consistency with the operational model of agonism

102

ACCEPTED MANUSCRIPT

The validity of a pharmacological model can be assessed by comparison with a validated model. The most widely-used pharmacological model of GPCR response is the operational model of agonism (Black and Leff, 1983). The model provides a quantitative description of efficacy. The foundation of the operational model is that the response is a hyperbolic function of

CR IP T

the concentration of agonist-occupied receptors. Here the kinetic model of GPCR response is compared with the operational model. It is shown that the operational model emerges from the kinetic model when the transduction potential is depleted (Models 3 and 4). In addition, the

AN US

operational model equation can be formulated based on the concept of transduction potential.

B.5.1. Emergence of the operational model from the kinetic model

In the operational model, response is typically measured at a single time point (Black and

M

Leff, 1983; Kenakin and Beek, 1980; Van der Graaf and Stam, 1999; Zernig et al., 1996) For the purposes of comparison, the kinetic model can be evaluated at a single time point. This is

ED

achieved by considering the time-independent term in the equations. This term is the multiplier of the bracketed term in the E vs t equations (see below).

PT

The models that assume depletion of transduction potential are Models 3 and 4. (In

CE

Model 3, E is recycled to EP whereas in Model 4 it is degraded.) It is evident from the equations for these models that the dependence of response on [RA] is hyperbolic. The E vs t equation for

AC

Model 3 is:

(

The time-independent term is, 103

)

ACCEPTED MANUSCRIPT

It is immediately evident that this expression is hyperbolic in terms of [RA]. (Formally, it is a

)

AN US

(

CR IP T

linear rational function of [RA]). Similarly, for Model 4 the E vs t equation is:

M

from which it is evident that the time-independent expression is hyperbolic in terms of [RA],

ED

It can also be shown that, at steady-state, Model 3 reduces to the operational model. The steady-state response in Model 3 is the horizontal asymptote, determined as the limit of the E vs t

AC

CE

PT

equation as t approaches infinity. For Model 3, this steady-state limit is,

The term EP(TOT) is the maximum amount of signal that can be generated within the system. This is the same parameter as Em in the operational model; Em is described as “The maximum effect in a system.” (Black and Leff, 1983) Consequently, we can substitute EP(TOT) for Em, giving,

104

ACCEPTED MANUSCRIPT

CR IP T

On rearranging we obtain,

The term kD / kE is a steady-state parameter because it does not contain time as a unit – the unit is

AN US

receptor concentration. kD / kE is the concentration of RA necessary to elicit a response that is 50% of Em. In the operational model, the same parameter is present in the equation; it is termed

ED

M

KE. Consequently, we can substitute kD / kE for KE, giving,

PT

This equation, from the kinetic model, is the same as the equation of the operational model

AC

CE

(Black and Leff, 1983), given below:

Consequently, Model 3 at steady-state reduces to the operational model. If we compare the derivations further it emerges that the transducer constant τ in the operational model is equivalent to the following function in kinetic Model 3:

105

ACCEPTED MANUSCRIPT

CR IP T

Alternatively, kτ can be expressed as a function of τ:

Unfortunately, the approach of removing the exponential terms by taking the limit as t

AN US

approaches infinity cannot be applied to Model 4 because this limit is zero.

B.5.2. Deriving the operational model from the concept of transduction potential The operational model was founded on an empirical observation and a deduction (Black

M

and Leff, 1983). The observation was that response is a hyperbolic function of agonist

ED

concentration. The deduction was that if the E vs A relationship is hyperbolic, the relationship between receptor occupancy and response must also be hyperbolic. Here it is shown that the

PT

operational model equation can be formulated from the concept of the transduction potential. Specifically, depletion of the transduction potential. This derivation is analogous to that of a

CE

particular form of the operational model, the “Hyperbolic followed by linear: simple receptor-

AC

transducer model” in (Black and Leff, 1983). In addition, it has been shown in a steady-state model that depletion of a transducer protein by RA gives rise to the operational model (Hall, 2006).

Applying the law of mass action, the response E is equal to the product of the transduction potential, agonist-occupied receptor concentration, and an equilibrium constant 1/K. This gives the equilibrium equation for E: 106

ACCEPTED MANUSCRIPT

CR IP T

The goal is an equation of the form used in the operational model, i.e. one that expresses E as a function of [RA] and constants. In the equilibrium equation above, E is a function of [RA] and a second variable, EP. EP is variable because it is depleted by formation of E. The formulation proceeds by replacing EP with an expression containing E and a constant. This is achieved by

AN US

applying the conservation of mass equation. The total amount of response material, EP(TOT), is

M

constant, and is defined by the conservation of mass equation,

ED

Solving this equation for EP and substituting into the equilibrium equation for E gives,

)

CE

PT

(

Solving this equation for E gives an equation of the desired form, one in which E is expressed as

AC

a function of [RA] and constants:

107

ACCEPTED MANUSCRIPT

This equation is analogous to the operational model equation. Em in the operational model is the maximum system effect. EP(TOT) is conceptually analogous – it is the total amount of transduction potential available for generation of response. KE in the operational model is the concentration of RA that elicits 50% of EM. Likewise, K in the equation above is the concentration of RA that

CR IP T

elicits a response that is 50% of EP(TOT). This finding shows the equation of the operational model

AC

CE

PT

ED

M

AN US

can be formulated from the concept of the transduction potential.

108

ACCEPTED MANUSCRIPT

Appendix C. Receptor reserve in the kinetic models Receptor reserve relates the fraction of receptor occupancy by agonist to the fraction of pharmacological response. In a system of receptor reserve, the fraction of maximal response to the agonist is greater than the fraction of maximal receptor occupancy by the agonist. Maximal

CR IP T

response can result from only partial receptor occupancy. The concentration of agonist necessary to elicit 50% of the response (the EC50) is less than the concentration required to occupy 50% of the receptors (the affinity, KA). By contrast, in a system without receptor reserve, the fraction of maximal response to agonist is equal to the fraction of receptor occupancy by agonist. Maximal

AN US

effect requires maximal occupancy, and the EC50 is equal to KA (i.e. potency equals affinity). The importance of receptor reserve (alternatively referred to as spare receptors) in pharmacological systems has been evaluated extensively (Kenakin, 1997). Here, two aspects are considered. First,

M

receptor reserve is a fundamental property of pharmacological response systems. It is a manifestation of an excess of response-generating capacity in the system. Second, receptor

ED

reserve complicates quantification of drug efficacy and affinity. With respect to efficacy, in a system of receptor reserve two agonists can have similar maximal effects but quite different

CE

PT

relative efficacies. These considerations are elaborated below.

C.1. Presence or absence of receptor reserve in the kinetic models

AC

Whether or not the kinetic models contain receptor reserve can be determined intuitively by defining this property in terms of the system as opposed to the receptor. Receptor reserve results from an excess of receptors relative to the response-generating capacity of the system. From the system perspective, receptor reserve results from insufficient response-generating capacity in the system for the receptor to generate its full response. This has been demonstrated

109

ACCEPTED MANUSCRIPT

for equilibrium pharmacological models (Hall, 2006; Strickland and Loeb, 1981). Consequently, receptor reserve can be viewed as system insufficiency. In the kinetic models, the responsegenerating capacity of the system is the transduction potential. The extent of system insufficiency is the extent of depletion of transduction potential. If there is abundant transduction

CR IP T

potential, it will not be depleted by generation of response. Reciprocally, if there is insufficient transduction potential, it will be by depleted by its transformation to the response. Viewed from this perspective, it can be seen that receptor reserve, i.e. system insufficiency, does not occur in Models 1 and 2 because there is no depletion of transduction potential in these models (Sections

AN US

3.1. and 3.2.). By contrast, receptor reserve is inherent in Models 3 and 4 because transduction potential is depleted in these models (Sections 3.3. and 3.4.).

The absence of receptor reserve is evident from the equations for Models 1 and 2

ED

M

(Appendix B.1. and B.2.):

CE

PT

Eq. 1

AC

Eq. 2

In both cases, response is directly proportional to the concentration of receptors, i.e. E = [RA]  constant, a relationship indicative of the absence of receptor reserve. In addition, expressing response as a fraction of the maximal response to agonist ( 110

) gives the following equation,

ACCEPTED MANUSCRIPT

CR IP T

Visual inspection of this equation indicates the fraction of response (the left-hand side) is equal

AN US

to the fraction of occupancy (the right-hand side). Finally, the equation defining EC50 is,

indicating potency is equal to affinity. This equivalence indicates there is no receptor reserve in the system for Models 1 and 2.

The presence of receptor reserve can be deduced from inspecting the equations for

M

Models 3 and 4. In the derivation, the E vs t equation is derived for [RA] (Appendix B.3 and

ED

B.4.). These equations are, for Models 3 and 4 respectively,

(

)

)

AC

CE

PT

(

These equations show that response is hyperbolically dependent on the agonist-occupied receptor concentration, a relationship which results in receptor reserve (Black and Leff, 1983;

111

ACCEPTED MANUSCRIPT

Trzeciakowski, 1999). For Model 3, receptor reserve is evident from the equation defining EC50

CR IP T

at steady-state (i.e. at infinite t). This equation is,

This equation indicates the higher the efficacy of the agonist (the higher the value of kτ), the lower the value of EC50, as predicted by models of receptor reserve (Black and Leff, 1983;

AN US

Trzeciakowski, 1999). (Unfortunately, for Model 4 an explicit equation defining EC50 cannot be derived).

With respect to detection of receptor reserve, in principle this can be done by plotting the

M

EC50 of the response vs the time of response measurement. However, this approach requires that the assumption of rapid agonist-receptor equilibration be correct. Making this assumption, in the

ED

absence of receptor reserve (Models 1 and 2) there is no change of agonist EC 50 over time (Fig. 1D and 3D). By contrast, in models that incorporate receptor reserve (Models 3 and 4) the EC50

PT

changes over time (Fig. 5D and 7D). For Model 3, the decrease of EC50 over time approaches a

CE

limit (Fig. 5D). By 20 min, the EC50 is reduced to 9.1 nM. From this value and the KA (100 nM) it is calculated that 50% effect is elicited by only 8% occupancy. As stated above, these

AC

considerations require agonist to rapidly equilibrate with the receptor. When the kinetics of agonist-receptor interaction are slow enough to impact the E vs t profile, the EC50 changes over time for all four Models (Fig. 18). Consequently, observing the change of EC50 over time cannot be used to infer the absence or presence of receptor reserve, unless it can be shown that agonist rapidly equilibrates with receptor (see Section 6).

112

ACCEPTED MANUSCRIPT

C.2. Impact of receptor reserve on quantifying agonist efficacy In single time point functional assays, measuring relative efficacy for a range of agonists is complicated in a system of receptor reserve. Typically, the only independent variable is the

CR IP T

agonist concentration, and the data are plotted as an agonist concentration-response curve. The empirical parameter from this curve that is used to quantify efficacy is the response at maximally-stimulating concentrations of agonist, termed maximal agonist effect (Neubig et al., 2003). Unfortunately, a series of agonists with a range of relative efficacy can elicit the same

AN US

maximal agonist effect. In steady-state models such as the operational model, this problem is solved by introducing a second independent variable, by introducing additional pharmacological response systems. For example, the receptor concentration is varied using irreversible

M

antagonists (Furchgott, 1966; Van der Graaf and Stam, 1999; Zernig et al., 1996) or molecular biology techniques (Wilson et al., 1996). Alternatively, the coupling efficiency is varied by

ED

measuring functional responses at different steps in a signaling pathway (Stott et al., 2016), for example in the cascade of G-protein activation, cAMP production, and cellular impedance

PT

(Peters et al., 2007). This approach yields estimates of efficacy parameters that account for the

CE

degree of receptor reserve by incorporating the receptor concentration, for example the transducer ratio τ in the operational model (τ = [R]TOT / KE).

AC

In the kinetic models, a second independent variable is already incorporated into the response measurement, this being time. Fortunately, adding this variable allows measurement of an efficacy parameter that accounts for receptor reserve. This parameter is the transduction rate constant, kτ, and its measurement for the four kinetic models is described in Section 3. kτ

113

ACCEPTED MANUSCRIPT

accounts for receptor reserve because it incorporates the receptor concentration, being equal to

AC

CE

PT

ED

M

AN US

CR IP T

EP(TOT)[R]TOTkE.

114

ACCEPTED MANUSCRIPT

Appendix D. Equations for kinetics of antagonist competition vs agonist response

Here equations are derived for measuring the kinetics of antagonist-receptor interaction in a functional assay. In the experiment, antagonist and agonist are applied according to various

CR IP T

schedules. The response to agonist is measured at various times. The response time course data are then fitted to E vs t equations to obtain estimates of the antagonist association and dissociation rate constants. Here the E vs t equations are derived first for Model 2: No EP depletion, E decays (the horizontal exponential E vs t profile, with t1/2 independent of [A], Fig.

AN US

3B, Section 3.2.). The scheme is below, where B is antagonist, and k3 and k4 are the antagonist association and dissociation rate constants, respectively. The other parameters are defined in

AC

CE

PT

ED

M

Table 1.

The E vs t equation for Model 2 is derived from the differential equations for E and [RB], using the Laplace transform method (Hoare, 2017; Mayersohn and Gibaldi, 1970). The derivation is divided into three steps: differential equation; Laplace transform; and taking the inverse Laplace

115

ACCEPTED MANUSCRIPT

transform to obtain the E vs t equation. Next, an equation is derived that assumes rapid antagonist-receptor equilibration. Finally, the E vs t equation is derived for Model 1, no EP depletion, no E decay. Note that equations incorporating antagonist binding kinetics were not derived for Models 3 and 4 because in these models the response generation mechanism is

CR IP T

second order with respect to time (see Appendix A.3.). The following assumptions apply: The kinetic model is appropriate for the response mechanism under study. Agonist is assumed to be at equilibrium with the receptor. Antagonist competitively inhibits agonist binding and binds the receptor according to a single-step interaction.

AN US

The equation used to fit the experimental data depends on the schedule of agonist and antagonist addition. Four different schedules are formularized here:

M

5. Schedule A: Simultaneous addition of agonist and antagonist. 6. Schedule B: Antagonist pre-incubation

ED

7. Schedule C: Antagonist pre-incubation and washout

CE

PT

8. Schedule D: Agonist pre-incubation

D.1. Schedule A: Simultaneous addition of agonist and antagonist

AC

In this schedule, antagonist and agonist are added simultaneously then the time course of response measured (Fig. 12A).

D.1.1. Differential equations The differential equation for E is: 116

ACCEPTED MANUSCRIPT

Since we are dealing with Model 2, the value of EP is independent of time. The time-dependent

CR IP T

variable is [RA]. It is time-dependent because [RB] is time-dependent and [RA] is linked to [RB] through [R]. Here, [RA] is substituted with an expression for [RB] and time-independent parameters, as follows:

Since A is at equilibrium with the receptor, the ratio of free R and RA is constant and can

where the constant

M

AN US

be represented by the following expression:

is fractional occupancy by A of accessible receptors, i.e. those not bound

ED

by B, and KA is the equilibrium dissociation constant for A. Next, this equation is rearranged as

CE

PT

follows:

AC

The bracketed term on the right-hand side, [R] + [RA] can be re-written in terms of [RB] and the constant [R]TOT using the conservation of mass equation for the receptor:

Rearranging,

117

ACCEPTED MANUSCRIPT

CR IP T

[R] + [RA] is now substituted with [R]TOT – [RB],

This is the desired expression: [RA] is expressed in terms of [RB] and time-independent

AN US

parameters (ρA and [R]TOT). This term is now substituted into the differential equation for E:

PT

ED

M

Next, we introduce the differential equation for antagonist binding to the receptor:

CE

This equation contains two time-dependent parameters, [R] and [RA]. We want to reduce it so that [RB] is the only time-dependent parameter. First, [R] is substituted with [RA], [RB] and

AC

[R]TOT using the conservation of mass equation:

Substituting into the differential equation,

118

ACCEPTED MANUSCRIPT

Simplifying yields, (

AN US

CR IP T

Next, [RA] is substituted with ρA([R]TOT – [RB]):

)

PT

ED

association rate of antagonist:

M

where kB(app) is the association rate of antagonist in the presence of agonist, i.e. the apparent

CE

This is the desired differential equation for [RB] because [RB] is the only time-dependent term.

AC

With the desired differential equations in hand, we move onto the next step, the Laplace transforms.

D.1.2. Laplace transforms The Laplace transforms for E is obtained as described in (Hoare, 2017; Mayersohn and Gibaldi, 1970): 119

ACCEPTED MANUSCRIPT

̅̅̅̅̅̅

̅

̅

The transform contains two time-dependent terms, E and [RB]. The goal is a transform in which

CR IP T

E is the only time-dependent term. This is achieved by taking the Laplace transform for [RB], solving it for [RB] and substituting it into the transform for E. The Laplace transform for [RB] is,

)

̅̅̅̅̅̅

AN US

(

̅̅̅̅̅̅

Solving for ̅̅̅̅̅̅ gives the expression that can be substituted into the transform for E:

(

)

(

)

ED

M

̅̅̅̅̅̅

PT

We now substitute the transform for [RB] into that for E and solve for E:

(

̅

)

AC

CE

(

)

This the desired equation – the Laplace transform for E in which E is the only time-dependent variable.

D.1.3. Inverse Laplace transform The inverse Laplace transform is now taken to obtain the E vs t equation: 120

ACCEPTED MANUSCRIPT

)

(

CR IP T

(

)

AN US

This equation is now simplified, to the final form used for curve fitting:

*

+

M

where E[B]=0, t∞ is the response to agonist at steady-state, i.e. the horizontal asymptote, in the

PT

ED

absence of antagonist:

AC

CE

The macroscopic parameters are,

(

)

121

ACCEPTED MANUSCRIPT

D.1.4. Rapid antagonist equilibration In Fig. S1D (Supplementary information), the effect of antagonist on agonist response kinetics is simulated for an antagonist that equilibrates rapidly with the receptor. Under these conditions, k4 is much greater than kD and the term

is effectively zero. Applying these limits to Eq. 5

AN US

CR IP T

gives the equation used to simulate the data in Fig. S1D:

where KB is the antagonist equilibrium dissociation constant.

D.1.5. Equation for Model 1, no EP depletion, no E decay

M

The differential equation for E for the linear E vs t profile (Model 1) is the same as that for the

ED

horizontal exponential profile (Model 2) except the EkD term is deleted because E is does not

CE

PT

degrade.

AC

The Laplace transform of this differential equation for E is,

̅̅̅̅̅̅

̅

122

ACCEPTED MANUSCRIPT

Following the same steps used for Model 2 above (Appendix D.1.3.), we arrive at the final Laplace transform in which E is the only time-dependent variable:

(

)

(

)

Taking the inverse Laplace transform yields the E vs t equation:

)

(

)

AN US

(

CR IP T

̅

M

This equation is now simplified, to the final form used for curve fitting:

(

)(

)+

PT

ED

*

AC

CE

where G[B]=0 is the gradient of the E vs t plot in the absence of antagonist, defined as,

D.2. Schedule B: Antagonist pre-incubation In this experimental schedule, antagonist is incubated with the receptor system for a certain period of time. Agonist is then added and then the time course of response measured (Fig. 12B). 123

ACCEPTED MANUSCRIPT

D.2.1. Differential equations The differential equations for E and [RB] are the same as those for Schedule A (Appendix

)

AN US

(

CR IP T

D.1.1.):

D.2.2. Laplace transforms

The Laplace transform for E is the same as that for Schedule A (Appendix D.1.2.): ̅̅̅̅̅̅

̅

ED

M

̅

PT

The Laplace transform of the differential equation for [RB] is,

(

)

̅̅̅̅̅̅

CE

̅̅̅̅̅̅

AC

where [RB]t=0 is the concentration of [RB] achieved at the end of the pre-incubation with antagonist, i.e. at the time of agonist addition (t = 0). Solving for ̅̅̅̅̅̅ gives,

(

̅̅̅̅̅̅

(

124

) )

ACCEPTED MANUSCRIPT

Substituting into the transform for E and solving for E gives,

CR IP T

where ρB(t=0) is fractional occupancy of the receptor by antagonist at the time of agonist addition:

(

̅

)

(

)

AN US

(

)

This the desired equation – the Laplace transform for E in which E is the only time-dependent

ED

D.2.3. Inverse Laplace transform

M

variable.

The Laplace transform for E is now converted to the E vs t equation by taking the inverse

AC

CE

PT

Laplace transform:

( (

)

(

)

This equation is simplified, to the final form used for curve fitting: 125

)

ACCEPTED MANUSCRIPT

*

+

CR IP T

Eq. 6

E[B]=0, t∞, ρA, kB(app), and N are defined in Schedule A (Appendix D.1.3.), and P defined as,

)

AN US

(

where,

)

M

(

PT

ED

where tPI is the time interval of the antagonist pre-incubation phase.

D.2.4. Rapid antagonist equilibration

CE

In Fig. S2D (Supplementary information), the effect of antagonist on agonist response kinetics is

AC

simulated for an antagonist that equilibrates rapidly with the receptor. Under these conditions, k4 is much greater than kD (so kD / k4 is effectively zero) and the term

is effectively zero.

Applying these limits to Eq. 6 gives the equation used to simulate the data in Fig. S2D:

126

ACCEPTED MANUSCRIPT

D.2.5. Equation for Model 1, no EP depletion, no E decay

CR IP T

The Laplace transform for E is the same as that for Schedule A (Appendix D.1.5.):

̅̅̅̅̅̅

̅

AN US

Following the same steps used for Model 2 above (Appendix D.2.3.), we arrive at the final Laplace transform in which E is the only time-dependent variable:

(

̅

)

(

)

M

(

)

PT

ED

Taking the inverse Laplace transform yields the E vs t equation:

CE

( )

(

)

AC

(

)

This equation is now simplified, to the final form used for curve fitting:

*

(

)(

127

)+

ACCEPTED MANUSCRIPT

where G[B]=0 is the gradient of the E vs t plot in the absence of antagonist, defined in Appendix

D.3. Schedule C: Antagonist pre-incubation and washout

CR IP T

D.1.5.

In this schedule, receptor is first exposed to antagonist. After a suitable incubation period to allow antagonist binding to the receptor, unbound antagonist is removed by washing under conditions in which the antagonist-receptor complex is preserved. Agonist is then added and the

AN US

response measured. (An alternative or complementary approach is to use a high concentration of agonist, sufficient to out-compete unbound antagonist from binding to the receptor.) This

D.3.1. Differential equations

M

approach is illustrated schematically in Fig. 12C.

CE

PT

ED

The differential equation for E is the same as that in Schedule A (Appendix D.1.1.):

AC

The differential equation for [RB] describes dissociation of B from RB:

D.3.2. Laplace transforms

128

ACCEPTED MANUSCRIPT

The Laplace transform for E is the same as that for Schedule A (Appendix D.1.2.):

The Laplace transform for [RB] is,

̅̅̅̅̅̅

AN US

̅̅̅̅̅̅

̅

CR IP T

̅̅̅̅̅̅

̅

where [RB]t=0 is [RB] at the moment of agonist addition. Solving for [RB] gives,

M

̅̅̅̅̅̅

ED

where ρB(t=0) is fractional occupancy of the receptor by antagonist at the time of agonist addition.

̅

AC

CE

PT

Substituting into the transform for E and solving for E gives,

This the desired equation – the Laplace transform for E in which E is the only time-dependent variable.

D.3.3. Inverse Laplace transform

129

ACCEPTED MANUSCRIPT

The Laplace transform for E is now converted to the E vs t equation by taking the inverse

CR IP T

Laplace transform:

This equation is now simplified, to the final form used for curve fitting:

AN US

Eq. 7

where E[B]=0, t∞, is defined in Schedule A (Appendix D.1.4) and Q is,

M

where ρB(t=0) is fractional receptor occupancy by antagonist at the moment of agonist addition,

ED

defined in Schedule B (Appendix D.2.2.).

In the experiment a control is run in which the antagonist is absent to determine E[B]=0,t∞

AC

CE

PT

and kD. The equation used to fit these data, derived from Eq. 7 when ρB(t=0) = 0, is,

D.3.4. Rapid antagonist equilibration When antagonist equilibrates rapidly with receptor, k4 is much greater than kD (so kD / k4 is effectively zero) and the term

is effectively zero. Applying these limits to Eq. 7 gives

the following equation: 130

ACCEPTED MANUSCRIPT

Note that this is the equation for the response to the agonist in the absence of antagonist (Eq. 2,

of the response by the agonist in Schedule C.

AN US

D.3.5. Equation for Model 1, no EP depletion, no E decay

CR IP T

Section 3.2.). This indicates that an antagonist that dissociates rapidly has no effect on generation

The Laplace transform for E is the same as that for Schedule A (Appendix D.1.5.):

̅̅̅̅̅̅

ED

M

̅

̅̅̅̅̅̅

CE

PT

The transform for [RB] is the same as that for Model 2 (Appendix D.3.2.):

AC

Following the same steps used for Model 2, we arrive at the final Laplace transform in which E is the only time-dependent variable:

̅

131

ACCEPTED MANUSCRIPT

CR IP T

Taking the inverse Laplace transform yields the E vs t equation:

This equation is now simplified, to the final form used for curve fitting:

]

AN US

[

where G[B]=0 is the gradient of the E vs t plot in the absence of antagonist, defined in Appendix D.1.5, and ρB(t=0) is fractional occupancy of receptor by antagonist at the moment of agonist

ED

M

addition, defined in Appendix D.2.2.

D.4. Schedule D: Agonist pre-incubation

PT

In this schedule, agonist is pre-incubated with receptor during which time a response will be

CE

produced. Antagonist is then added and the time-course of response measured, during which time the response to agonist will decrease owing to inhibition of agonist binding to the receptor

AC

by antagonist. This approach is illustrated schematically in Fig. 12D. Here t = 0 is the time of antagonist addition, in contrast to Schedules A, B and C in which t = 0 is the time of agonist addition.

D.4.1. Differential equations

132

ACCEPTED MANUSCRIPT

The differential equations for E and [RB] are the same as those for Schedule A (Appendix

(

CR IP T

D.1.1.):

)

D.4.2. Laplace transforms

AN US

In taking the Laplace transform for E, we include the amount of E formed during the preincubation, i.e. the level of E at the time of antagonist addition:

̅̅̅̅̅̅

̅

ED

M

̅

̅̅̅̅̅̅

( (

) )

CE

PT

The Laplace transform for [RB] is the same as that in Schedule A (Appendix D.1.2.):

AC

Substituting this transform for [RB] into the transform for E gives the transform in which E is the only time-dependent variable:

(

̅

(

133

) )

ACCEPTED MANUSCRIPT

D.4.3. Inverse Laplace transform

)

(

AN US

(

CR IP T

The inverse Laplace transform is now taken to obtain the E vs t equation:

)

M

This equation is re-arranged for curve fitting as follows. First, Et=0 is expressed in terms of the

PT

ED

model parameters and the pre-incubation time, tPI:

AC

CE

Now this expression is substituted into the E vs t equation:

(

)

) 134

(

ACCEPTED MANUSCRIPT

Now the equation is re-arranged to:

+

CR IP T

*

AN US

It is recommended that the preincubation step be long enough for the horizontal asymptote to be closely approached, in order to simplify the data analysis. Under these conditions, the term e-kDtPI

ED

[

M

is zero and the equation above reduces to the final form used for curve fitting:

]

Eq. 8

PT

E[B]=0, t∞, ρA, kB(app), M and N are defined in Schedule A (Appendix D.1.3). The equation for the

CE

control (absence of antagonist) for the post-antagonist addition can be derived by setting [B]

AC

equal to zero (so kB(app) = k4), yielding,

D.4.4. Rapid antagonist equilibration

135

ACCEPTED MANUSCRIPT

In Fig. S4D (Supplementary information), the effect of antagonist on agonist response kinetics is simulated for an antagonist that equilibrates rapidly with the receptor. Under these conditions, k4 is much greater than kD (so kD / k4 is effectively zero) and the term

is effectively zero.

]

AN US

[

CR IP T

Applying these limits to Eq.8 gives the equation used to simulate the data in Fig. S4D:

D.4.5. Equation for Model 1, no EP depletion, no E decay

In taking the Laplace transform for E, we include the amount of E formed during the pre-

̅̅̅̅̅̅

ED

̅

M

incubation, i.e. the level of E at the time of antagonist addition:

PT

Following the same steps used for Model 2 (Appendix D.4.3.), we arrive at the final Laplace

AC

CE

transform in which E is the only time-dependent variable:

(

̅

(

Taking the inverse Laplace transform yields the E vs t equation:

136

) )

ACCEPTED MANUSCRIPT

(

)

(

)

This equation is re-arranged for curve fitting as follows. First, E(t=0) is expressed in terms of the

CR IP T

model parameters and the pre-incubation time, tPI:

AN US

Now this expression is substituted into the E vs t equation:

)

(

)

ED

M

(

PT

Now the equation is re-arranged to the final form used for curve fitting:

*

(

CE

,

AC

In the absence of antagonist, the response is given by,

137

)(

)+-

ACCEPTED MANUSCRIPT

where G[B]=0 is the gradient of the E vs t plot in the absence of antagonist, defined in Appendix

AC

CE

PT

ED

M

AN US

CR IP T

D.1.5.

138

ACCEPTED MANUSCRIPT

Appendix E. Equations for kinetic response models with slow agonist-receptor equilibration

In Section 7, the kinetic response models are extended to incorporate slow agonist-receptor

CR IP T

equilibration (Models 5-8). In Models 1-4, [RA] is independent of time because RA equilibrates rapidly relative to the rate of response generation. In Models 5-8, [RA] is a time-dependent variable because RA formation is sufficiently slow that it delays response generation. This can be handled by extending the derivation used for Models 1-4. The general differential equation for E

AN US

is the same (Appendix B). An additional differential equation is required, that for [RA]. This expression is incorporated into the equation for E using the Laplace transform method (Hoare, 2017; Mayersohn and Gibaldi, 1970). The Laplace transform for E and for [RA] is taken, then the

M

transform for [RA] substituted into that for E. The inverse Laplace transform is then taken to yield the E vs t equation. This was done for Models 5 and 6, those without depletion of

ED

transduction potential. E vs t equations were not obtained for Models 7 and 8, those with depletion of transduction potential, because the response generation mechanism is second order

PT

with respect to time (Appendix A.3.). Specifically, both variables in the response generation

CE

mechanism, [RA] and EP, are time-dependent in Models 7 and 8. The model is described by the following scheme, where k1 is the association rate constant

AC

for agonist, and k2 the dissociation rate constant:

139

AN US

CR IP T

ACCEPTED MANUSCRIPT

E.1. Model 5: No EP depletion, no E decay, slow agonist equilibration This model is the analogue of Model 1; there is no depletion of response source and the response

ED

M

does not decay (Section 3.1.). The differential equation for E is that given in Appendix B.1:

PT

This equation contains two time-dependent terms, E and [RA]. The goal is an equation with one time-dependent term, the one being measured in the experiment, E. This is done as follows. First,

AC

CE

the differential equation for [RA] is written:

This equation contains two time-dependent terms, [R] and [RA]. [R] can be substituted with an expression for [RA], derived from the conservation of mass equation for the receptor:

140

ACCEPTED MANUSCRIPT

CR IP T

Solving for [RA]:

AN US

Substituting into the differential equation for [RA] and rearranging gives,

where,

̅̅̅̅̅̅

̅

PT

ED

M

Next, the Laplace transform of the differential equation for E and for [RA] is taken.

̅̅̅̅̅̅

CE

̅̅̅̅̅̅

AC

It is now straightforward to obtain an expression for E in which E is the only time-dependent variable. First, the transform for [RA] is solved for ̅̅̅̅̅̅:

̅̅̅̅̅̅

(

141

)

ACCEPTED MANUSCRIPT

We now substitute the transform for [RA] into that for E and solve for E:

(

)

CR IP T

̅

This is the desired form of the expression for E, one in which E is the only time-dependent variable. Now the E vs t equation is obtained by taking the inverse Laplace transform, as

AN US

described in (Hoare, 2017):

(

)

ED

M

Substituting EP[R]TOTkE with kτ and rearranging gives the final equation used for fitting the data:

(

)

CE

PT

Eq. 9

E.2. Model 6: No EP depletion, E decay, slow agonist equilibration

AC

This model is the analogue of Model 2; there is no depletion of response source and the response decays (Section 3.2.). The differential equation for E is given in Appendix B.2:

142

ACCEPTED MANUSCRIPT

We now follow the same steps used for Model 5 (see above, Appendix E.1.). The final Laplace transform for E is,

̅

)

CR IP T

(

The E vs t equation is obtained by taking the inverse Laplace transform:

)

Eq. 10

M

AN US

(

ED

E.3. Model 7: EP depletion, E recycles to EP, slow agonist equilibration This model is the analogue of Model 3; there is depletion of response source and the response

PT

recycles (Section 3.3.). In this model and Model 8, an E vs t equation was not obtained (except at saturating agonist concentrations – see below) because the response generation process is second

CE

order with respect to time (see Appendix A.3.). Instead, E vs t data were simulated by numerical solution of the differential equations, as described in [(Vauquelin et al., 2001), Appendix A.2.].

AC

The Euler method equations used are,

143

ACCEPTED MANUSCRIPT

When a saturating concentration of agonist is used, it can be assumed occupancy is achieved sufficiently quickly that [RA] is independent of time. The differential equation for this scenario

CR IP T

is,

AN US

The E vs t equation is obtained by following the steps in Appendix B.3:

(

)

)

Eq. 11

ED

M

(

PT

E.4. Model 8: EP depletion, E decay, slow agonist equilibration This model is the analogue of Model 4; there is no depletion of response source, and the response

CE

decays. As described for Model 7 above, the analytic E vs t equation was not obtained (except at

AC

saturating agonist concentrations – see below). Instead E vs t data were simulated by numerical solution of the differential equations [(Vauquelin et al., 2001), Appendix A.2.). The Euler method equations are:

144

ACCEPTED MANUSCRIPT

As described for Model 7, it is possible to obtain an E vs t equation for a saturating concentration

AN US

CR IP T

of agonist. The differential equation for E is,

The E vs t equation is obtained by following the steps in Appendix B.4.:

)

Eq. 12

PT

ED

M

(

E.5. Test for rapid vs slow agonist-receptor binding kinetics

CE

In Fig. 17, an experimental paradigm is shown for testing for rapid vs slow agonist-receptor

AC

binding kinetics in the functional response. In this experiment, receptor is exposed to agonist for a certain period during which the response is generated (the stimulation phase). A large excess of competitive antagonist is then added, which blocks further agonist association with the receptor. The response continues to be recorded after antagonist addition (the antagonist block phase). Here equations are derived for defining the response during the antagonist block phase.

145

ACCEPTED MANUSCRIPT

For rapid agonist-receptor binding kinetics, it is assumed the dissociation rate of the agonist-receptor complex is sufficiently large that [RA] immediately drops to zero upon addition of the antagonist. Consequently, there is no further generation of the response. Mathematically, this is represented by setting the value of EP[RA]kE to zero in the differential equations for

CR IP T

Models 1 to 4. The E vs t equations are obtained by following the derivations for Models 1 to 4 in Appendix B. Here, t is time after antagonist addition, and Et=0 is the level of response at the time of antagonist addition.

AN US

Model 1:

ED

M

Model 2, 3 and 4:

For slow agonist-receptor binding kinetics (Models 5-8), [RA] declines according to the

PT

dissociation rate of the agonist-receptor complex. For Models 5 and 6, these equations are

CE

derived in the section on measuring agonist residence time (Appendix E.6.):

AC

Model 5:

where,

(

146

)

ACCEPTED MANUSCRIPT

Model 6: (

)

CR IP T

For slow agonist-receptor binding kinetics in Models 7 and 8, data in Fig. 17 were generated by numerical solution of the differential equation using the Euler method [(Vauquelin et al., 2001), see Appendix A.2.]. This was done by setting EP[RA]kE to zero in the differential equations in

AN US

Appendix E.3. and E.4., giving:

ED

M

Model 7:

AC

CE

PT

Model 8:

E.6. Equations for measuring agonist residence time

147

ACCEPTED MANUSCRIPT

In this experimental paradigm, agonist is exposed to receptor for a certain time resulting in accumulation of response, then antagonist is added and the time course of the response measured (see diagram in Fig 18). A saturating concentration of antagonist is used, to preclude any further stimulation of signaling by the agonist-occupied receptor. The decay of the response subsequent

CR IP T

to antagonist addition is a consequence of dissociation of the agonist-receptor complex. By formularizing this mechanism, it is possible to measure the dissociation rate constant of the

PT

ED

M

AN US

agonist, k2. The scheme is,

The equations are derived first for Model 6, no EP depletion, E decay. Next, the equations are

CE

derived for Model 5, no EP depletion, no E decay. Analytical equations were not derived for Models 7 and 8 because in these models the mechanism is second order with respect to time

AC

(Appendix A.3.).

E.6.1. Differential equations The differential equations used in the derivation are,

148

CR IP T

ACCEPTED MANUSCRIPT

where t is the time after antagonist addition. Note that in the equation for [RA] it is assumed that [B] is in very large excess, such that any free receptor resulting from dissociation of RA or RB is

present in the differential equation for [RA].

M

E.6.2. Laplace transforms

AN US

instantaneously bound by B. Consequently, the term for association of A with R ([R][A]k1) is not

ED

The Laplace transforms are,

̅̅̅̅̅̅̅̅̅̅̅̅̅

PT

̅

̅̅̅̅̅̅̅̅̅̅̅̅̅

CE

̅̅̅̅̅̅̅̅̅̅̅̅̅

̅

AC

Note that the transforms incorporate the amount of the species at the time of antagonist addition, Et=0 and [RA]t=0. Solving the [RA] transform gives,

̅̅̅̅̅̅̅̅̅̅̅̅̅

149

ACCEPTED MANUSCRIPT

This expression is now substituted into the transform for E:

̅

̅

CR IP T

Solving for E gives,

AN US

̅

E.6.3. Inverse Laplace transform

ED

M

The E vs t equation is now obtained, by taking the inverse Laplace transform for E:

The data analysis is simplified if the response at the moment of antagonist addition is at the

PT

horizontal asymptote, i.e. at steady state. Under these conditions, agonist binding has reached

AC

CE

equilibrium and the term EP[RA]t=0kE is defined as,

Since EP[R]TOTkE is kτ this reduces to,

150

ACCEPTED MANUSCRIPT

This equation can be written in terms of the response to agonist. If antagonist is added at the horizontal plateau, response to agonist at the moment of antagonist addition is defined by the

CR IP T

limit of Eq. 2 as time approaches infinity (Appendix B.2.):

AN US

This equation is rearranged as follows:

M

Note the right hand-side of this equation is the same as the right-hand side of the equation for

PT

ED

EP[RA]t=0kE above. Consequently, if antagonist is added at the horizontal plateau, it follows that,

AC

CE

The right-hand side is substituted for EP[RA]t=0kE in the E vs t equation above, giving,

Rearranging gives the final E vs t equation for curve fitting:

151

ACCEPTED MANUSCRIPT

(

)

E.6.4. Equation for Model 5, no EP depletion, no E decay

CR IP T

Equation 13

This derivation follows the same steps as that for Model 5 above. The only difference is that the

AN US

expression for E decay is not present in the differential equation or Laplace transform for E, since in Model 5 E does not decay (Section 7.1.). The differential equation and Laplace

̅̅̅̅̅̅̅̅̅̅̅̅̅

PT

ED

̅

M

transform are, respectively,

̅

AC

CE

The final Laplace transform for E is,

Taking the inverse transform gives the E vs t equation:

152

ACCEPTED MANUSCRIPT

The term EP[RA]t=0kE is handled as follows. At the moment of antagonist addition, i.e. at t = 0, [RA] is defined by the ligand association equation:

)

CR IP T

(

where k1(obs) = [A]k1 + k2 and tPI is the time period of the stimulation phase. Substituting into the

AN US

E vs t equation gives the final equation for curve-fitting:

(

ED

M

where,

)

AC

CE

PT

Equation 14

153

ACCEPTED MANUSCRIPT

References

AC

CE

PT

ED

M

AN US

CR IP T

Aranyi, P., 1980. Kinetics of the hormone-receptor interaction. Competition experiments with slowly equilibrating ligands. Biochim Biophys Acta 628, 220-7. Aranyi, P., 1982. Dependence of rate constants of the glucocorticoid hormone-receptor interaction on steroid structure. J Steroid Biochem 17, 137-41. Arunlakshana, O., Schild, H. O., 1959. Some quantitative uses of drug antagonists. Br J Pharmacol Chemother 14, 48-58. Baker, J. G., 2010. The selectivity of beta-adrenoceptor agonists at human beta1-, beta2- and beta3-adrenoceptors. Br J Pharmacol 160, 1048-61, doi:10.1111/j.14765381.2010.00754.x. Baker, J. G., Hall, I. P., Hill, S. J., 2004. Temporal characteristics of cAMP response elementmediated gene transcription: requirement for sustained cAMP production. Mol Pharmacol 65, 986-98, doi:10.1124/mol.65.4.986. Benet, L. Z., Turi, J. S., 1971. Use of general partial fraction theorem for obtaining inverse laplace transforms in pharmacokinetic analysis. J Pharm Sci 60, 1593-4. Black, J. W., Leff, P., 1983. Operational models of pharmacological agonism. Proc R Soc Lond B Biol Sci 220, 141-62. Bot, I., Ortiz Zacarias, N. V., de Witte, W. E., de Vries, H., van Santbrink, P. J., van der Velden, D., Kroner, M. J., van der Berg, D. J., Stamos, D., de Lange, E. C., Kuiper, J., AP, I. J., Heitman, L. H., 2017. A novel CCR2 antagonist inhibits atherogenesis in apoE deficient mice by achieving high receptor occupancy. Sci Rep 7, 52, doi:10.1038/s41598-01700104-z. Colquhoun, D., 1998. Binding, gating, affinity and efficacy: the interpretation of structureactivity relationships for agonists and of the effects of mutating receptors. Br J Pharmacol 125, 924-47. Copeland, R. A., 2016. The drug-target residence time model: a 10-year retrospective. Nat Rev Drug Discov 15, 87-95, doi:10.1038/nrd.2015.18. Copeland, R. A., Pompliano, D. L., Meek, T. D., 2006. Drug-target residence time and its implications for lead optimization. Nat Rev Drug Discov 5, 730-9, doi:10.1038/nrd2082. Cusack, K. P., Wang, Y., Hoemann, M. Z., Marjanovic, J., Heym, R. G., Vasudevan, A., 2015. Design strategies to address kinetics of drug binding and residence time. Bioorg Med Chem Lett 25, 2019-27, doi:10.1016/j.bmcl.2015.02.027. De Lean, A., Stadel, J. M., Lefkowitz, R. J., 1980. A ternary complex model explains the agonist-specific binding properties of the adenylate cyclase-coupled beta-adrenergic receptor. J Biol Chem 255, 7108-17. Dowling, M. R., Charlton, S. J., 2006. Quantifying the association and dissociation rates of unlabelled antagonists at the muscarinic M3 receptor. Br J Pharmacol 148, 927-37, doi:10.1038/sj.bjp.0706819. Ferrandon, S., Feinstein, T. N., Castro, M., Wang, B., Bouley, R., Potts, J. T., Gardella, T. J., Vilardaga, J. P., 2009. Sustained cyclic AMP production by parathyroid hormone receptor endocytosis. Nat Chem Biol 5, 734-42, doi:10.1038/nchembio.206. Fleck, B. A., Hoare, S. R., Pick, R. R., Bradbury, M. J., Grigoriadis, D. E., 2012. Binding kinetics redefine the antagonist pharmacology of the corticotropin-releasing factor type 1 receptor. J Pharmacol Exp Ther 341, 518-31, doi:10.1124/jpet.111.188714.

154

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

AN US

CR IP T

Frace, A. M., Mery, P. F., Fischmeister, R., Hartzell, H. C., 1993. Rate-limiting steps in the betaadrenergic stimulation of cardiac calcium current. J Gen Physiol 101, 337-53. Fredriksson, R., Lagerstrom, M. C., Lundin, L. G., Schioth, H. B., 2003. The G-protein-coupled receptors in the human genome form five main families. Phylogenetic analysis, paralogon groups, and fingerprints. Mol Pharmacol 63, 1256-72, doi:10.1124/mol.63.6.1256. Frolik, C. A., Black, E. C., Cain, R. L., Satterwhite, J. H., Brown-Augsburger, P. L., Sato, M., Hock, J. M., 2003. Anabolic and catabolic bone effects of human parathyroid hormone (1-34) are predicted by duration of hormone exposure. Bone 33, 372-9. Furchgott, R. F., 1966. The use of -haloalkylamines in the differentiation of receptors and in the determination of dissociation constants of receptor-agonist complexes. Advances in Drug Research 3, 21-55. Gilman, A. G., 1987. G proteins: transducers of receptor-generated signals. Annu Rev Biochem 56, 615-49, doi:10.1146/annurev.bi.56.070187.003151. Guo, D., Mulder-Krieger, T., AP, I. J., Heitman, L. H., 2012. Functional efficacy of adenosine A(2)A receptor agonists is positively correlated to their receptor residence time. Br J Pharmacol 166, 1846-59, doi:10.1111/j.1476-5381.2012.01897.x. Guo, D., Hillger, J. M., AP, I. J., Heitman, L. H., 2014. Drug-target residence time-a case for G protein-coupled receptors. Med Res Rev 34, 856-92, doi:10.1002/med.21307. Hall, D. A., 2006. Predicting dose-response curve behavior: Mathematical models of allosteric receptor-ligand interactions. In: Bowery, N. G., (Ed.), Allosteric receptor modulation in drug targeting. CRC Press, Boca Raton, Florida, pp. 39-77. Hill, A. V., 1910. The possible effects of the aggregation of the molecules of haemoglobin on its dissociation curves. J. Physiol. 40, 4-7. Hoare, S. R. J., 2017. Receptor binding kinetics equations: Derivation using the Laplace transform method. J Pharmacol Toxicol Methods, doi:10.1016/j.vascn.2017.08.004. Hodgkin, A. L., Huxley, A. F., 1952. A quantitative description of membrane current and its application to conduction and excitation in nerve. J Physiol 117, 500-44. Hoffmann, C., Castro, M., Rinken, A., Leurs, R., Hill, S. J., Vischer, H. F., 2015. Ligand Residence Time at G-protein-Coupled Receptors - Why We Should Take Our Time To Study It. Mol Pharmacol 88, 552-60, doi:10.1124/mol.115.099671. Hothersall, J. D., Brown, A. J., Dale, I., Rawlins, P., 2016. Can residence time offer a useful strategy to target agonist drugs for sustained GPCR responses? Drug Discov Today 21, 90-6, doi:10.1016/j.drudis.2015.07.015. Hothersall, J. D., Guo, D., Sarda, S., Sheppard, R. J., Chen, H., Keur, W., Waring, M. J., AP, I. J., Hill, S. J., Dale, I. L., Rawlins, P. B., 2017. Structure-Activity Relationships of the Sustained Effects of Adenosine A2A Receptor Agonists Driven by Slow Dissociation Kinetics. Mol Pharmacol 91, 25-38, doi:10.1124/mol.116.105551. Kenakin, T., 2002. Drug efficacy at G protein-coupled receptors. Annu Rev Pharmacol Toxicol 42, 349-79, doi:10.1146/annurev.pharmtox.42.091401.113012. Kenakin, T., 2015. Gaddum Memorial Lecture 2014: receptors as an evolving concept: from switches to biased microprocessors. Br J Pharmacol 172, 4238-53, doi:10.1111/bph.13217. Kenakin, T., Jenkinson, S., Watson, C., 2006. Determining the potency and molecular mechanism of action of insurmountable antagonists. J Pharmacol Exp Ther 319, 710-23, doi:10.1124/jpet.106.107375.

155

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

AN US

CR IP T

Kenakin, T. P., 1997. Stiumulus-response mechanisms. Pharmacologic analysis of drug-receptor interaction. Lippincott-Raven Publishers, Philadelphia, pp. 78-105. Kenakin, T. P., Beek, D., 1980. Is prenalterol (H133/80) really a selective beta 1 adrenoceptor agonist? Tissue selectivity resulting from differences in stimulus-response relationships. J Pharmacol Exp Ther 213, 406-13. Klein Herenbrink, C., Sykes, D. A., Donthamsetti, P., Canals, M., Coudrat, T., Shonberg, J., Scammells, P. J., Capuano, B., Sexton, P. M., Charlton, S. J., Javitch, J. A., Christopoulos, A., Lane, J. R., 2016. The role of kinetic context in apparent biased agonism at GPCRs. Nat Commun 7, 10842, doi:10.1038/ncomms10842. Lacourcière, Y., Asmar, R., 1999. A comparison of the efficacy and duration of action of candesartan cilexetil and losartan as assessed by clinic and ambulatory blood pressure after a missed dose, in truly hypertensive patients: a placebo-controlled, forced titration study. Candesartan/Losartan study investigators. Am J Hypertens 12, 1181-7. Langmuir, I., 1918. The adsorption of gases on plane surfaces of glass, mica and platinum. J. Am. Chem. Soc. 40, 1361-1403. Liu, A. M., Ho, M. K., Wong, C. S., Chan, J. H., Pau, A. H., Wong, Y. H., 2003. Galpha(16/z) chimeras efficiently link a wide range of G protein-coupled receptors to calcium mobilization. J Biomol Screen 8, 39-49, doi:10.1177/1087057102239665. Lohse, M. J., Nikolaev, V. O., Hein, P., Hoffmann, C., Vilardaga, J. P., Bunemann, M., 2008. Optical techniques to analyze real-time activation and signaling of G-protein-coupled receptors. Trends Pharmacol Sci 29, 159-65, doi:10.1016/j.tips.2007.12.002. Louvel, J., Guo, D., Agliardi, M., Mocking, T. A., Kars, R., Pham, T. P., Xia, L., de Vries, H., Brussee, J., Heitman, L. H., Ijzerman, A. P., 2014. Agonists for the adenosine A1 receptor with tunable residence time. A Case for nonribose 4-amino-6-aryl-5-cyano-2thiopyrimidines. J Med Chem 57, 3213-22, doi:10.1021/jm401643m. Malany, S., Hernandez, L. M., Smith, W. F., Crowe, P. D., Hoare, S. R., 2009. Analytical method for simultaneously measuring ex vivo drug receptor occupancy and dissociation rate: application to (R)-dimethindene occupancy of central histamine H1 receptors. J Recept Signal Transduct Res 29, 84-93, doi:10.1080/10799890902721339. Manglik, A., Kobilka, B., 2014. The role of protein dynamics in GPCR function: insights from the beta2AR and rhodopsin. Curr Opin Cell Biol 27, 136-43, doi:10.1016/j.ceb.2014.01.008. Mayersohn, M., Gibaldi, M., 1970. Mathematical methods in pharmacokinetics. I. Use of the Laplace transform for solving differential rate equations. Vol. 34, American Journal of Pharmaceutical Education, pp. 608-614. Mayersohn, M., Gibaldi, M., 1971. Mathematical methods in pharmacokinetics. II. Solution of the Two Compartment Open Model. Vol. 35, American Journal of Pharmaceutical Education, pp. 19-28. Michaelis, L., Menten, M. L., Johnson, K. A., Goody, R. S., 2011. The original Michaelis constant: translation of the 1913 Michaelis-Menten paper. Biochemistry 50, 8264-9, doi:10.1021/bi201284u. Motulsky, H. J., 2017. Analysis checklist: Comparing nonlinear fits. Vol. 2017. Motulsky, H. J., 2018a. How Prism generates random numbers. Accessed 14 January 2018. Motulsky, H. J., 2018b. Simulating a XY data table. Accessed 14 January 2018. Motulsky, H. J., Mahan, L. C., 1984. The kinetics of competitive radioligand binding predicted by the law of mass action. Mol Pharmacol 25, 1-9. 156

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

AN US

CR IP T

Motulsky, H. J., Mahan, L. C., Insel, P. A., 1985. Radioligand, agonists and membrane receptors on intact cells: data analysis in a bind. Trends in Pharmacological Sciences 6, 317-319. Mould, R., Brown, J., Marshall, F. H., Langmead, C. J., 2014. Binding kinetics differentiates functional antagonism of orexin-2 receptor ligands. Br J Pharmacol 171, 351-63, doi:10.1111/bph.12245. Neubig, R. R., Spedding, M., Kenakin, T., Christopoulos, A., International Union of Pharmacology Committee on Receptor, N., Drug, C., 2003. International Union of Pharmacology Committee on Receptor Nomenclature and Drug Classification. XXXVIII. Update on terms and symbols in quantitative pharmacology. Pharmacol Rev 55, 597-606, doi:10.1124/pr.55.4.4. Paton, W. D. M., Rang, H. P., 1966. A kinetic approach to the mechanism of drug action. In: Harper, N. J., Simmonds, A. B., Eds.), Advances in Drug Research. Academic Press, New York, pp. 57-80. Peters, M. F., Knappenberger, K. S., Wilkins, D., Sygowski, L. A., Lazor, L. A., Liu, J., Scott, C. W., 2007. Evaluation of cellular dielectric spectroscopy, a whole-cell, label-free technology for drug discovery on Gi-coupled GPCRs. J Biomol Screen 12, 312-9, doi:10.1177/1087057106298637. Pohl, S. L., Birnbaumer, L., Rodbell, M., 1969. Glucagon-sensitive adenyl cylase in plasma membrane of hepatic parenchymal cells. Science 164, 566-7. Princen, K., Hatse, S., Vermeire, K., De Clercq, E., Schols, D., 2003. Evaluation of SDF1/CXCR4-induced Ca2+ signaling by fluorometric imaging plate reader (FLIPR) and flow cytometry. Cytometry A 51, 35-45, doi:10.1002/cyto.a.10008. Ramsey, S. J., Attkins, N. J., Fish, R., van der Graaf, P. H., 2011. Quantitative pharmacological analysis of antagonist binding kinetics at CRF1 receptors in vitro and in vivo. Br J Pharmacol 164, 992-1007, doi:10.1111/j.1476-5381.2011.01390.x. Rask-Andersen, M., Masuram, S., Schioth, H. B., 2014. The druggable genome: Evaluation of drug targets in clinical trials suggests major shifts in molecular class and indication. Annu Rev Pharmacol Toxicol 54, 9-26, doi:10.1146/annurev-pharmtox-011613-135943. Riccobene, T. A., Omann, G. M., Linderman, J. J., 1999. Modeling activation and desensitization of G-protein coupled receptors provides insight into ligand efficacy. J Theor Biol 200, 207-22, doi:10.1006/jtbi.1999.0988. Riddy, D. M., Valant, C., Rueda, P., Charman, W. N., Sexton, P. M., Summers, R. J., Christopoulos, A., Langmead, C. J., 2015. Label-Free Kinetics: Exploiting Functional Hemi-Equilibrium to Derive Rate Constants for Muscarinic Receptor Antagonists. Mol Pharmacol 88, 779-90, doi:10.1124/mol.115.100545. Rodbell, M., Lin, M. C., Salomon, Y., 1974. Evidence for interdependent action of glucagon and nucleotides on the hepatic adenylate cyclase system. J Biol Chem 249, 59-65. Rohatgi, A., 2018. WebPlotDigitizer. Web based tool to extract data from plots, images, and maps. Accessed 4 February 2018. Rosethorne, E. M., Bradley, M. E., Gherbi, K., Sykes, D. A., Sattikar, A., Wright, J. D., Renard, E., Trifilieff, A., Fairhurst, R. A., Charlton, S. J., 2016. Long Receptor Residence Time of C26 Contributes to Super Agonist Activity at the Human beta2 Adrenoceptor. Mol Pharmacol 89, 467-75, doi:10.1124/mol.115.101253. Samama, P., Cotecchia, S., Costa, T., Lefkowitz, R. J., 1993. A mutation-induced activated state of the beta 2-adrenergic receptor. Extending the ternary complex model. J Biol Chem 268, 4625-36. 157

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

AN US

CR IP T

Shea, L., Linderman, J. J., 1997. Mechanistic model of G-protein signal transduction. Determinants of efficacy and effect of precoupled receptors. Biochem Pharmacol 53, 519-30. Slack, R. J., Hall, D. A., 2012. Development of operational models of receptor activation including constitutive receptor activity and their use to determine the efficacy of the chemokine CCL17 at the CC chemokine receptor CCR4. Br J Pharmacol 166, 1774-92, doi:10.1111/j.1476-5381.2012.01901.x. Sriram, K., Insel, P. A., 2018. GPCRs as targets for approved drugs: How many targets and how many drugs? Mol Pharmacol, doi:10.1124/mol.117.111062. Stott, L. A., Hall, D. A., Holliday, N. D., 2016. Unravelling intrinsic efficacy and ligand bias at G protein coupled receptors: A practical guide to assessing functional data. Biochem Pharmacol 101, 1-12, doi:10.1016/j.bcp.2015.10.011. Strange, P. G., 2010. Use of the GTPgammaS ([35S]GTPgammaS and Eu-GTPgammaS) binding assay for analysis of ligand potency and efficacy at G protein-coupled receptors. Br J Pharmacol 161, 1238-49, doi:10.1111/j.1476-5381.2010.00963.x. Strickland, S., Loeb, J. N., 1981. Obligatory separation of hormone binding and biological response curves in systems dependent upon secondary mediators of hormone action. Proc Natl Acad Sci U S A 78, 1366-70. Sungkaworn, T., Jobin, M. L., Burnecki, K., Weron, A., Lohse, M. J., Calebiro, D., 2017. Singlemolecule imaging reveals receptor-G protein interactions at cell surface hot spots. Nature 550, 543-547, doi:10.1038/nature24264. Sykes, D. A., Dowling, M. R., Charlton, S. J., 2009. Exploring the mechanism of agonist efficacy: a relationship between efficacy and agonist dissociation rate at the muscarinic M3 receptor. Mol Pharmacol 76, 543-51, doi:10.1124/mol.108.054452. Sykes, D. A., Parry, C., Reilly, J., Wright, P., Fairhurst, R. A., Charlton, S. J., 2014a. Observed drug-receptor association rates are governed by membrane affinity: the importance of establishing "micro-pharmacokinetic/pharmacodynamic relationships" at the beta2adrenoceptor. Mol Pharmacol 85, 608-17, doi:10.1124/mol.113.090209. Sykes, D. A., Dowling, M. R., Leighton-Davies, J., Kent, T. C., Fawcett, L., Renard, E., Trifilieff, A., Charlton, S. J., 2012. The Influence of receptor kinetics on the onset and duration of action and the therapeutic index of NVA237 and tiotropium. J Pharmacol Exp Ther 343, 520-8, doi:10.1124/jpet.112.194456. Sykes, D. A., Riddy, D. M., Stamp, C., Bradley, M. E., McGuiness, N., Sattikar, A., Guerini, D., Rodrigues, I., Glaenzel, A., Dowling, M. R., Mullershausen, F., Charlton, S. J., 2014b. Investigating the molecular mechanisms through which FTY720-P causes persistent S1P1 receptor internalization. Br J Pharmacol 171, 4797-807, doi:10.1111/bph.12620. Traynor, J. R., Clark, M. J., Remmers, A. E., 2002. Relationship between rate and extent of G protein activation: comparison between full and partial opioid agonists. J Pharmacol Exp Ther 300, 157-61. Trzeciakowski, J. P., 1999. Stimulus amplification, efficacy, and the operational model. Part I-binary complex occupancy mechanisms. J Theor Biol 198, 329-46, doi:10.1006/jtbi.1999.0919. Van der Graaf, P. H., Stam, W. B., 1999. Analysis of receptor inactivation experiments with the operational model of agonism yields correlated estimates of agonist affinity and efficacy. J Pharmacol Toxicol Methods 41, 117-25.

158

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

AN US

CR IP T

Vauquelin, G., Van Liefde, I., 2006. Slow antagonist dissociation and long-lasting in vivo receptor protection. Trends Pharmacol Sci 27, 356-9, doi:10.1016/j.tips.2006.05.001. Vauquelin, G., Morsing, P., Fierens, F. L., De Backer, J. P., Vanderheyden, P. M., 2001. A twostate receptor model for the interaction between angiotensin II type 1 receptors and nonpeptide antagonists. Biochem Pharmacol 61, 277-84. Vilardaga, J. P., Bunemann, M., Feinstein, T. N., Lambert, N., Nikolaev, V. O., Engelhardt, S., Lohse, M. J., Hoffmann, C., 2009. GPCR and G proteins: drug efficacy and activation in live cells. Mol Endocrinol 23, 590-9, doi:10.1210/me.2008-0204. Waelbroeck, M., 2001. Activation of guanosine 5'-[gamma-(35)S]thio-triphosphate binding through M(1) muscarinic receptors in transfected Chinese hamster ovary cell membranes; 1. Mathematical analysis of catalytic G protein activation. Mol Pharmacol 59, 875-85. Waelbroeck, M., Tastenoy, M., Camus, J., Christophe, J., 1991. Binding kinetics of quinuclidinyl benzilate and methyl-quinuclidinyl benzilate enantiomers at neuronal (M1), cardiac (M2), and pancreatic (M3) muscarinic receptors. Mol Pharmacol 40, 413-20. Weiss, J. M., Morgan, P. H., Lutz, M. W., Kenakin, T. P., 1996. The cubic ternary complex receptor-occupancy model. I. Model description. J Theor Biol 178, 158-167, doi:10.1006/jtbi.1996.0139. Whalen, E. J., Rajagopal, S., Lefkowitz, R. J., 2011. Therapeutic potential of beta-arrestin- and G protein-biased agonists. Trends Mol Med 17, 126-39, doi:10.1016/j.molmed.2010.11.004. Wilson, S., Chambers, J. K., Park, J. E., Ladurner, A., Cronk, D. W., Chapman, C. G., Kallender, H., Browne, M. J., Murphy, G. J., Young, P. W., 1996. Agonist potency at the cloned human beta-3 adrenoceptor depends on receptor expression level and nature of assay. J Pharmacol Exp Ther 279, 214-21. Wittmann, H. J., Seifert, R., Strasser, A., 2011. Influence of the N-terminus and the E2-loop onto the binding kinetics of the antagonist mepyramine and the partial agonist phenoprodifen to H(1)R. Biochem Pharmacol 82, 1910-8, doi:10.1016/j.bcp.2011.09.005. Woodroffe, P. J., Bridge, L. J., King, J. R., Hill, S. J., 2009. Modelling the activation of Gprotein coupled receptors by a single drug. Math Biosci 219, 32-55, doi:10.1016/j.mbs.2009.02.003. Zafar, M. A., Droege, C., Foertsch, M., Panos, R. J., 2014. Update on ultra-long-acting beta agonists in chronic obstructive pulmonary disease. Expert Opin Investig Drugs 23, 1687701, doi:10.1517/13543784.2014.942730. Zernig, G., Issaevitch, T., Woods, J. H., 1996. Calculation of agonist efficacy, apparent affinity, and receptor population changes after administration of insurmountable antagonists: comparison of different analytical approaches. J Pharmacol Toxicol Methods 35, 223-37.

159

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

AN US

CR IP T

Graphical Abstract

160