Analyzing kinetic signaling data for G-protein-coupled receptors

In classical pharmacology, bioassay data are fit to general equations (e.g. the dose response equation) to determine empirical drug parameters (e.g. EC50 and Emax), which are then used to calculate chemical parameters such as affinity and efficacy. Here we used a similar approach for kinetic, time course signaling data, to allow empirical and chemical definition of signaling by G-protein-coupled receptors in kinetic terms. Experimental data are analyzed using general time course equations (model-free approach) and mechanistic model equations (mechanistic approach) in the commonly-used curve-fitting program, GraphPad Prism. A literature survey indicated signaling time course data usually conform to one of four curve shapes: the straight line, association exponential curve, rise-and-fall to zero curve, and rise-and-fall to steady-state curve. In the model-free approach, the initial rate of signaling is quantified and this is done by curve-fitting to the whole time course, avoiding the need to select the linear part of the curve. It is shown that the four shapes are consistent with a mechanistic model of signaling, based on enzyme kinetics, with the shape defined by the regulation of signaling mechanisms (e.g. receptor desensitization, signal degradation). Signaling efficacy is the initial rate of signaling by agonist-occupied receptor (kτ), simply the rate of signal generation before it becomes affected by regulation mechanisms, measurable using the model-free analysis. Regulation of signaling parameters such as the receptor desensitization rate constant can be estimated if the mechanism is known. This study extends the empirical and mechanistic approach used in classical pharmacology to kinetic signaling data, facilitating optimization of new therapeutics in kinetic terms.

Pharmacological data analysis provides the drug activity parameters used to optimize the biological activity of new therapeutics and to identify mechanisms of receptor-mediated signal transduction for G-protein-coupled receptors (GPCRs). In the typical analysis process, activity is measured at various drug concentrations and the data fit to a concentration-response equation by curve fitting, for example to a sigmoid curve equation [1][2][3] . This yields empirical drug parameters, usually the potency (EC 50 ) and a measure of the maximal signaling capacity (E max ). This analysis is described as "Model-free" because it makes minimal assumptions regarding the mechanism of signal transduction. These parameters can then be used to calculate mechanistic drug parameters that define the effect of the drug on the system in chemical terms (e.g. efficacy and macroscopic affinity) [3][4][5][6][7][8] . These chemical parameters aid the translation of drug effect measurement from in vitro test systems to animal models and human disease, for example by minimizing the effect of tissue-specific parameters on the quantification of drug activity. Biased agonism measurement is an example of this process currently in use by many laboratories. Signaling assay data are fit to concentration-response equations and the empirical parameters EC 50 and E max used to calculate mechanistic parameters of biased signaling, such as transducer coefficients and log efficacy ratios (reviewed in 9 ).
The manner in which GPCR signaling changes over time, i.e. the kinetics/dynamics of response, is currently of considerable interest in quantitative pharmacology. The time course of GPCR response has been measured since the earliest studies of muscle contraction 1 , through the discovery of G-protein-mediated signaling 10 , the application of Ca 2+ indicator dyes 11 , and the most recent advances in biosensor technology that enable high-throughput signaling kinetic analysis [12][13][14] . Measuring the time course of receptor signaling revealed fundamental mechanisms of GPCR function. The classic rise-and-fall time course of cAMP in S49 cells indicated desensitization of β 2 adrenoceptor signaling 15 . Sustained signaling after agonist washout by parathyroid hormone, sphingosine 1-phosphate and thyroid-stimulating hormone receptors revealed signaling by internalized receptors [16][17][18] , a new spatiotemporal paradigm of GPCR signaling of potential therapeutic utility [19][20][21][22][23] . Signaling kinetics can also affect quantification of drug activity. Of potential concern, drug effect can be dependent on the time point at which the response is measured. For example, biased signaling of the D 2 dopamine receptor was shown to be timedependent 24 . For the CB 1 receptor, the potency for internalization increased over time 25 . For AT 1 angiotensin receptor-mediated arrestin recruitment the potency increased over time, as did the E max of partial agonists 26 . Finally, differences of signaling kinetics can be a manifestation of structural differences of receptor-ligand interaction. For example, the nonpeptide GLP-1 receptor agonist TT-OAD2 protrudes outside of the receptor core, interacting with lipid, potentially explaining the unusually slow signaling kinetics of the ligand 27 .
A data analysis framework for time course data of GPCR signaling would enable the quantification of efficacy in kinetic terms. One way to do this is to measure the initial rate of signaling, analogous to the initial rate of enzyme activity 26,[28][29][30][31][32] , an approach recently applied to receptor receptor-arrestin interaction 26 . Kinetic analysis could also enable measurement of regulation of signaling parameters, such as the rate constant for receptor desensitization. At present, time course data for GPCR signaling are rarely analyzed by curve fitting to estimate kinetic pharmacological parameters (in contrast to the near-universal application of curve fitting to concentration-dependence data). The absence of curve fitting to time course data might be resulting in lost opportunities, particularly for biosensor assay modalities in which the time course data are collected by default. Kinetic insights into GPCR function might be being missed. Potential benefits of the kinetic dimension to ligand optimization could be going unrealized. The goal of this study was to develop a straightforward data analysis framework for quantifying the kinetics of GPCR signaling, enabling investigators to measure kinetic drug parameters useful for practical pharmacological application.

Results
In this study, a data analysis framework for curve fitting of time course data was developed for GPCR signaling kinetics/dynamics. First, we identified the types of time course curve shape by conducting a comprehensive literature survey. It was discovered that the curve shapes, more precisely the temporal profiles, conform to a limited number of shapes (four). This survey also showed how the time course is shaped by regulation of signaling mechanisms (e.g. receptor desensitization and response degradation). The shapes were defined by simple equations that can be used to analyze time course data using familiar curve-fitting software (e.g. Prism, Graph-Pad Software, Inc., San Diego CA). This analysis yields a kinetic drug parameter, the initial rate of signaling by agonist-occupied receptor, that defines efficacy as the rate of signal generation before it is impacted by regulation of signaling mechanisms. This parameter, termed k τ , provides a biologically meaningful and intuitive metric of the kinetics of signal transduction. Finally, mechanistic models are applied to the time course data using a kinetic mechanistic model of agonism developed previously 33 , and extended here to incorporate receptor desensitization and sustained signaling by internalized receptors. This analysis demonstrated the model-free equations emerge as general forms of the mechanistic equations, providing a mechanistic foundation for the analysis approach.
Literature survey of time course shapes. The first step in analyzing biological activity data is visual inspection of the data. In pharmacological data analysis, this led to the realization that concentration-response data are usually described by a sigmoid curve (when the x axis is the logarithm of ligand concentration) 1,2 . Here we surveyed the GPCR signaling literature for time course data. The survey was designed to be comprehensive, spanning (1) The full range of heterotrimeric G-protein classes (G s , G i , G q/11 , G 12/13 ); (2) response times of milliseconds to hours; (3) Proximity to the receptor, from direct receptor interaction (G-protein and arrestin) to downstream signals (e.g. cell contraction/relaxation and gene expression); (4) Types of response, including chemical (second messengers), mechanical (cell structure change), electrical (ion channel currents) and genetic (gene expression). This was done for experiments when the receptor was exposed continuously to the agonist, i.e. agonist washout experiments were not considered.
Within this survey, the large majority of time course profiles conformed to one of four shapes (Fig. 1). Each shape is defined by a corresponding equation (Eqs. 1-4 in Methods). The equations are amenable to straightforward curve-fitting analysis in familiar software. The shapes are, in order of increasing complexity: 1. The straight line (Fig. 1a). Response increases continuously over time at a constant rate (defined by Eq. 1). 2. The association exponential curve (Fig. 1b, Eq. 2). Initially, response increases rapidly, being almost linear.
The response then slows then finally approaches a plateau at which the response level remains constant over time. 3. The rise-and-fall to baseline curve (Fig. 1c, Eq. 3). The first phase resembles the association exponential curve -response increases rapidly at first, then slows. The response then reaches a peak level. Subsequently, the response level decreases and finally falls back to the baseline level, i.e. the level before addition of ligand. 4. The rise-and-fall to steady-state curve (Fig. 1d, Eq. 4). Response rises and falls as for the rise-and-fall to baseline curve. However, the response declines to a steady-state level, above that of the baseline level before addition of ligand.
Straight line time course profile. The straight line profile was evident in second messenger responses under a specific condition-when signaling was unregulated, i.e. when regulation of signaling mechanisms were blocked (Fig. 2). GPCR signaling is regulated over the short term by two primary mechanisms: receptor desensitization (involving receptor phosphorylation and subsequent arrestin binding) [34][35][36] and degradation of the response (for example metabolism of second messengers) [37][38][39] . The straight line profile was observed when both mechanisms were blocked. For example, AT 1 angiotensin receptor-stimulated inositol phosphates (IP) production was linear when desensitization was blocked (using arrestin knock out cells) and response degradation was blocked (using www.nature.com/scientificreports/ Li + to block IP breakdown) (Fig. 2a) 40 . The same result was obtained for the proteinase-activated receptor PAR1 receptor (see Fig. 2a in 41 ). A second example is provided by the GnRH 1 receptor, which lacks a C-terminal tail and so does not interact with arrestin 42 . A linear time course was observed, in the presence of Li + (Fig. 2b) 43 . (See also Fig. 1a,b in 44 and Fig. 3 in 45 .) A third example is provided by the glucagon receptor. A linear time course of cAMP accumulation was observed in hepatic membranes 46 . The membranes lacked cAMP phosphodiesterase activity 47 , and likely lacked receptor desensitization components owing to extensive washing of the preparation (Fig. 2c). Linear time course data were also observed for long duration responses, far downstream, at the level of DNA, specifically gene expression and DNA synthesis ( Supplementary Fig. S1) 48,49 . In these responses there was a delay before the onset of the response. This has been rationalized for gene expression responses by a need for a build-up of signal transduction intermediates to a threshold level that then initiates the process 48 .
Association exponential time course profile. The association exponential profile was the most commonlyobserved time course shape, especially for second messenger molecules such as cAMP and inositol phosphates, but also for a variety of other signals, from upstream events such as arrestin recruitment to downstream cellular functions such as muscle cell relaxation (Fig. 3, Supplementary Fig. S2). Regarding regulation of signaling, the profile was evident when one of the two mechanisms was operative and one was blocked (the mechanisms being receptor desensitization and response degradation). Examples include AT 1 receptor-stimulated IP production, with receptor desensitization in operation (arrestin wild-type cells) but without response degradation (IP metabolism blocked by Li + ) (Fig. 3a) 40 ; GnRH 1 receptor-simulated IP production with response degradation in www.nature.com/scientificreports/ operation (no Li + present) but without receptor desensitization (the GnRH receptor lacking a C-terminal tail) (Fig. 3b) 45 ; and β 2 adrenoceptor-stimulated cAMP generation with response degradation (cAMP metabolism) but without receptor desensitization (arrestin knockout cells) (Fig. 3c) 53 . The association exponential profile was observed in direct receptor-transducer interaction assays, including receptor-G-protein interaction and receptor-arrestin interaction ( Supplementary Fig. S2a,b) 26,54 . This was expected; mechanistically, the response is likely resultant from a bimolecular interaction and the association exponential equation is a general form of the familiar bimolecular association kinetics equation. G-protein activation was also described by the profile, measured by [ 35 S]GTPγS binding, or by a G-protein activation biosensor ( Supplementary Fig. S2c,d) 54,55 . cAMP production responses also conform to the profile, both stimulation (Supplementary Fig. S2e) and inhibition ( Supplementary Fig. S2f) 56,57 . A specialized signaling pathway, β-catenin stabilization via the Wnt-Frizzled system, was also described by the profile 58 (Supplementary Fig. S2g). Electrical signaling was also described by the association exponential profile, as shown by GIRK channel gating in Xenopus In all cases, regulation of signaling activity (receptor desensitization and response degradation) was minimized as described in the text boxes. (a) Inositol phosphates production via the AT 1 receptor stimulated by 100 nM AngII in fibroblasts from arrestin knock-out mice (from Fig. 3c of 40 ). (b) Inositol phosphates production via the GnRH 1 receptor stimulated by 1 μM GnRH in HEK293 cells (from Fig. 4a of 43 ). (c) cAMP production in hepatic membranes stimulated by 2 μM glucagon (data from Fig. 1 of 46 , basal response in the absence of glucagon subtracted). Note the membrane preparation used lacks cAMP phosphodiesterase activity 47 . Data were fit to the straight line equation (Eq. 1) using Prism 8 50 . The Slope value on the panels, which is equal to the initial rate and k τ , is the fitted value ± the fit SEM 51 Fig. S2h) 59 . In the final example, a downstream cellular response was found to conform to the profile, relaxation of human airway smooth muscle cells via the β 2 adrenoceptor ( Supplementary Fig. S2i) 60 .
Rise-and-fall to baseline time course profile. The rise-and-fall to baseline profile is a classic curve shape in pharmacology, leading to the definition of "Fade" (decline in the response to a continuous application of agonist 62 ) and the discovery of the underlying mechanism (regulation of signaling processes, especially receptor desensitization) 15  (c) cAMP production via the β 2 adrenoceptor stimulated by 1 μM isoproterenol in fibroblasts from arrestin knock-out mice (data from Fig. 5d of 53 ). Note a slight decline was observed at later time points (> 10 min), which have been excluded from the figure. The fitted values SSR (steady-state response) and k (the rate constant) are the fitted value ± the fit SEM 51,52 . The initial rate, equal to k τ , was calculated as the steady state response (SSR) multiplied by the rate constant k.  Fig. 4, diacylglycerol production via the AT 1 receptor (Fig. 4a) 65 and cAMP production via the β 2 adrenoceptor 53 (Fig. 4b). In both studies the mechanisms underlying the curve shape were investigated. The rise-and-fall to baseline mechanism was evident when both regulation of signaling mechanisms were operative, i.e. when there were no experimental manipulations of receptor desensitization or degradation of the response. When the regulation mechanisms were manipulated, the shape of the curve changed. When receptor desensitization was blocked or attenuated, the resulting curve shape approached the association exponential curve, as shown in Fig. 2a of the original study for the AT 1 receptor 65 and Fig. 3c for the β 2 adrenoceptor. Response degradation was in operation, demonstrated for the β 2 adrenoceptor by the effect of phosphodiesterase inhibition 53 and assumed for the AT 1 receptor because diacylglycerol is typically cleared rapidly by diacylglycerol kinases 38 .
The rise-and-fall to baseline response is also well known in calcium signaling, representing the change of cytoplasmic Ca 2+ concentration on stimulation by GPCR agonists (usually when there is no extracellular Ca 2+ in the assay). An unusual regulation mechanism is involved in Ca 2+ signaling via intracellular stores. The amount of cytoplasmic Ca 2+ that can be obtained is limited by the amount of Ca 2+ in the intracellular stores 66,67 . Depletion of Ca 2+ in the store as it is being released into the cytoplasm limits the amount and rate of cytoplasmic Ca 2+ release. This regulation mechanism can be described as depletion of response precursor 33 . The response is also regulated by clearance of cytoplasmic Ca 2+ out of the cell 68,69 . This can be described as response degradation. These two regulation processes result in the rise-and-fall to baseline profile as described previously 33 . An example is provided by GnRH-stimulated Ca 2+ mobilization in pituitary gonadotrophs, as shown in Fig. 4c 70 . Other  Fig. 1a, absence of extracellular Ca 2+ , of 70 ). Note in (b) the fitting procedure involved minimization of outlier contribution to the fit ("Robust regression") 74 , which precluded estimation of the fitted parameter error and R 2 values (see Supplementary file "Curve fit results" for details). The initial rate, equal to k τ , is the fitted parameter C. www.nature.com/scientificreports/ receptors giving this profile include the BB1 receptor 71 and AT 1 receptor 26 . It is important to note that other time course shapes are observed for cytoplasmic Ca 2+ mobilization, including the rise-and-fall to steady-state curve (considered separately below), and calcium oscillations and waves 72,73 . In addition, rise-and-fall profiles are observed in Ca 2+ responses which are not fit well by the rise-and-fall equations used in this study (e.g. in Supplementary Fig. S3).
Rise-and-fall to steady-state time course profile. A second rise-and-fall curve shape is encountered in GPCR signaling in which the response, after rising to the peak, falls to a plateau level of response which is above baseline, i.e. response declines to a steady-state after peaking. Examples of this shape are shown in Fig. 5. In this study, an equation was identified that describes these data, termed the rise-and-fall to steady-state equation (Eq. 4). This equation was identified as a general form of explicit GPCR signaling mechanism equations, as described in Appendix 2. These mechanisms are more complex than the simplest mechanisms considered previously 33 and include recently-discovered mechanisms, such as persistent signaling by internalized receptors 16,17 . These mechanisms are considered below ("More complex models"). We have loaded the equation into a custom Prism template, "Signaling kinetic model-free equations" available in the supplementary files. The most commonly-encountered example was calcium mobilization, the increase of cytoplasmic Ca 2+ upon application of the GPCR agonist. A representative example is shown in Fig. 5a, cytoplasmic Ca 2+ concentration stimulated by GnRH via the GnRH 1 receptor 70 . It is well known that the cytoplasmic Ca 2+ level usually reaches a plateau after sustained application of the agonist, particularly when Ca 2+ is included in the extracellular medium. Under this condition, Ca 2+ can re-enter the cell and this can result in a steady-state being reached between the numerous processes controlling cytoplasmic Ca 2+ concentration 68,72,73,75 . This mechanism can be represented in the context of the regulation of signaling mechanisms as a reformation of the response from the response decay product, in this case reappearance of cytoplasmic Ca 2+ from the extracellular space. This mechanism was formulated in Appendix 2.3.3 (see supplementary files). This mechanism is described in general form by the rise-andfall to steady-state equation (Eq. 3) and in explicit form by Eq. (15). Data are analyzed by the general form (see "Measuring the initial rate from curve fit parameters" below) and explicit form ("Measuring model parameters").
The rise-and-fall to steady-state profile is also commonly observed in numerous downstream signaling responses, examples of which are shown in Fig. 5: ERK1/2 phosphorylation by lysophosphatidic acid 76 , and protein kinase C and Rho activation via the AT 1 receptor 77 . In some cases, the data fit the rise-and-fall to steadystate equation well (high R 2 value) but the fitted parameters were ambiguous, a scenario encountered when the two rate constant values were almost equal. An example is provided in Supplementary Fig. S4. This observation requires further investigation, for example using structural identifiability analysis.
initial rate measurement using model-free analysis. In order for time course analysis to be useful for pharmacological investigations, pharmacological parameters need to be extracted from the data that capture the temporal dimension of activity. One kinetic parameter used routinely for another target class (enzymes) and occasionally for GPCRs is the initial rate of activity. This is the rate at the earliest part of the time course, when the response increases linearly over time. This parameter offers several benefits. (1) It is a biologically-meaningful descriptor of signaling activity. (2) In principle, the initial rate is effectively a pure efficacy parameter, being the rate of signal generation before it is affected by regulation of signaling mechanisms, which affect the shape of the later part of the time course. This is formally demonstrated using the kinetic pharmacological model below.
(3) From a practical perspective, the initial rate provides a single parameter, as opposed to multiple parameters, of drug efficacy, simplifying application to ligand optimization in medicinal chemistry projects. (4) The initial rate can be determined regardless of the shape of the time course, enabling comparison of a ligand's efficacy between responses with different temporal profiles (see below). We recently applied the initial rate approach to arrestin-receptor interaction, a special case in which there was no signal transduction 26 . Here it is applied universally to GPCR signaling responses. It is important to note the initial rate is system-dependent because it is dependent on system parameters such as the receptor and effector concentrations (see "Comparison of models with model-free analysis" below).
Measuring the initial rate from curve fit parameters. For curved time course data, the initial rate is often measured by assessing which data points lie on the linear portion at the start of the curve, then fitting a straight line equation to those points. Here a more efficient method is developed. The entire time course is fit to the equation that defines the curve, then the initial rate calculated from the fitted parameters. This requires an additional equation that defines the initial rate ( IR ) in terms of the fitted parameters. Rise-and-fall to steady-state, from Eq. (4): www.nature.com/scientificreports/ This method was used to quantify the initial rate of signaling for the responses in Figs. 2, 3, 4 and 5 and Supplementary Figs. S1, S2. Data were fit to the appropriate time course equation using Prism 8. The fitted parameter values and the fitted standard error are shown in the panels of the figures. For linear and association exponential profiles data were analyzed using built-in equations provided with the software. For the rise-and-fall profiles, user-entered custom equations were used. (A ready-to-use Prism template containing these equations, "Signaling kinetic model-free equations, " is provided in the supplementary files.) In some cases, the signal initiation time was a variable in the equation (see Eqs. [5][6][7][8][9], allowing for delay between addition of ligand and initiation of signal (e.g. in gene expression assays, Supplementary Fig. S1), or accommodating the scenario where the ligand addition time is not precisely defined. The data presented in this study were fit well by the equations; in all but one case, the correlation coefficient was > 0.95, and the standard error of the fit parameters was typically < 10% of the fitted value (see "Curve fit results" Excel file in Supporting Material). With two exceptions, data were fit well using the default fitting method of Prism 8 (least squares regression with medium convergence criteria 78 ). The exceptions were the rise-and-fall fit in Fig. 4b (cAMP production via the β 2 adrenoceptor without phosphodiesterase inhibition) and Fig. 5b (LPA-stimulated ERK phosphorylation). In these cases, the default method  Fig. 3a,b, respectively, of 77 , fit to Eq. 8). Note in (b) the fitting procedure involved minimization of outlier contribution to the fit ("Robust regression") 74 , which precluded estimation of R 2 values. Fitted parameter values are given in the file "Curve fit results" in Supporting Material. Initial rate was calculated using the equation, Initial rate = SSR × (Dk 1 − (D − 1)k 2 ) (See Appendix 1.1). Data in (a) were also analyzed using the mechanism equation for precursor depletion and response degradation to steady-state (Eq. 17 in Methods). The fitted curve overlies that of the rise-and-fall to steady-state fit. www.nature.com/scientificreports/ yielded an ambiguous fit 79 (see "Curve fit results" Excel file in Supporting Material), which was resolved using the "Robust regression" option which minimizes the contribution of outliers to the fit 74 . The initial rate value was then calculated from the fitted parameter estimates using the equations above. Note the simplicity of obtaining the initial rate in most cases. For example, for the association exponential curve, the initial rate is simply the rate constant multiplied by the steady-state response. For the rise-and-fall to baseline curve, it is simply the value of the fitted parameter C . The calculated initial rate values are shown in the figures (Figs. 2, 3, 4 and 5, Supplementary Figs. S1, S2). A maximally-stimulating concentration of agonist was employed in all the examples provided. Under this condition, the initial rate is the efficacy parameter of the agonist in kinetic terms, here termed IR max .
It is also possible to perform global fitting of the signaling kinetic data with time and ligand concentrations as independent variables. This can require knowledge of which parameters in the model are dependent on the ligand concentration, which in turn can require knowledge of the underlying mechanism. The methods used for this approach are described in ref 33 .
Comparison with single time point analysis. Historically, pharmacological responses have been measured for multiple concentrations of ligand at a single time point and ligand activity described by the potency (e.g. EC 50 ) and maximal response (E max ). Here we compared this single time point concentration-response analysis with a concentration-response analysis of the initial rate of signaling. For this purpose, we measured cAMP accumulation and arrestin recruitment via the V 2 vasopressin receptor. This receptor is activated by two cognate endogenous ligands, vasopressin and oxytocin 80 . The receptor is well known to couple stably to arrestin when bound by vasopressin (a so-called class B arrestin binding profile 34 ). The receptor is activated by two endogenous ligands, vasopressin and oxytocin, that differ in their spatiotemporal mechanisms of signaling 81 . Geneticallyencoded biosensors were used to quantify these responses in real time in live cells expressing the receptor 26,82 , as described in Methods. Figure 6 shows the time courses for a range of concentrations of the two ligands for cAMP accumulation and arrestin recruitment. In all cases, the time course data were fit well by the association exponential equation (Eq. 6, Fig. 6, see also "Curve fit results" Excel file in Supporting Material). From the fitted parameters, the initial rate was calculated. This was done by multiplying the steady-state response by the rate constant as described above (see "Curve fit results" Excel file for values). The initial rate for maximally-stimulating concentrations of ligand is shown graphically by the dashed straight lines in Fig. 6c,f. The concentration-response of the initial rate was then evaluated, as shown in Fig. 7a. The concentration response was fit well by a standard sigmoid curve equation, from which the efficacy and potency of the ligand could be determined ( Table 1). The efficacy is the maximal initial rate ( IR max ), i.e. the upper plateau of the curve. The potency is the midpoint concentration of the ligand, specifically the concentration of ligand giving an initial rate half of the maximum value. Here this is termed L 50 ( Table 1).
The concentration-response analysis indicated oxytocin was less active for recruiting arrestin than vasopressin (Fig. 7a, Table 1), consistent with a previous report 81 . This was manifest as a 6.2-fold lower potency for oxytocin compared with vasopressin (580 nM vs 93 nM respectively, − log L 50 values significantly different, Table 1). The maximal initial rate of arrestin recruitment ( IR max ) for oxytocin was value 70% of the vasopressin value but the IR max values were not significantly different (Table 1). These results indicate the weaker arrestin recruitment by oxytocin can be accounted for by a lower affinity of oxytocin compared with vasopressin for the arrestin-bound receptor. By contrast the two ligands were effectively equivalent for stimulation of cAMP production (Fig. 6c, Table 1). These findings imply biased agonism of oxytocin relative to vasopressin (G-protein bias) in terms of the initial rate.
Next, data were analyzed using a single time point. A time late in the time course was selected when the responses for high concentrations was approaching steady-state (18 min after ligand addition, 22 min after the read start, Fig. 6). This reflects a modality in which the assay is left long enough for the maximal signal to accumulate, often referred to as an "Endpoint" modality. The concentration-response curves are shown in Fig. 7b and results of the sigmoid curve analysis in Table 2. Using this metric, the parameters were slightly different from the initial rate analysis (Fig. 7a, Table 1). Specifically, the maximal response of oxytocin was similar to that of vasopressin for both arrestin and cAMP responses. In addition, rather than being uquipotent in the cAMP assay, vasopressin was more potent than oxytocin (L 50 of 0.089 and 0.36 nM, − log L 50 values significantly different, Table 2). This means that the biased agonism in terms of potency observed in the initial rate analysis for oxytocin was not observed in the endpoint modality (since the potency difference for arrestin, 4.0-fold, was similar to the potency difference for cAMP, 6.9-fold, Table 2). One possible explanation for this difference is that the initial rate for vasopressin in generating cAMP was limited by slow receptor-ligand equilibration, a scenario discussed in theoretical terms previously 33 and in the Discussion below. Another difference was that the potency for arrestin recruitment was higher using the endpoint approach compared with the initial rate approach (by 4.4-fold for vasopressin and 3.2-fold for oxytocin, compare Tables 1 and 2). This difference was also observed for the AT 1 angiotensin receptor and reflects the kinetic mechanism of arrestin recruitment, as described in ref. 26 .

Mechanistic pharmacological model of GpcR signaling and regulation kinetics.
Recently a mechanistic kinetic pharmacological model was developed to quantify GPCR signaling activity in kinetic terms 33 that is being applied to understanding signaling efficacy in kinetic terms 8 . The model comprises a minimal number of parameters to enable parameter estimation by curve fitting. The approach is macroscopic-ligand signaling activity is defined at the level of the whole system, like the operational model of agonism 7 . (Numerous microscopic systems-biology type models have been developed to simulate GPCR signaling dynamics that quantify each interaction in the signaling cascade 8 www.nature.com/scientificreports/  Table 1. V 2 initial rate dose response for cAMP generation and arrestin recruitment stimulated by vasopressin and oxytocin via the V 2 vasopressin receptor. The initial rate of the responses were measured using geneticallyencoded biosensors 26,82 as described in Methods. Time course data were analyzed to determine the initial rate as described in Fig. 6. The initial rate was then plotted against the ligand concentration (Fig. 7a) and the data fit to a sigmoid curve equation 83 to determine the maximal initial rate (IR max , the upper plateau, equivalent to k τ ) and L 50 (concentration of ligand giving an initial rate 50% of IR max   Fig. 6 for raw data. The data were fit to a sigmoid curve equation 83 to determine the maximal initial rate (IR max , the upper plateau, equivalent to k τ ) and L 50  www.nature.com/scientificreports/ they contain too many parameters to be useful in medicinal chemistry campaigns.) The model is illustrated in Fig. 8. Signaling activity is quantified using a single, readily-measurable parameter, k τ , which is the initial rate of signaling by the agonist-bound receptor in enzymatic terms 33 . The model is extended to incorporate known regulation of signaling mechanisms. GPCR signaling is regulated to limit signaling, preventing over-stimulation of the cell. We consider the properties of the model, how it relates to the model-free analysis, and how to estimate the efficacy and regulation parameters by curve fitting. The model gives rise to the four curve shapes observed experimentally. It emerged that the specific shape is simply dependent on the number of regulatory mechanisms (0, straight line; 1, association exponential; 2, rise-and-fall to baseline). It is shown the model-free analysis emerges from the mechanistic model, the equations of the former being general forms of the equations for the latter. Estimating efficacy is straightforward and does not require knowledge of the mechanism-it is shown k τ is equal to IR max . The regulation parameters can be estimated if the mechanism is known, for example the receptor desensitization rate constant.
Model description. The kinetic mechanistic model describes GPCR signaling in terms of enzyme activity. A signal results from conversion of a signal precursor (analogous to the substrate) to the signal (the product) by the agonist-bound GPCR (the enzyme) (Fig. 8), as described previously 33 . The rate of ligand efficacy for generating the response is quantified as the initial rate of signaling, termed k τ , analogous to the initial rate of enzyme activity. k τ is the rate of response generation by the agonist-occupied receptor, defined as the product of the total receptor concentration ( [R] TOT ), total signal precursor ( E P ) and the response-generation rate constant ( k E ): (By comparison, the initial rate of enzyme activity is the product of the enzyme concentration, substrate concentration and the catalytic rate constant.) The model as formulated in this study assumes receptor-ligand occupancy does not change over time, that the equilibrium level of occupancy is rapidly achieved and is not rate-liming. This condition is likely met for maximally-stimulating ligand concentrations (used to quantify efficacy), and for all concentrations of lower potency ligands, scenarios which result in rapid association of ligand with receptor (see Discussion). Note the model framework is amenable to extension to incorporate the kinetics of receptor-ligand binding, as described previously 33 . The model can be extended to incorporate regulation of signaling. The canonical mechanisms are response degradation and receptor desensitization. Response degradation is the process by which the signal, once generated, is cleared over time. Examples include breakdown of second messenger molecules, such as cAMP by phosphodiesterases 37 , and de-activation of G-protein by hydrolysis of bound GTP to GDP by the intrinsic GTPase activity of the G-protein 89 . Some signals are decreased by clearance of the signaling species from the relevant compartment, for example efflux of cytosolic Ca 2+ ions 69 . In the model, response degradation (pink region of Fig. 8) is represented simply by exponential decay, governed by k D , the response degradation rate constant. This component of the model was described previously 33 . Figure 8. Kinetic mechanistic model of agonism incorporating regulation of signaling. Both mechanisms of short-term signaling regulation are included -receptor desensitization (yellow) and response degradation (pink). The response generation process is in green. In this process, the response precursor ( E P ) is converted to the response ( E ) by the agonist-occupied receptor ( RA ). This proceeds at a rate defined by k τ , the transduction rate constant, which is the initial rate of response generation by the agonist-occupied receptor. Receptor desensitization is represented by transformation of RA into the inactive receptor R 0 A that does not generate a response (because it can't couple to E P ). The rate constant for desensitization is k DES . Degradation of response is represented by decay of E to D , governed by the degradation rate constant k D . This model is an extension the original kinetic mechanistic model 33  www.nature.com/scientificreports/ Receptor desensitization model. In the present study, the model is extended to incorporate the second canonical regulation process, receptor desensitization ( Fig. 8 yellow region). Receptor desensitization typically results from receptor phosphorylation and subsequent arrestin binding, which inhibits G-protein interaction [34][35][36] . This process is represented simply as an exponential decay of the agonist-bound receptor concentration that can couple to the signaling machinery. This is governed by the desensitization rate constant, k DES . The basic model is formulated in Appendix 2.1 and is represented by Scheme 1 below: The time course shape predicted by this model is an association exponential curve (Fig. 9b). This is consistent with experimental results; systems in which the sole regulation mechanism is receptor desensitization yield an association exponential time course. Examples include IP production by the AT 1 (Fig. 3a) 40 , PAR1 (see Fig. 2b in 41 ) and C-terminally-extended GnRH receptor (see Fig. 3c in 43 ) (when response degradation is blocked by Li + ). This curve shape makes sense intuitively. At early times response is generated rapidly. The response then slows, because response generation is attenuated by the loss of active receptor. Ultimately the response approaches a limit (the plateau). At this limit the response level does not change because no new response is being generated and the existing response is not degraded.
This model can now be extended to incorporate both regulation mechanisms, receptor desensitization and response degradation (Appendix 2.2), represented by Scheme 2 below: The time course shape for this model is a rise-and-fall to baseline curve (Fig. 9e). This makes sense intuitively. The response rises rapidly then slows as receptor desensitization starts to slow the rate of response generation. The response becomes further limited owing to response degradation. The response reaches an upper limit (the peak) when the rate of response generation equals the rate of degradation. After this, response degradation predominates over response generation. Less and less new response is generated because the active receptor concentration is declining, and ultimately no new response will be generated because the active receptor concentration will decline to zero. Ultimately the response level falls to zero once all existing response has been degraded. This profile is evident in systems where both mechanisms have been shown to be in operation. Examples include diacylglycerol production via the AT 1 receptor (Fig. 4a) 65 , cAMP accumulation via the β 2 adrenoceptor in the absence of phosphodiesterase inhibition (Fig. 4b) 53 , and IP production via the CCK1 receptor in the absence of Li + (Fig. 10a of 90 ).
More complex models. More complex regulation and signaling mechanisms discovered experimentally can be represented using the model formulation, as described in Appendix 2.3. Receptor can resensitize and then contribute to signaling. This model is formulated in combination with response degradation in Appendix 2.3.1. The receptor, once desensitized and internalized, can continue to signal from internal compartments [16][17][18] . This model is formulated in Appendix 2.3.2. A third complex model describes calcium signaling; response is degraded (cleared from the cytoplasm) but the response can reform (from calcium that re-enters the cell). This process, in combination with depletion of the response precursor, is formulated in Appendix 2.3.3. In all three models the time course is described by a rise-and-fall to steady-state equation (see Appendix 2.3).
Comparison of models with experimental data. The receptor desensitization model completes the set of kinetic mechanistic signaling models developed for GPCRs, initiated in ref. 33 . This allows comparison of a complete kinetic model of GPCR signaling and regulation with experimental data, and with the model-free analysis. The model mechanisms are shown schematically in Fig. 9 and are: no regulation; receptor desensitization; response degradation; response recycling; receptor desensitization and response degradation; precursor depletion and response degradation; and the three more complicated models described above and in Appendix 2.3. Simulators The time course curve shapes predicted by the model are those observed experimentally -the straight line, the association exponential, and the two rise-and-fall curves (Fig. 9). It emerges that the shape is dependent on the number of regulatory mechanisms. With no regulatory mechanisms, the time course is a straight line (Fig. 9a). With one mechanism, an association exponential curve results (receptor desensitization, response degradation, or response recycling, Fig. 9b-d). With two mechanisms, a rise-and-fall to baseline curve results (Fig. 9e,f). More precisely, when an input regulation mechanism (desensitization or precursor depletion) is coupled with an output mechanism (response degradation). With the more complicated mechanisms, a rise-and-fall to baseline profile results (Fig. 9g-i). www.nature.com/scientificreports/ Comparison of models with model-free analysis. The four curve shapes of the mechanistic model are the same as those of the model-free analysis. This is shown graphically by comparing Fig. 9 with Fig. 1. It is also evident mathematically by inspection of the equations; the model equations listed in Supplementary Tables S1-S4 conform to the forms of the general equations (Eqs. [1][2][3][4]. This finding indicates the model-free analysis emerges from known mechanisms of GPCR signaling and regulation. The maximal initial rate from the model-free analysis ( IR max ) is equal to k τ from the mechanistic model. This is evident by taking the initial rate of the equations, i.e. the limit as time approaches zero. By definition this limit is IR max for the model-free analysis. For the mechanistic equations, this limit is k τ (for a maximally-stimulating concentration of agonist, Appendix 1.2). This supports the hypothesis that the initial rate is purely an efficacy term, being the rate of response generation before it is impacted by regulation of signaling mechanisms, because k τ contains efficacy terms but no regulation terms. It also provides a mechanistic foundation of IR max , indicating it is analogous to the initial of an enzyme's activity (formally, the product of the receptor and precursor concentrations and a rate constant). Finally, this formulation demonstrates the maximal initial rate is system-dependent, being dependent on the receptor and precursor levels, which can differ between systems.
Measuring model parameters. The models contain parameters for signal generation and for regulation of signaling. These parameters can be estimated by curve fitting by two methods. The time course data can be fit directly to equations that explicitly describe the model. These equations are listed in Supplementary Information and are provided in a custom Prism template in the supplementary files ("Signaling kinetic mechanism equations"). Alternatively, the data can be fit to general time course equations (Eqs. 1-4) and the fitted parameter values used to calculate the model parameters. These approaches are used here to estimate efficacy and regulation parameters, which can be done using just a maximally-stimulating concentration of agonist. Note in all the literature examples in this study, a maximally-stimulating concentration was used.
The general equation fitting approach is used here for the linear, association exponential, and rise-and-fall to baseline examples. For the linear examples, the only parameter to be estimated is k τ and this is equal to the slope of the line for a maximally-stimulating concentration of agonist. The resulting k τ values are the Slope values in Fig. 2 and Supplementary Fig. 1.
For the association exponential time course, there are two parameters, the efficacy parameter k τ , and a regulation parameter specified by the model (for example, the desensitization rate constant k DES ). k τ can be estimated using the same method as that used in the model-free analysis to calculate IR max , since these parameters are equivalent (see "Measuring the initial rate from curve fit parameters" above). The rate constant k is multiplied by the steady-state level of response, for a maximally-stimulating concentration of agonist. This resulting k τ values are equal to the initial rate values given in Fig. 3 and Supplementary Fig. 2, and the IR max values in Table 1. Note in all the model variants that give rise to the association exponential profile, k τ is the product of k and the steadystate response, i.e. k τ is calculated in the same way. This means k τ can be estimated without knowing the specific mechanism underlying the curve shape. The regulation parameter of the association exponential time course is dependent on the mechanism. In all cases, the regulation parameter is given by the rate constant parameter k of the general equation (Supplementary Table S2). If the mechanism is known, the value of k can be ascribed to a regulatory process. This is shown using the studies cited here where the regulation of signaling mechanism was evaluated. In Fig. 3a, IP generation by the AT 1 receptor 40 , the mechanism was most likely receptor desensitization since deletion of arrestin changes the shape to a straight line (Fig. 2a) and response degradation was blocked by Li + . The rate constant k then represents the desensitization rate constant k DES . The fitted value of 0.040 min −1 indicated a receptor desensitization half time of 17 min. In Fig. 3b,c, the mechanism is most likely response degradation since receptor desensitization was minimized, the GnRH receptor lacking a C-terminal tail 45 (Fig. 3b), the β 2 adrenoceptor-expressing cells devoid of arrestin 40 (Fig. 3c), and in both cases inhibitors of degradation excluded. The k value then represents the response degradation rate constant k D . The fitted values, 0.12 min −1 for the GnRH 1 receptor and 1.1 min −1 for the β 2 adrenoceptor, indicated response degradation halftimes of 5.8 min and 38 s, respectively.
For the rise-and-fall to baseline mechanisms, there are three parameters -the efficacy parameter k τ , and two regulation parameters. k τ can again be estimated using the general equation without knowing the specific mechanism underlying the time course profile. k τ is equal to the parameter C of the general equation (Eq. 3) when a maximally-stimulating agonist concentration is employed, and these values are shown for the literature examples in Fig. 4. The rate constant values can be ascribed to the mechanism rate constants if the mechanism is known. In Fig. 4a, diacyglycerol production by the AT 1 receptor, the mechanism is most likely receptor desensitization and response degradation, since blocking arrestin recruitment changes the shape to an association exponential (see Fig. 2a of ref. 65 ) and diacyclglycerol is rapidly degraded by diacylglycerol kinases 38 . It is then necessary to determine which of the general rate constants ( k 1 and k 2 ) corresponds to which of the mechanism rate constants ( k DES and k D ). In this case, it is likely the faster of the rates ( k 1 ) corresponds to k D , since blocking arrestin recruitment gives an association exponential profile with rate similar to k 1 (see Fig. 2a of ref. 65 ). This approach gives a k D ( k 1 ) value of 5.5 min −1 and a k DES ( k 2 ) value of 0.80 min −1 , corresponding to degradation and desensitization half times of 7.6 s and 52 s, respectively. In Fig. 4c, cytoplasmic Ca 2+ elevation via the GnRH receptor, the mechanism is most likely precursor depletion and response degradation, based on the known mechanism of calcium responses (see "Rise-and-fall to baseline time course profile" above). The more rapid of the two rate constants is most likely precursor depletion since this is the first event in the signaling cascade. This logic gives a k DEP ( k 1 ) value of 0.46 s −1 and a k D ( k 2 ) value of 0.072 s −1 , corresponding to depletion and degradation half times of 1.5 s and 9.6 s, respectively. Figure 4b shows cAMP generation by the β 2 adrenoceptor. In this example, likely resulting from receptor desensitization and response degradation based on the experiments in ref. 53  www.nature.com/scientificreports/ values are close to one another (1.5 and 0.96 min −1 ), preventing the ascribing of the constants to the regulatory processes but demonstrating the rate of the processes is similar in this system. For the rise-and-fall to steady-state models, k τ can be estimated using the general equation and a maximallystimulating concentration of agonist. k τ is equal to IR max , calculated from the parameters as described above. For estimating the regulation parameters, we recommend using the mechanistic model equation rather than the general equation. For this curve shape, the general equation parameters correspond to complicated combinations of the mechanism parameters, which introduces constraints of the parameter values that are not incorporated into the general model equations (Appendix 2.3). This approach was used to analyze the Ca 2+ mobilization data in Fig. 5a. In this experiment, Ca 2+ was included in the extracellular medium, which results in a steady-state being reached between cytoplasmic Ca 2+ entry and export. This is the "Response degradation to steady-state with precursor depletion" model (Appendix 2.3.3). The equation (Eq. (17), see Methods) fitted the data well (R 2 value 0.995, see curve in Fig. 5a, standard error of the fitted model parameters less than 10%-see "Curve fit results" Excel file in Supporting Material). The regulation rate constant values were 0.37 s −1 for k DEP , 0.052 s −1 for k D and 0.0033 s −1 for k R , corresponding to half times for precursor depletion, response degradation and response reformation of 1.9 s, 13 s and 3.5 min, respectively. As expected, the k τ value from the mechanism equation fit was the same as the value from the model-free analysis (240 nM s −1 ).

Discussion
New technologies have enabled efficient measurement of the kinetics of GPCR signaling using continuous-read, real-time bioassay modalities [12][13][14] . Methods are now required to extract pharmacological parameters from the time course data. Such methods have been developed for specific responses that are immediately proximal to the receptor (internalization 25 and arrestin-receptor interaction 26 ). Here methods were designed for universal application to GPCR signaling and regulation responses. Curve-fitting methods are described for analyzing time course signaling data, designed for routine use by pharmacologists, enabling application to medicinal chemistry and receptor research. These models could account for a broad array of historical time course signaling data. The equations are simple and built into commonly-used curve-fitting programs, or are provided in ready-to-use templates. The efficacy for signaling is quantified kinetically by an intuitive, familiar and biologically meaningful parameter, the initial rate of signaling, which is obtained from the curve fit parameters. The underlying theoretical framework is familiar, based on the concepts of enzyme kinetics and the known mechanisms of regulation of signaling. Resources are provided in Supporting Material to facilitate implementation and application of the models, including Prism templates containing the equations, and a time course data simulator for the mechanistic models. With these new analysis tools, the model can be tested prospectively using new data, testing the effects of regulation of signaling, evaluating how the initial rate characterization compares with historical single time point pharmacological methods, and investigating the effect of receptor reserve (see below). This approach has practical benefits for optimization of new molecules and for GPCR signaling research. In project workflow, it allows representation of the time course data set by a minimal number of informative parameters (e.g. k τ , k DES ). This enables tabulation of kinetic signaling data for a series of ligands, necessary for medicinal chemists to establish kinetic structure-activity relationships. If the response at a specific time point is required, for example when comparing in vivo and in vitro activity, this can be calculated from the curve fit parameters, avoiding the necessity of measuring precisely the same time points in every experiment. The approach solves the time-dependency problem, in which pharmacological parameter values can be dependent on the time point at which they are measured, complicating quantification of efficacy and biased agonism [24][25][26] . This is because: (1) Activity is quantified as a rate constant value, which is constant over time; (2) Differences of curve shape can be accommodated because the same efficacy parameter, the initial rate, can be extracted from all four of the commonly-encountered time course curve shapes. Finally, the method separates efficacy and regulation, allowing these parameters to be quantified separately (e.g. k τ and k DES ), enabling independent structure-activity optimization of these processes.
The mechanistic model yields insights that explain the time course curve shapes and provide mechanistic meaning to the parameters from the model-free analysis. The model and experimental data show that regulation of signaling defines the time course curve shape. Regulation changes the shape from a straight line (where response is generated indefinitely), to an association exponential curve or a rise-and-fall curve, with the type of curve dependent on the number and nature of the regulation mechanisms (Fig. 9). The underlying mechanism can be interrogated by blocking receptor desensitization and/or response degradation. The model supports the hypothesis that the initial rate of activity is purely an efficacy term, unaffected by regulation of signaling: It is shown that k τ , devoid of regulation terms, is the limit of the equations as time approaches zero, the formal definition of the initial rate (Appendix 1). The model translates the empirical parameters of the model-free analysis into biologically-meaningful parameters of GPCR signaling mechanisms. IR max is equivalent to k τ , the initial rate of signaling generation by the agonist occupied receptor, analogous to the initial rate of an enzyme's activity. The rate constants k , k 1 and k 2 correspond simply to the rate constants of the regulation processes, e.g. receptor desensitization and response degradation. If the mechanism is known, the rate constant value can be assigned to a regulatory activity, enabling regulation of signaling to be quantified in simple, kinetic terms.
The relationship between the kinetic analysis and existing pharmacological analysis is now considered. Established pharmacological analysis quantifies signaling in terms of efficacy and affinity, usually at a single time point at which equilibrium between receptor and ligand is assumed [3][4][5][6][7][8] . The efficacy term in the kinetic approach ( IR max or k τ ) is a direct analogue of established efficacy parameters (E max , and τ of the operational model 7 ). Consequently, the rank order of efficacy should be the same using both approaches. Affinity is more difficult to measure using the kinetic approach because of the equilibration issue. Properly measuring affinity from a concentration-response requires the assay to be incubated long enough for equilibrium between receptor and Scientific RepoRtS | (2020) 10:12263 | https://doi.org/10.1038/s41598-020-67844-3 www.nature.com/scientificreports/ ligand to be closely approached 91 . By contrast, the initial rate in the kinetic method is quantified using the earliest time points after the addition of ligand and under these conditions it cannot be generally assumed the ligand and receptor are at equilibrium 5,33,92,93 . The impact can be evaluated by considering typical values of receptor-ligand binding rate constants, using the mass-action receptor-ligand association kinetics equation 94 . For low affinity ligands, the equilibration time is short, potentially enabling accurate estimation of affinity when using responses that proceed over several minutes. A compound with 10 μM affinity, binding with an association rate constant ( k on ) of 10 8 M −1 min −1 , reaches 97% of its equilibrium occupancy within 0.1 s at its K d concentration (i.e. 10 μM). However, for high affinity ligands the equilibration time can be long relative to the duration of the response. A compound with 1 nM affinity and a k on of 10 8 M −1 min −1 takes 18 min to reach 97% of its equilibrium receptor occupancy when applied at its K d concentration. The problem is magnified when the response is very rapid, as previously described 92 , particularly in Ca 2+ signaling, which proceeds in seconds. This issue can be circumvented when receptor and ligand can be pre-incubated prior to the initiation of the response, for example in GTPγS binding assays, but this option is rarely available in whole cell assays, the most commonly used assay modality. Importantly, this issue does not affect estimation of efficacy. IR max and k τ are estimated using maximally-effective concentrations which, due to mass action, bind more rapidly. For example, when applied at 10 μM concentration, the 1 nM K d , 10 8 M −1 min −1 k on compound reaches 97% of equilibrium occupancy within 0.2 s. Finally, the effect of signal amplification, manifest as receptor reserve, in the kinetic mechanistic model requires further investigation. While the model can incorporate receptor reserve intrinsically (as a depletion of response precursor) 33 , the rate of signaling in a multistep pathway is ultimately a function of the rates of the individual steps and the stoichiometric relationships between them. For example, the initial rate could be susceptible to receptor reserve effects if there is a rate-limiting step upstream of the response being measured. This issue requires experimental investigation.
In conclusion, this study introduces straightforward data analysis methods to quantify the kinetics of GPCR signaling. The simplicity of the analysis procedures, the intuitive nature of the parameters, and the mechanistic foundation in known and emerging GPCR signaling paradigms will facilitate pharmacological discovery and optimization in kinetic terms. where SSR is the steady-state response, specifically the ligand-specific response as time approaches infinity. Note the y asymptote value as time approaches infinity is SSR + Baseline. k is the observed rate constant in units of t −1 . The corresponding equation in Prism 8 is named "One-phase association, " 61 in which Span corresponds to SSR, Y0 corresponds to Baseline, and Plateau corresponds to SSR + Baseline. The rise-and-fall to baseline equation is Eq. (3):

Methods
where C is a fitting constant in units of y-units.t −1 (which is the initial rate of signaling, see Appendix 1.1) and k 1 and k 2 are observed rate constants in units of t −1 . In the analysis, k 1 is the larger of the two rate constant values, constrained to be greater than k 2 (see "Fitting procedures" below). This equation has been loaded into a Prism template as a user-defined equation named "Rise-and-fall to baseline time course" (see "Signaling kinetic modelfree equations" file available in the supplementary files). The rise-and-fall to steady-state equation is Eq. (4): where, as above, SSR is the steady-state response, specifically the ligand-specific response as time approaches infinity. Note the y asymptote value as time approaches infinity is SSR + Baseline. D is a unitless fitting constant, and k 1 and k 2 are observed rate constants in units of t −1 . In the analysis, k 1 is the larger of the two rate constant values, constrained to be greater than k 2 (see "Fitting procedures" below). This equation has been loaded into a Prism template as a user-defined equation named "Rise-and-fall to steady state time course" in which the www.nature.com/scientificreports/ parameter SteadyState corresponds to SSR (see "Signaling kinetic model-free equations" file available in the supplementary files). In certain circumstances it is desirable to float the initiation of signaling time in the analysis. This enables the analysis to accommodate slight uncertainty as to the precise time point of ligand addition, especially for rapidly-generated signals. It also allows for signaling mechanisms where there is a delay between ligand-receptor binding and initiation of detectable signaling, often observed in gene expression assays, and thought to be due to a necessary build-up of signal transduction intermediates to a threshold level 48 . When floating the start time in the analysis it is necessary to establish the baseline response level prior to the start time. Equations 1-4 above can be adapted to incorporate a fitted start time, the corresponding equations being Eqs. (5-8) below: where t 0 is the signal initiation time and t is time uncorrected for the signal initiation time. Equation (6) is equivalent to the equation named "Plateau followed by one phase association" in Prism 8, in which Span corresponds to SSR and Y0 corresponds to Baseline 95 . Equations 5,7 and 8 have been loaded into a Prism template as User-defined equations named "Baseline then straight line time course", "Baseline then rise-and-fall to baseline time course" and "Baseline then rise-and-fall to steady state time course" respectively (see "Signaling kinetic model-free equations" Prism file available in the supplementary files). In Prism the equations are coded with an IF statement to fit the baseline before t 0 . For example, Eq. (7) is written in Prism as, In some cases, stimulation of signaling results in a decrease of response signal. In other words, the detected response signal decreases over time. This can result from technical considerations, for example, certain biosensors decrease in signal upon binding the response analyte (so called "Downward" sensors). It can also result from the signaling mechanism, for example stimulation of G i is often detected by the resulting inhibition of cAMP accumulation. This can be handled by taking the inverse of the detected signal and plotting this versus time (creating an "Upward" time course, see for example Supplementary Fig. S2d). Alternatively, downward data can be fit to the downward analogue of the equations. Examples of this approach are used in this study (inhibition of cAMP accumulation, see Supplementary Fig. S2f, and GIRK channel currents, see Supplementary Fig. S2h); data were analyzed with the downward analogue of the association exponential equation, which is the exponential decay equation shown below (Eq. 9): Note the response at t = 0 is equal to SSR + Baseline. The corresponding equation in Prism 8 is named "Plateau followed by one phase decay, " in which Span corresponds to SSR, Plateau corresponds to Baseline, and Y0 corresponds to SSR + Baseline 96 . calculation of initial rate. The initial rate of signaling was determined by first fitting the data to the time course equations, then using the fitted parameter values to calculate the initial rate ( IR ). This calculation utilized an equation defining the initial rate in terms of the fitted parameters, the equation obtained by taking the limit of the time course equation as time approaches zero (the formal definition of the initial rate condition). These equations were derived in Appendix 1: Straight line: Association exponential: Rise-and-fall to baseline: Rise-and-fall to steady-state: Equations for kinetic mechanistic model analysis and simulation. The kinetic mechanistic model was used to simulate and analyze time course data using the model equations, derived in Ref. 33 and Appendix 2. where, fitting procedures. Data were fit by nonlinear regression to the equations using Prism 8. In almost all cases, the default fitting settings were used, specifically least-squares nonlinear regression with medium convergence criteria, no special handling of outliers and no weighting 78 . In two cases, Figs. 4b and 5b, the default settings yielded an ambiguous fit, in that some of the parameter values were not resolved (see "Curve fit results" Excel file in Supporting Material) 79 . In these cases, unambiguous parameter values were obtained using the "Robust regression" option which minimizes the contribution of outliers to the fit 74 . In the rise-and-fall fits (Eqs. 3,4,7,8), which contain two rate constants k 1 and k 2 , k 1 was assumed to be the larger of the rate constant values and this was handled by constraining k 1 to be greater than k 2 . In all cases, rate constant values were constrained to be greater than zero. The fitted values, the standard error and the correlation coefficient R 2 are given in the "Curve fit results" Excel file in Supporting Material. Values were given to four significant figures. When the signal initiation time was floated (Eqs. 5-9), the initial value (X0) was entered manually. In order to obtain an estimate of the error of the fitted parameters, the standard error was computed as described 51 , assuming a symmetrical confidence interval. This required changing the default "Calculate CI or parameters" from asymmetrical to symmetrical in the "Confidence" dialogue 97 .
The baseline response, that in the absence of ligand, was handled by incorporating a baseline term into the equations and fitting this parameter as part of the analysis. A commonly-used alternative, subtracting baseline from the data set and excluding the baseline term from the equation, was found to distort the fit significantly for the rise-and-fall to baseline equation, unless the subtracted baseline value was in very close agreement with that determined by incorporating the baseline into the curve-fitting procedure.
For the V 2 vasopressin receptor experiments, the biosensor fluorescence intensity data were normalized to baseline. Specifically, response was quantified as the fluorescence after ligand addition divided by the mean baseline signal before addition, a metric termed F/ΔF. (The mean baseline signal is the mean of the fluorescence intensity measurements before ligand addition.) The data was handled in this way because it provides a convenient intra-well control, minimizing error due to any slight well-to-well differences in the amount of sensor or number of cells. The arrestin sensor data was inverted by calculating the value 1 − F/ΔF. This was done because the arrestin sensor is a downward sensor (a decrease in fluorescence resulted from receptor interaction) 26 . In the fitting procedure, each technical replicate was considered an individual point.
Concentration-response data for the initial rate were fit to a sigmoid curve equation, the "Log(agonist) vs. response − Variable slope" equation in Prism 83 : The "Bottom" parameter was constrained to zero. (Note in the Prism formulation, L 50 is written as EC 50 . The EC 50 term is not used here because it has an explicit pharmacological definition 62 .)

Measurement of cAMp generation and arrestin recruitment using biosensors.
Geneticallyencoded biosensors were used to measure cAMP generation and arrestin recruitment via the V 2 vasopressin receptor. The sensors have been described previously (red cADDis for cAMP 82 and see ref. 26 . for arrestin recruitment), and are packaged in the BacMam vector for delivery to cells. The cDNA for the V 2 vasopressin receptor was obtained from the cDNA Resource Center (Bloomsburg University, Bloomsburg, PA). The experiments were conducted in HEK293T cells transfected with the receptor and each sensor individually (i.e. one batch of cells with receptor and the cAMP sensor and a separate batch with receptor and arrestin sensor). HEK 293T cells were cultured in Eagle's minimum essential media (EMEM) supplemented with 10% fetal bovine serum (FBS) and penicillin-streptomycin at 37 °C in 5% CO 2 . One day before the transfection, HEK293T cells were seeded at a density of 27,000 cells/well in a 96-well plate (Greiner Cellcoat #655946). Approximately 16-20 h later, cells were transfected with 50 ng of the V 2 receptor plasmid per well, using Lipofectamine 2000 from Invitrogen (Waltham, MA USA). After an approximately 4 h incubation period at 37 °C and 5% CO 2 , the transfection mix was replaced with 100uL complete cell culture media and the cells were allowed to recover for approximately 1.5 h under normal growth conditions. The cells were then transduced with either the red cADDis or green arrestin sensor BacMam stocks. To prepare the transduction mixture, the BacMam containing cADDis or the arrestin sensor, 2 mM sodium butyrate, and EMEM were combined in a final volume of 50 μL. For each experiment, 6.18 × 10 8 viral genes of cADDis virus or 4.24 × 10 8 viral genes of arrestin sensor virus were added to each well. The 50 uL www.nature.com/scientificreports/ transduction mixture was then added to the 96-well plate (50 uL/well) and incubated for approximately 24 h at 37 °C in 5% CO 2 . Thirty minutes prior to fluorescence plate reader experiments, the media in each well was replaced with 150 μL of Dulbecco's phosphate buffered saline (DPBS) supplemented with Ca 2+ (0.9 mM) and Mg 2+ (0.5 mM). Fluorescence plate reader experiments were performed on the Synergy Mx reader (BioTek, Winooski, VT). The green fluorescence detection was recorded using 488/20 nm excitation and 525/20 nm fluorescence emission, while red fluorescence detection was recorded using 565/20 nm excitation wavelength and 603/20 nm fluorescence emission. After acquiring the baseline fluorescence for several minutes, drug was added manually with a multichannel pipette in a volume of 50 µL at the indicated time points.
Oxytocin and vasopressin were obtained from Cayman Chemical (Ann Arbor, MI USA). All working concentration of drugs were dissolved in DPBS and added manually to the HEK 293T cells at the indicated concentrations and time points.
Statistical analysis. Concentration-response data for V 2 vasopressin receptor responses to vasopressin and oxytocin were compared statistically by paired, two-tailed t-test using Prism 8. Comparisons were performed for maximal effect and L 50 (Tables 1, 2).