A kinetic operational model of agonism incorporating receptor desensitization for G-protein-coupled receptors

Pharmacological responses are modulated over time by regulation of signaling mechanisms. The canonical short-term regulation mechanisms are receptor desensitization and degradation of the response. Here for the first time a pharmacological model for measuring drug parameters is developed that incorporates short-term mechanisms of regulation of signaling. The model is formulated in a manner that enables measurement of drug parameters using familiar curve fitting methods. The efficacy parameter is kτ, which is simply the initial rate of signaling before it becomes limited by regulation mechanisms. The regulation parameters are rate constants, kDES for receptor desensitization and kD for response degradation. Efficacy and regulation are separate parameters, meaning these properties can be optimized independently of one another in drug discovery. The parameters can be applied to translate in vitro findings to in vivo efficacy in terms of the magnitude and duration of drug effect. When the time course data conform to certain shapes, for example the association exponential curve, a mechanism-agnostic approach can be applied to estimate agonist efficacy, without the need to know the underlying regulatory mechanisms. The model was verified by comparison with historical data and by fitting these data to estimate the model parameters. This new model for quantifying drug activity can be broadly applied to the short-term cell signaling assays used routinely in drug discovery and to aid their translation to in vivo efficacy, facilitating the development of new therapeutics. Highlights Regulation of signaling impacts measurement of drug effect Receptor desensitization is incorporated here into a kinetic model of signaling Drug effect and signaling regulation can now be measured independently The analysis framework is designed for signaling assays used in drug discovery These new analysis capabilities will aid development of new therapeutics


Introduction
Pharmacological theory has provided the concepts and analytical framework now used routinely in modern drug discovery and development (Rang, 2006). In lead optimization, potency (EC50) and maximal effect (Emax) are used to establish structure-activity relationships and predict the in vivo effectiveness of new molecules. In pharmacological models, these empirical measures are converted into chemical drug parameters such as affinity and efficacy (Kenakin, 2009b). Most pharmacological theory and much drug development are based upon G-protein-coupled receptors (GPCRs). From bacteria to higher organisms these receptors have evolved to respond to a remarkable diversity of signals, including photons, ions, neuromodulators, hormones and enzymes (Fredriksson et al., 2003). They have proved to be highly tractable targets for small molecule therapeutics (Sriram and Insel, 2018).
GPCR signaling is regulated in order to prevent over-stimulation of cellular responses (Hausdorff et al., 1990;Krupnick and Benovic, 1998). These mechanisms act to attenuate signaling on persistent or repeated exposure of the receptor to the agonist. Regulation occurs at two levels in the signal transduction process, at the level of signal generation by the receptor (receptor desensitization), and at the level of signal degradation downstream of the receptor. In the receptor desensitization mechanism, agonist-bound receptor is phosphorylated on intracellular residues by kinase enzymes (Inglese et al., 1993;Stadel et al., 1983). The phosphorylated receptor then interacts with arrestin (Krupnick and Benovic, 1998;Lohse et al., 1990). Arrestin binding prevents access of G-protein to the receptor, blocking G-protein signaling. The signal, once generated, can be degraded by regulatory processes in the signaling pathway. These include metabolic processes that reduce the level of second messenger molecules such as cAMP and inositol phosphates 4 (Berridge, 1993;Chang, 1968;Naccarato et al., 1974), and export mechanisms such as the efflux of calcium ions (Carafoli, 1991).
Regulation of signaling is not included in the pharmacological models used routinely to measure agonist activity in drug development. [Numerous models have been developed to simulate system behavior, revealing the manifestation of mechanisms in signaling data (Leff, 1986;Riccobene et al., 1999;Violin et al., 2008;Xin et al., 2008). These models are not routinely used to measure ligand activity.] Notably, there is no pharmacological parameter for regulation of signaling in common use. Incorporating regulation of signaling would be useful for a number of reasons. It could aid prediction of in vivo activity; agonists that differ in their regulation of signaling would be anticipated to differ in their in vivo activity, especially for prolonged responses (Hothersall et al., 2016) or on repeated dosing (Tay et al., 2018). In analyzing SAR, it could tease apart effects of chemical substitutions on signaling efficacy versus receptor desensitization.
Currently, these processes are effectively combined in the response measurement and this can result in time-dependence of drug parameter estimates (Leff, 1986;Navratilova et al., 2007;Riccobene et al., 1999). This effect complicates the comparison between compounds and between signaling assays, especially in biased agonism assessment (Klein Herenbrink et al., 2016;Lane et al., 2017); in principle, differences at early time points can be largely eliminated at later time points.
A pharmacological model that incorporates regulation of signaling must incorporate time because regulation of signaling is dynamic. Regulation controls how the response evolves over time, for example fade (the response declines following a peak) and leveling off (the response reaches a plateau). Recently a dynamic pharmacological model for G-protein-coupled receptor signaling was presented for quantifying agonist efficacy in routine drug discovery pharmacology 5 . This model included one of the two regulation of signaling axes, degradation of the response. In the current study, the model is extended to incorporate receptor desensitization. This provides a model that can be broadly applied to the short-term cell signaling assays typically used in drug discovery. The aims off this study were: 1. To formulate the model 2. Determine its manifestation in the shape of time course data 3. Determine the drug parameters that can be measured by curve fitting 4. Develop a readily-adoptable data analysis framework using generic, familiar equations found in commercial curve-fitting software (e.g. GraphPad Prism, Sigmaplot, XLfit).
5. Measure drug parameters by applying the model to experimental data.

Framework
7 by phosphodiesterases (Chang, 1968), and de-activation of G-protein by hydrolysis of bound GTP to GDP by the intrinsic GTPase activity of the G-protein (Gilman, 1987). Some signals are decreased by clearance of the signaling species from the relevant compartment, for example efflux of cytosolic Ca 2+ ions (Carafoli, 1991). Response degradation is incorporated into the model as an exponential decay of the response (Fig. 1, pink region), as described previously . Signal transduction can also be dampened over time by a third mechanism, depletion of the response precursor. The most familiar example is in calcium signaling, in which depletion of Ca 2+ from intracellular stores limits the amount of Ca 2+ that can be released into the cytosol (Yu and Hinkle, 2000). This property can be incorporated into the model as a decrease of response precursor over time, as described previously .
It is well known that additional processes effect GPCR responses, especially over longer timeframes. The receptor can resensitize and so re-participate in signaling (Ferguson, 2001).
Recently, it has emerged that regulatory processes can result in prolonged signaling. For example, receptor internalization can result in a long-lived signaling complex in intracellular vesicles (Ferrandon et al., 2009;Hothersall et al., 2016). These mechanisms can be incorporated by extending the model and such models are currently being evaluated.

Formulating the model
The model was formulated in a manner that yields measurable parameters of ligand efficacy and response regulation that can be applied to drug discovery, for example to establish structure-activity relationships in lead optimization (Kenakin, 2009a).

8
In the response generation step, agonist-bound receptor (RA) interacts with a precursor of the response (EP) and converts it to the response (E). The rate of response generation is defined by mass action, i.e. is a function of the concentration of the interacting species ([RA] and EP) and a rate constant (kE). The rate is termed the transduction rate constant (kτ). It is defined as: where EP(TOT) is the total response precursor available in the system and [R]TOT the total concentration of receptors. kτ can be described as the initial rate of response generation. It is the rate of response generation before the response is appreciably affected by regulation processes. It is analogous to the initial rate of enzyme activity, [ ] [ ] , where S is substrate, E enzyme and kCAT the catalytic rate constant.
kτ is the efficacy of the agonist. It is the macroscopic efficacy of the agonist in the system being measured and as such can be viewed as a kinetic analogue of traditional efficacy parameters, such as τ of the operational model (Black and Leff, 1983). kτ is determined by the agonist because it is dependent on kE, the microscopic efficiency of the agonist for generating the response. In other words, agonists with different efficacies manifest different values of kτ for the response being measured because they differ in the value of kE. kτ is also defined by system parameters that are independent of the agonist, specifically [R]TOT and EP (TOT). As a result, the value of kτ for a single agonist can differ between response systems. This can occur for three reasons: differing receptor concentration (different [R]TOT); differing levels of response precursor (different EP(TOT); and differing types of response precursor that the receptor is coupling to (different kE).
Experimentally, kτ can be measured in routine pharmacological assays using familiar curve fitting methods and equations, as described in Section 3. A maximally-stimulating concentration of agonist can be used to measure kτ. This property enables measurement of agonist efficacy regardless of the kinetics of agonist binding to the receptor, as described below (Section 2.2.3).
Importantly, kτ can be measured independently of the rate(s) of response regulation, as described in Section 2.2.2.

Response regulation
The canonical signaling regulation mechanisms in the model are receptor desensitization and response degradation (Fig. 1). In the desensitization mechanism, the active receptor RA is desensitized, forming the inactive receptor, R0A, which cannot couple to the response precursor ( Fig. 1, yellow region). Receptor desensitization is modeled as an exponential decay of the active receptor concentration. It is assumed to be an exponential decay process because arrestin binding to the receptor has been shown to proceed according to an exponential process (for example, see (Berglund et al., 2003)). The rate is defined by the receptor desensitization rate constant kDES. The half-time for receptor desensitization is 0.693 / kDES. kDES is an agonist-dependent parameter because it is dependent on the nature of the agonist bound to the receptor. Regarding response degradation, the mechanism is assumed to be an exponential decay process, as described previously   (Fig. 1, pink region). The rate is defined by kD, the response degradation rate constant, with the half-time for degradation being 0.693 / kD. This rate is independent of the agonist because response degradation occurs downstream of the receptor.
It is critical to note that the agonist efficacy parameter kτ is independent of the system regulation parameters (kDES and kD) in this model. This means that agonist efficacy can be 10 measured in a system in which the response is regulated, as shown in Section 3. This emerges from the microscopic features of the model. At the level of a single receptor, the process of response generation is independent from the processes of receptor desensitization and response degradation.
Mathematically, the capacity of a single agonist-bound receptor to generate the response is governed by kE (see green panel in Fig. 1), not by kDES or kD. This is because mechanistically: 1) receptor desensitization only affects the number of active receptors, not the capacity of a single active receptor to generate the response. 2) Response degradation (kD) does not affect response generation because it is downstream of the receptor. From a macroscopic perspective, the efficacy parameter measured in the experiment, kτ, is independent of the regulation parameters for two reasons: 1) It is the initial rate of signal transduction, i.e. the rate before the response is affected by the regulation processes.

Agonist binding
In classical pharmacological models, it is implicitly assumed that agonist-receptor binding is at equilibrium at the time point at which the response is measured (Bdioui et al., 2018;Black and Leff, 1983;Colquhoun, 1998;Rang, 2006;Slack and Hall, 2012). In a kinetic model, it cannot be generally assumed that agonist is at equilibrium with the receptor at all time points of response measurement (unless a maximally-effective concentration of agonist is employed -see below).
The equilibrium assumption is probably reasonable for low potency ligands when the response progresses over timescales of at least several minutes, for example second messenger generation and gene expression (Bdioui et al., 2018). Such ligands (EC50 in the micromolar range) equilibrate rapidly even at the low concentrations required to span the concentration-response curve 11 (equilibration time < 1 min) . However, high potency agonists (EC50 in the low nM range) equilibrate more slowly, especially at the low concentrations required to span the concentration-response curve (at least several minutes). In addition, lack of equilibration is likely to be an issue for responses that are extremely rapid (Bdioui et al., 2018;Charlton and Vauquelin, 2010). For example, the rise of intracellular Ca 2+ occurs within a few seconds of agonist application. This problem can be solved by incorporating the kinetics of agonist binding into the model. This approach is applied for the new models in this study (see Appendix A, Supplementary   Fig. S2, Supplemental information). This approach was applied to the original kinetic operational model  and a variant of the operational model that assumes agonist binding is rate-limiting (Slack and Hall, 2012).
Fortunately, agonist efficacy can be determined in a way that minimizes the equilibration problem. This is by virtue of the fact that a saturating concentration of agonist can be used to measure agonist efficacy (see Section 3). The higher the concentration of ligand, the more rapidly it associates with the receptor, due to mass action. For example, let us consider a slowlyequilibrating agonist with an equilibrium dissociation constant of 10 nM and a t1/2 for association, at 10 nM, of 10 min. If this agonist is applied at the maximal concentration typically employed in in vitro signaling assays (10 µM), association with the receptor is greatly accelerated (association t1/2 of 1 second).

Equation format
In a kinetic functional assay, the response to agonist is measured over multiple time points and the response value is plotted against time. Equations were derived for analyzing response vs time data to obtain estimates of the model parameters (Appendix A). The equations are of a form that can be used by curve fitting software commonly-employed in drug discovery and receptor research, for example Prism (GraphPad Software, Inc.) XFfit (ID Business Solutions Ltd.) and SigmaPlot (Systat Software, Inc.). These equations take the analytic form, y = f(t), where y is response and f(t) is a function of time and pharmacological parameters. This is referred to as an E vs t equation in this study. [Note that certain programs allow fitting of data to differential equations rather than their analytic solutions, such as Dynafit from BioKin Ltd, (Kuzmic, 2009).] Some of the models resulted in systems of nonlinear differential equations, which are often difficult to solve analytically. Surprisingly, analytical solutions were obtained for most of these systems. Importantly, this allowed extension of the model to incorporate agonist binding kinetics (for example, Appendix A.5) and to accommodate simultaneous receptor desensitization and precursor depletion (Appendix A.9). 13

Results
The manifestation of the model in response time course data is first evaluated using simulated data generated using the model equations. The models are then verified by comparison with historical experimental data. In these historical studies, regulation mechanisms were manipulated (effectively deleted), enabling a piecewise verification of the model. The model is then applied to the experimental data to obtain estimates of the model parameters, using standard curve-fitting procedures. In the course of this investigation, a pattern emerged in which the shape of the time course was determined by the number of regulation processes. No regulation processes results in a straight line (Section 3.1). A single regulation process results in an "association" exponential curve, here termed a horizontal exponential curve (Section 3.2). Two processes (an input process and an output process) results in a rise-and-fall exponential curve (Section 3.3). This property enabled a generalized approach to be used to analyze these data (Sections 3.2.4 and 3.3.4).
In this approach, agonist efficacy (kτ) could be measured without knowledge of the underlying regulation mechanism. More complex curves resulted from multiple mechanisms acting in concert.
The models were considered in order of increasing complexity: 1. An unregulated response, giving a straight line E vs t profile (Section 3.1).
2. A single regulation mechanism, giving a horizontal exponential curve (Section 3.2).
3. Receptor desensitization or precursor depletion (input process), and response degradation (output process), resulting in a rise-and-fall curve (Section 3.3).
4. More complex mechanisms resulting in more complex curve shapes (Section 3.4).
14 Some of the models were described in the original study -Models 1,3,4 and 8 .
Here they are considered again in the context of regulation of signaling and in order to enable comparison with the new models.
In this section the models are considered for a maximally-stimulating concentration of agonist, which enables measurement of efficacy (kτ) and the regulation parameters. Data for multiple concentrations of agonist are given in Supplementary information.

Unregulated response -straight line profile (Model 1)
In the absence of any regulation mechanism, the model predicts a linear increase of the response over time. In other words the E vs t profile is a straight line. This makes sense because there is nothing to stop the signal being generated, and there is no degradation of the signal, so the response accumulates continuously. This mechanism is shown schematically in Fig. 1 (greenshaded region) and in Scheme 1 below:

Scheme 1
This model has been described previously, being the minimal mechanism of the kinetic operational model . Response is generated from response precursor, at a rate governed by  (Naccarato et al., 1974). In this response there is minimal precursor depletion because the precursor, phosphatidylinositol 4,5-bisphosphate, is continuously generated (Xu et al., 2003). Consequently, an inositol phosphates assay run in the presence of Li + can be reasonably assumed to be minimally impacted by two of the regulation mechanisms, response degradation and precursor depletion. This leaves us with the remaining mechanism, receptor desensitization.
Blocking receptor desensitization in systems with minimal degradation and depletion results in a linear time course of the response, consistent with the model. This is shown in the following examples. Receptor desensitization in the inositol phosphates pathway has been blocked by genetic ablation of β-arrestin. Embryonic fibroblasts were isolated from arrestin knock-out mice and GPCR signaling dynamics evaluated in these cells (Kohout et al., 2001). In this system, the IP time course was a straight line. This was shown for angiotensin II via the AT1 angiotensin receptor (Kohout et al., 2001), and for thrombin via the proteinase-activated PAR1 receptor (Paing et al., 2002). The AT1 receptor data are reproduced in Fig. 3. A second example of diminished receptor desensitization is provided by the mammalian gonadotropin releasing hormone GnRH1 receptor.
This receptor lacks the C-terminal tail (Eidne et al., 1992). This region of the receptor is a critical determinant of canonical receptor desensitization because it contains the sites of receptor phosphorylation and is a determinant of arrestin binding (Krupnick and Benovic, 1998). A linear time course of inositol phosphates generation has been observed for the mammalian GnRH1 receptor expressed endogenously and heterologously in numerous cell lines (Davidson et al., 1994;Davis et al., 1986;Heding et al., 1998). (Data from two of these studies are reproduced in Figs. 4 and 5.) [Longer durations of GnRH exposure (typically > 20 min) results in attenuation of the response, owing to arrestin-independent receptor internalization (Heding et al., 2000).] The kτ value can be estimated by fitting experimental data to Eq. 1. This can be done by fitting the E vs t data for maximally-stimulating concentration of agonist to a generic straight line equation: = where m is the gradient of the line. kτ is equal to m (evident from Eq. 1). This analysis is applied here to the arrestin knock out and GnRH1 receptor data . For angiotensin II-stimulated inositol phosphates accumulation via the AT1 receptor in arrestin knockout mouse embryonic fibroblasts (Kohout et al., 2001), the kτ value is 0.13-fold over basal.min -1 (Fig 3). For GnRH-stimulated inositol phosphates accumulation via the GnRH1 receptor in HEK293 cells (Heding et al., 1998), the kτ value is 2.2% of total added radioactivity.min -1 (Fig. 4). For GnRH-stimulated IP3 accumulation via the GnRH1 receptor endogenously-expressed in rat granulosa cells (Davis et al., 1986), the kτ value is 580 cpm IP3.min -1 (Fig. 5).

Single regulation mechanism -horizontal exponential profile
In this section models are considered in which a single regulation mechanism is in operation -receptor desensitization, response degradation or depletion of response precursor. In all cases the E vs t profile is a horizontal exponential curve. (The horizontal exponential curve is defined by the equation = Plateau × (1 − − . ). A familiar example is the ligand binding association curve.) This commonality enables a mechanism-agnostic approach to be applied to analyze horizontal exponential time course data (Section 3.2.4.).

Receptor desensitization (Model 2)
Receptor desensitization is represented in the model as a reduction of the active receptor concentration over time ("active" meaning the capacity to generate a response). The basic model of receptor desensitization is shown schematically in Fig. 1 (green and yellow shaded regions) and in Scheme 2 below:

Scheme 2
The time course profile predicted by this model is a horizontal exponential curve ( Fig 2B).
This 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 horizontal plateau). At this limit the response level does not change because no new response is being generated and the existing response is not degraded. The equation defining response over time for the desensitization model, for a maximally-stimulating concentration of agonist, is Eq. 2, derived in Appendix A.2: The effect of receptor desensitization on the shape of the time course can be tested by introducing receptor desensitization into a system that is unregulated. This can be done using the systems described in Section 3.1 (Figs. 3 and 4). In the arrestin knock out cell system the effect of desensitization can be tested by comparing wild-type cells that express β-arrestin with the knockout cells that do not. For the AT1 inositol phosphate response (Kohout et al., 2001), the shape of the E vs t profile is a horizontal exponential curve when arrestin is expressed (Fig. 3), implying receptor desensitization results in a horizontal exponential curve. The same finding has been reported for the PAR1 receptor (Paing et al., 2002). For the GnRH1 receptor system, receptor desensitization has been introduced by incorporating a C-terminal tail into the receptor, from a receptor that is desensitized (the thyrotropin-releasing hormone TRH1 receptor) (Heding et al., 1998). This results in a horizontal exponential curve (Fig. 4). The horizontal exponential time course in these receptor desensitization systems provides experimental verification for the desensitization model.
The model parameters can be estimated by fitting experimental time course data to Eq. 2.
A simple way to do this is to analyze the data using a generic horizontal exponential equation.
First, the E vs t data for a maximally-stimulating concentration of agonist is fitted to the following generic horizontal exponential equation: "Plateau" is the response at the horizontal asymptote (i.e. E as t → ∞) and K is the observed rate constant. This can be done in Prism using the built-in "One phase association" equation (Motulsky, 2019a (Fig. 3).
Multiplying Plateau by K gives the kτ value 0.10-fold.min -1 . If this analysis is valid, the kτ value should be equal to the value in the non-desensitized system (arrestin knock-out cells), assuming the only difference between the systems is receptor desensitization. The data are in reasonable agreement with this assumption; the kτ value in arrestin knockout cells was 0.13-fold.min -1 (Fig.   3). The kDES value is equal to K. The kDES value was 0.042 min -1 . This translates to a half-time of desensitization of 17 min for the angiotensin II-occupied AT1 receptor in this system. For the GnRH1 receptor engineered with a C-terminal tail, the Plateau value was 1.8 % and the K value was 0.86 min -1 (Fig 4). The calculated kτ value is 1.5 %.min -1 , close to the value for the nondesensitized receptor of (2.2 %.min -1 ). The kDES value, equal to K, was 0.86 min -1 . This translates to a half-time of desensitization of 0.81 min for the C-terminally-extended GnRH1 receptor in this system.

Response degradation (Models 3 and 4)
Response degradation was incorporated into the model as an exponential decay of the response, as described previously . The model, Model 3, is shown schematically in Fig. 1 (green and pink shaded regions) and in Scheme 3 below:

Scheme 3
The E vs t profile predicted by this model is a horizontal exponential curve (Fig. 2C). This can be rationalized as follows. At first, the response level increases rapidly because the response generation process predominates. As time progresses, the rate of increase slows, as response degradation increases. Finally, a plateau of the response approached, at which there is a steadystate where the rate of response generation equals the rate of response degradation. The equation for this model, at a maximally-stimulating concentration of agonist, is Eq 8, derived in Appendix A.3: Eq. 3

22
A variant of this model incorporates recycling of the response back to the response precursor . It is shown in Scheme 4 below:

Scheme 4
The equation The effect of response degradation on the shape of the time course can be investigated by leaving out the response degradation inhibitor. For the inositol phosphates response this can be done by omitting Li + , the inhibitor of the terminal step of inositol dephosphorylation (Naccarato et al., 1974). Using the GnRH1 receptor enables degradation to be examined in isolation from receptor desensitization because this receptor does not desensitize acutely (Section 3.1). Data from this system are reproduced in Fig. 5 (Davis et al., 1986). Omitting Li + results in a horizontal exponential curve (Fig. 5), implying response degradation results in a horizontal exponential curve.
This finding is consistent with the model.
A limitation of the degradation inhibitor approach is noted here. The most familiar example of response degradation is breakdown of cAMP by phosphodiesterase (PDE) enzymes. This process is inhibited, and cAMP elevated, by PDE inhibitors. This treatment is used routinely in measuring cAMP signaling mediated by GPCRs. However, standard methods do not completely inhibit PDE activity, such that some degradation of cAMP remains. This is because of solubility limits of the most-commonly used PDE inhibitor, isobutylmethylxanthine. A concentration commonly applied is 0.5 mM, which inhibits degradation by only 75% (Schulz and Mailman, 1984). Subtype-selectivity is an issue with another commonly used inhibitor, rolipram. This compound is selective for PDE4 (Houslay et al., 2005) but cell lines used routinely for the study of GPCRs express multiple PDE subtypes and appreciable cAMP levels can be detected in the presence of subtype-selective PDE inhibitors (Motte et al., 2017). Incomplete inhibition of cAMP breakdown can complicate the interpretation of the effect of response degradation on cAMP signaling (Gimenez et al., 2015;Xin et al., 2008). In the context of the kinetic model in this section, 24 application of a response degradation inhibitor at a partially effective concentration does not linearize the E vs t profile but instead results in a decrease of the rate constant and an increase of the plateau (data not shown).
The model parameters can be estimated by fitting experimental time course data for a maximally-stimulating agonist concentration to the generic horizontal exponential equation: kτ can be determined by multiplying the Plateau by K, as is evident from Eqs. 8 and 11. This method can be applied whether or not the response is recycled, evident from inspecting the equations. This analysis was applied to the GnRH signaling data in Fig. 5. The Plateau value was 2,100 cpm IP3 and the K value 0.12 min -1 . The resulting value of kτ was 250 cpm IP3.min -1 . This value was in the same range as that for the unregulated response (580 cpm IP3.min -1 , for the response in the presence of Li + , Section 3.1).
It is not possible to determine kD using this approach unless it is known whether or not response is recycled, i.e. whether Scheme 4 or Scheme 5 applies. This is because the definition of K is different in these two scenarios. If response is not recycled, K is equal to kD. If response is recycled, K is equal to kτ / EP(TOT) + kD.

Depletion of response precursor (Model 5)
In this mechanism, the only process of response limitation is depletion of response precursor. The precursor is depleted because it is converted to the response by the response-25 generation process. This mechanism results in a horizontal exponential profile. This model, Model 5, is rarely encountered and is presented here for the purpose of completeness.
This model is represented by Scheme 5:

Scheme 5
The equation defining this mechanism, for a maximally-stimulating concentration of agonist is Eq. 5, derived in Appendix A.5: This is a horizontal exponential equation. Analysis using the generic horizontal exponential equation (Section 3.2.1) gives a Plateau value corresponding to EP(TOT) (the total amount of response precursor) and a K value corresponding to kτ / EP(TOT). kτ can be determined by multiplying the Plateau value by the K value, as is done for the other horizontal exponential models. An example from the literature can be found in early studies of glucagon-stimulated cAMP accumulation in hepatic membranes (Pohl et al., 1971). This system is minimally regulated, lacking receptor desensitization machinery owing to the use of washed membranes, and lacking phosphodiesterase activity (Pohl et al., 1969). Depletion of response precursor results from the absence of an ATP-regenerating system, ATP being the precursor of cAMP. The resulting E vs t profile is a horizontal exponential curve (see Fig. 4B of (Pohl et al., 1971)), consistent with the kinetic model.

Mechanism-agnostic analysis of horizontal exponential time course data for estimating agonist efficacy and affinity
Time course data from all of the horizontal exponential models considered here can be analyzed the same way to determine agonist efficacy for generating the response, i.e. kτ. This is done by multiplying the Plateau by the rate constant of the horizontal exponential fit for a maximally-stimulating concentration of agonist. The theoretical basis of this approach is evident from the general form of all the model equations (Eqs. 5,8,11 and 14): where K is the observed rate constant. From this equation, it is evident that kτ can be determined by multiplying the Plateau by the rate constant: This means that agonist efficacy (kτ) can be determined without knowing the regulatory mechanism underlying the horizontal exponential profile, assuming the profile results from one of the models described in this study. This analysis method can be described as agnostic regarding the regulatory mechanism. The observed rate constant K can be described as the response regulation rate constant, here termed kREG. This rate constant governs the time it takes for the response to level off.
This agnostic approach is used here to re-analyze historical data on agonist activity at the µ-opioid receptor, from the study of . In this study, the time course data for the response, [ 35 S]-GTPγS binding, conform to a horizontal exponential curve and the data were analyzed in the original study using the generic equation to determine the plateau and the rate constant (termed Bmax and k, respectively, in the original article) . The parameter values from the study are shown in Table 1. From these values, kτ can be calculated, by multiplying Bmax by k. The kτ values indicate a broad range of agonist activity, spanning 20-fold, from 0.0024 fraction.min -1 for nalbuphine to 0.048 fraction.min -1 for DAMGO. The results are consistent with the literature, since: 1) Butorphanol and nalbuphine, partial agonists in terms of the kτ value, are partial agonists in other assays of µ opioid receptor activity (Hoskin and Hanks, 1991). 2) Morphine has a slightly lower kτ value than DAMGO (Table 1), and has slightly lower Emax than DAMGO in [ 35 S]-GTPγS binding assays (Traynor and Nahorski, 1995).
Agonist affinity can also be determined using the agnostic approach. This analysis assumes agonist rapidly equilibrates with the receptor. [If it does not, the E vs t profile deviates progressively from a horizontal exponential curve as the agonist concentration is decreased 28 (Supplemental Fig. S2).] This analysis is based on the common general form of all the horizontal exponential model equations for non-saturating concentrations of agonist (Eqs. 4, 7, 10 and 13): Multiplying the Plateau by the rate constant gives kτ: This can be done in Prism using the "log(agonist) vs. response --Variable slope (four parameters)" equation, with the parameter "Bottom" set to zero (Motulsky, 2019b). The A50 value is KA and the y value at maximally-effective agonist concentrations ("Top" in the equation above) is kτ. The Hill slope parameter accommodates effects that deviate the concentration range exposed to the receptor from the concentration range added to the assay, for example precision of serial dilution or biological processes that affect local concentration in the vicinity of the receptor. The Hill slope can also allow cooperativity of agonist binding to be incorporated.

Two regulation mechanisms: Rise-and-fall exponential profiles
In unmodified response systems, GPCR signaling is typically regulated at both levels of the regulation system, the response generation process (receptor desensitization or, in the case of calcium signaling, precursor depletion) and the response degradation process. This scenario results in a rise-and-fall time course profile, commonly observed for unmodified GPCR responses, as shown below. The kτ value can be determined by fitting time course data to a generic equation and this can be done without knowing the mechanism that regulates the response generation process (Section 3.3.3). The mechanisms that result in this profile are presented below. They incorporate response degradation with a second regulation process operating on the response input, either receptor desensitization or precursor depletion.

Receptor desensitization and response degradation (Model 6)
First, the model with receptor desensitization and response degradation is considered. The model combines the yellow, green and pink regions of Fig. 1 and is shown schematically in Scheme 6 below:

Scheme 6
This model predicts a rise-and-fall to zero curve (Fig. 2D). 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. The response level declines because it is being degraded and ultimately the response level falls to zero once all existing response has been degraded.
The equation defining the response over time is, for a maximally-stimulating concentration of agonist, Eq. 6 (derived in Appendix A.6): This equation is of the same form as a generic rise-and-fall to zero exponential equation, familiar in pharmacokinetics as the equation defining drug concentration over time after oral dosing: The model was validated by comparison with experimental data. Fig. 6 shows data for the IP3 response to cholecystokinin (CCK) via the CCK1 receptor in pancreatic acinar cells (Streb et al., 1985). The mechanism of this response conforms to Scheme 6: It is regulated by response degradation [shown by sensitivity to Li + , Fig. 6, (Streb et al., 1985)] and by receptor desensitization (Klueppelberg et al., 1991). The shape of the time course is a rise-and-fall to zero exponential curve (solid squares, Fig. 6), consistent with the kinetic model. A second example is diacylglycerol signaling by the AT1 angiotensin receptor (Violin et al., 2006) (Fig. 7). This response is modulated by receptor desensitization; it is sensitive to removal of phosphorylation sites from the receptor (Violin et al., 2006) (Fig. 7). It is presumably also modulated by response degradation because diacylglycerol is cleared by diacylglycerol kinases (Merida et al., 2008). The shape of the time course is a rise-and-fall to zero exponential curve (Fig. 7, sold squares), consistent with the kinetic model.
The data in Fig. 6 and 7 can be analyzed using the model to determine agonist efficacy for generating the response, in other words the kτ value. A simple way to do this is to fit the data to the generic rise-and-fall to zero equation: By comparison with Eq. 6, it is evident that C is equal to kτ. This method is applied to the experimental data in Figs. 6 and 7. (Note the analysis assumes a maximally-stimulating agonist concentration, employed in these experiments. (0.14 min -1 ). Therefore, in this case, K2 was ascribed to kDES. This means, by a process of elimination, K1 corresponds to kD. Consequently, kD was estimated to be 0.39 min -1 . In simple terms, this means the half-life for CCK1-receptor desensitization was 5.0 min and the half-life for IP3 degradation was 1.8 min.
The same approach was applied to the AT1 diacylglycerol response. In this case, receptor desensitization was blocked using a phosphorylation-deficient receptor (Fig. 7). This results in a horizontal exponential curve (Fig. 7), as anticipated for a response regulated solely by response degradation (Section 3.2.2.). Analysing these data as described in Section 3.2.2. gave a kD value of 6.1 min -1 . This value was close to that of K1 of the wild-type response (4.1 min -1 ) so K1 was ascribed to kD. As a result, by a process of elimination K2 was assigned to kDES, giving a kDES value of 0.74 min -1 . From this analysis it is concluded the half-life for AT1 receptor desensitization in this system was 0.94 min and the half-life of diacylglycerol clearance was 0.17 min.

Precursor depletion and response degradation -calcium signaling (Model 8)
In this model, the response precursor concentration is depleted by being converted to the response, and the response, once generated, is completely cleared. The model has been presented previously . It is represented by Scheme 8 below:

Scheme 8
This model predicts a rise-and-fall to zero curve (Fig 8). This makes sense intuitively. The response rises rapidly at first as response generation predominates. Then response generation becomes limited because the precursor becomes depleted, and response is further limited because it degrades. A peak is reached at which the rate of response generation equals the rate of response degradation. Ultimately all the accessible precursor becomes depleted so no new response is generated. The response ultimately declines to zero because the existing response is completely degraded and no new response is generated because the precursor is depleted.
This model can describe calcium signaling via GPCRs in certain circumstances, specifically Ca 2+ entry into the cytoplasm from the endoplasmic reticulum (ER). The response precursor is Ca 2+ stored in the ER. Response generation is the increase of cytoplasmic Ca 2+ concentration resulting from opening of the IP3-gated Ca 2+ channel in the ER membrane (Berridge, 1993;Berridge et al., 2003). Precursor depletion is the decrease of Ca 2+ in the ER resulting from Ca 2+ release into the cytoplasm (Hofer et al., 1998;Miyawaki et al., 1997;Yu and Hinkle, 1997;Yu and Hinkle, 2000). Response degradation is the clearance of cytoplasmic Ca 2+ into the extracellular space by the plasma membrane Ca 2+ pump (Carafoli, 1991;Hofer et al., 1998).
Receptor desensitization likely does not greatly impact the temporal profile of a single Ca 2+ response for rapid Ca 2+ responses; receptor desensitization proceeds over minutes, whereas a rapid Ca 2+ response peaks in a few seconds and is largely cleared within a minute (exemplified in Fig.   8). (Receptor desensitization likely operates to limit the Ca 2+ response on a second application of the agonist, or to limit sustained Ca 2+ responses.) The model does not accommodate the refilling of the ER with Ca 2 (which would represent response recycling). ER refilling can be blocked by removing extracellular calcium, since the ER is refilled primarily by Ca 2+ that enters into the cell (Hofer et al., 1998;Tsien and Tsien, 1990;Yu and Hinkle, 2000).
The equation defining the response over time is, for a maximally-stimulating concentration of agonist, Eq. 8, derived in Appendix A.8: This equation introduces a parameter that represents the rate of precursor depletion. This parameter is the depletion rate constant, termed kDEP. It is defined as the product of the receptor concentration and the response generation rate constant: The half-life of response precursor is determined from this parameter, as 0.693 / kDEP.

36
The model was applied to experimental Ca 2+ mobilization data to estimate agonist efficacy and the rate constants. This was done using data for the GnRH1 receptor (Merelli et al., 1992). [In this case, we can assume the response is unaffected by receptor desensitization because this receptor does not undergo short-term desensitization (Section 3.2.1).] This can be done by fitting the data to the generic rise-and-fall equation introduced in Section 3.3.1: Agonist efficacy (kτ) can be determined from this fit -it is equal to C. Fitting the equation to Ca 2+ mobilized by the GnRH1 receptor in α3T-1 clonal pituitary gonadotrophs (Merelli et al., 1992) gives a kτ value of 290 nM.sec -1 (Fig 8). The fit also provides estimates of K1 and K2. In principle it isn't possible to determine which of these parameters can be ascribed to which of the model rate constants (kDEP or kD) (see Section 3.3.1.). However, in practice it is likely that the faster of the two fitted rate values is kDEP, because this is the rate of emptying of the Ca 2+ store, which presumably proceeds more rapidly than clearance of Ca 2+ from the cytoplasm (defined by kD).
Applying this logic to the fit to the GnRH data gives a kDEP value of 0.49 sec -1 and kD value of 0.071 sec -1 (Fig. 8), corresponding to half-times for store emptying and cytoplasmic clearance of 1.4 sec and 9.8 sec, respectively.

Mechanism-agnostic analysis of rise-and-fall time course data for estimating agonist efficacy and affinity
37 Data from the rise-and-fall mechanisms in this study can be analyzed using a generic equation to estimate agonist efficacy (kτ). When using a maximally-stimulating concentration of agonist, the data are fit to the equation below and the fitted parameter C is equal to kτ: A Prism template is provided in the Supplementary material in which this equation, termed " [Pharmechanics] Rise-and-fall to zero time course" is pre-loaded.
Agonist affinity, KA, can also be determined using this agnostic analysis if the agonist equilibrates rapidly with the receptor. [If it does not, the E vs t profile deviates progressively from a horizontal exponential curve as the agonist concentration is decreased (Supplemental Fig. S2).] For non-saturating concentrations of agonist the equations conform to the following general equation:

Other mechanisms
More complex mechanisms involve multiple regulation processes acting in concert. These mechanisms result in equations that are more complex than for the mechanisms described thus far.

Receptor desensitization and response degradation with recycling (Model 7)
This model is an extension of Model 6. The receptor desensitizes and the response degrades but in this case the response degrades back to the response precursor, i.e. it recycles. The model is represented by Scheme 7 below:

Scheme 7
39 An E vs t equation could not be found for this model (see Appendix A.7. for details). The E vs t profile was simulated by numerical solution of the differential equations and is shown in Fig. 9A. The profile is a modified form of the rise-and-fall profile in which there is a lag in the decline phase. This results in a bulge on the right-hand, downward part of the curve (Fig. 9A). This results because the recycling regenerates the precursor, delaying the decline of the response that results from degradation of the response.

Receptor desensitization and response precursor depletion (Model 9)
In this model, the receptor desensitizes and the precursor is depleted as it is being converted to the response. The model is represented by Scheme 9:

Scheme 9
40 This model is included in part to demonstrate that analytical solutions can be obtained even though the system of differential equations is nonlinear. [Another example is the slow agonist equilibration variant of Model 5 (Appendix A.5.).] This is because the equation defining E contains the product of two time-dependent terms, EP and [RA]:

= [ ]
The derivation method is shown in Appendix A.9. and the E vs t equation for a maximallystimulating concentration of agonist is Eq. 9: This is the equation for a Gompertz curve, originally derived to describe human mortality and then applied generally to modeling population biology (Kirkwood, 2015). It resembles the horizontal exponential equation but the resulting curve ascends more slowly at the later time points. This is shown in the simulated E vs t data in Fig. 9B.

41
In this model all three mechanisms are in operation. The model is represented by Scheme 10 below:

Scheme 10
An analytical equation for response could not be found [although a solution for response precursor could be found (Eq. A.14, Appendix A.10, Fig. S3B)]. E vs t data were simulated by numerical solution of the differential equations (Appendix A.10.). A rise-and-fall type profile results (Fig. 9C). It is evident by visual comparison that this profile appears similar in shape to rise-and-fall profiles fitted by the generic rise-and-fall equation, for example Fig 2D). Even though mathematically the model does not conform to the generic rise-and-fall equation, simulated data generated using the model were very well, although not precisely, fit by the generic equation (Fig.   9C). This makes sense intuitively because the effects balance out. The more depletion occurs, the less influence receptor desensitization has on limiting the response, and vice-versa. Conceptually the model contains two macroscopic regulation processes -regulation of the input that generates 42 the response (receptor desensitization and precursor depletion), and regulation of the output (response degradation). We ran data simulations with a diverse array of model parameter values for a maximally-stimulating concentration of agonist and fit the data to the generic rise-and-fall equation in Section 3.3.3. In all cases the fitted kτ value was within 15% of the simulation value (data not shown).

Discussion
In this study a kinetic pharmacological model of GPCR signaling is presented that incorporates regulation of signaling. Regulation mechanisms control the magnitude of the signal over time and so impact the measurement of pharmacological activity of agonist ligands. The model incorporates the two canonical mechanisms of short-term regulation, receptor desensitization and response degradation, providing a comprehensive model that can be broadly applied to signaling assays used in drug discovery. The model was verified using historical data from studies in which the regulation processes were carefully manipulated. Meaningful drug parameters were obtained using familiar curve fitting methods applied to time course data. The parameters quantify the initial rate of signal generation by the agonist (kτ, a measure of efficacy) and the rate(s) of signal regulation. These parameters are designed for determining structureactivity relationships of agonist efficacy and signaling regulation, especially in the leadoptimization stage of drug discovery.
The drug efficacy term that emerges from the model is kτ . This is simply the initial rate of signaling stimulated by the agonist, defined as the product of the precursor concentration, receptor concentration and a rate constant (Section 2.2.1). As an efficacy parameter, kτ has a number of benefits. It is biochemically meaningful and, being analogous to the initial rate of enzyme activity, conceptually familiar to most investigators. From a pharmacological theory perspective, kτ has the benefit of being independent of receptor reserve, as described previously . Consequently, it might be a better predictor in in vivo efficacy than efficacy measurements that are dependent on receptor reserve. kτ is more straightforward to measure than other efficacy parameters that take into account receptor reserve. For example, measuring τ of the operational model usually requires system manipulations and/or grouped analyses (Black and Leff, 44 1983;Kenakin et al., 2012;Klein Herenbrink et al., 2016;Slack and Hall, 2012;Van der Graaf and Stam, 1999;Zernig et al., 1996). All that is required for measuring kτ is a time course for a maximally-stimulating concentration of agonist (see below).
The regulation parameters in the model are the rate constants for receptor desensitization (kDES) and response degradation (kD value is simply the fitted value of the rate constant. Remarkably, it is not necessary to know the regulatory mechanism in order to determine kτ, providing the data conform to a horizontal exponential or the rise-and-fall curve. This means kτ can be determined using a mechanismagnostic analysis; the data are fit to the generic equation and kτ determined from the fitted generic parameters. These fits also provide an agnostic regulation parameter, kREG. This parameter defines the rate at which the signal is dampened by the regulatory process. If the regulatory mechanism is subsequently determined, the generic rate constants can be ascribed to kDES, kD or kDEP. The notable limitation of the model is the technical difficulty of measuring equilibrium binding affinity of the agonist for the receptor. It cannot be generally assumed that agonist binding is at equilibrium with the receptor throughout the time course of response measurement. It can take time for the concentration of the agonist-receptor complex to rise and closely approach its equilibrium value (Section 2.2.3). While likely not an issue for low potency ligands (affinity in the µM range), lack of equilibration is likely to affect the time course of response for high affinity ligands, leading to erroneous estimates of ligand affinity. This problem is exacerbated for rapid responses, such as calcium mobilization (Bdioui et al., 2018;Charlton and Vauquelin, 2010). In principal, the problem can be avoided by incorporating receptor binding kinetics of the agonist into the model equations. Such equations were derived for most of the models (Appendix A).
Simulated data indicated a lag in the initial rise of response when agonist equilibrates slowly with receptor (Supplemental Fig. S2). This feature is potentially diagnostic of slow agonist equilibration. Slow equilibration can also be detected in agonist washout or blockade experiments . These considerations notwithstanding, the quantification of agonist binding kinetics using the functional models here is likely to be challenging because the equations are highly parameterized.
In drug discovery, a simple application of the model is as follows. A biosensor is used to measure response over time. Biosensors enable real-time measurement of signaling using a single well/plate for all the time points, greatly improving the efficiency time course measurement (Lohse et al., 2012;Marullo and Bouvier, 2007 the efficacy value; kREG, the regulation parameter; and the conventional EC50 and Emax parameters measured at a single time point. The data from such a paradigm could be used as follows in a lead optimization project. As done conventionally, the EC50 at a single time point gives an estimate of potency, albeit one affected by receptor reserve and potentially affected by differences of efficacy.
Efficacy can now be assessed using kτ. This provides a better estimate of efficacy than Emax from a concentration-response curve at a single time point because it is independent of receptor reserve.
This potentially aids selection of compounds for in vivo testing, and rationalization of results from these studies. The regulation rate constant might also be useful translationally. Since it defines the 47 rate of blunting of the response, it might contribute to defining the duration of efficacy in vivo, depending on the relevance of the regulation mechanisms in the signaling assay to the duration of the response in vivo.
In summary, this study reports a new pharmacological model for measuring drug parameters in conventional signaling assays that, for the first time, considers regulation of signaling. The model will facilitate the development of new therapeutics by enabling distinction of agonist efficacy from regulation in lead optimization and by improving the translation of pharmacological data from in vitro assays to in vivo efficacy. 48

Appendix A: Derivation of equations
In this study, equations were derived that define the response over time for the kinetic models. The equations are of the analytic form, y = f(t), where y is response and f(t) is a function of time and pharmacological parameters. The analytic form is accommodated by curve-fitting software commonly used by pharmacologists. The model equations were derived as follows. First, the model mechanism was drawn out schematically in chemical terms (see Fig. 1 for general scheme). Next, differential equations describing the time-dependent processes were derived based on first principles, specifically the law of mass action for interactions between species (e.g. response precursor and receptor), and exponential decay for the decline of a single species over time (e.g. response degradation). Next, the differential equations were re-arranged to enable straightforward solution to the analytic form. In some cases, this step involved the Laplace transform method (Hoare, 2017;Mayersohn and Gibaldi, 1970). Finally, the analytic solution was obtained. The derivation assumes the RA-EP complex is sufficiently transient that it does not contribute appreciably to the mass balance of receptor or response precursor.
The equations were used to simulate model behavior and analyze experimental data (Section 3). In addition, a side-by-side comparison of E vs t data from the models is presented in Supplemental information. The following models have been described previously -1,3,4 and 8 .

A.1. Unregulated response (Model 1)
In this model the response is not regulated -the receptor does not desensitize, the response does not degrade and there is no precursor depletion. This is Model 1 of the original kinetic model . The model is represented by Scheme 1:

Scheme 1
The E vs t equation for a rapidly-equilibrating agonist was derived in  At maximally-effective agonist concentrations, fractional occupancy is unity, i.e., The E vs t equation for a saturating concentration of agonist is then (Eq. 1), ,[ ]→∞ =

Eq. 1
Note that this is the equation for a straight line intersecting the origin, where the gradient is kτ.
An equation has been derived that assumes agonist binding is not at equilibrium with the receptor [Eq. 9 of ]. Here k1 is the agonist-receptor association rate constant and k2 the dissociation rate constant. The equation is,

A.2. Receptor desensitization (Model 2)
Receptor desensitization is incorporated into the model as a decline of the agonist-bound receptor concentration over time, shown in Scheme 2:

Scheme 2
RA is the receptor species that can couple to the response precursor to generate the response, referred to as the active receptor. R0A is the desensitized receptor, which cannot couple to the response precursor and so is inactive with respect to signaling. kDES is the receptor desensitization rate constant. In this particular model, the response does not decay and the response precursor concentration remains constant over time, i.e. it is not depleted by generation of the response.
The differential equation for the response is, The goal is an equation in which E is the only time-dependent variable. This can be achieved using The E vs t equation is obtained by taking the inverse Laplace transform:

A.3. Response degradation (Model 3)
In this model the response is degraded over time. There is no receptor desensitization or precursor depletion. This is Model 2 of the original kinetic model . The model is represented by Scheme 3:

Scheme 3
55 The E vs t equation for this mechanism has been derived previously  and is given here as Eq. A.5: At a saturating concentration of agonist, the equation reduces to Eq. 3: When agonist is not at equilibrium with the receptor, the equation, derived as Eq 10 in , is:

A.4. Response degradation and recycling (Model 4)
In a variant of the response degradation model (Model 3), the response degradation process results in the reformation of the response precursor. In other words, the response recycles to the 56 response precursor. This is Model 3 of the original kinetic model  and is represented by Scheme 4:

Scheme 4
Note in this variant the response precursor becomes partially depleted. The E vs t equation here, Eq. A.7 below, is a rearranged form of the equation in the original study (Eq. 3 of (Hoare et al.,

A.5. Depletion of response precursor (Model 5)
In this model, the only response-limiting process is depletion of response precursor. There is no receptor desensitization or response degradation. The response precursor is depleted by its conversion to the response during the response generation process. This model is represented by Scheme 5:

Scheme 5
Two species are time-dependent in this model, E and EP. The derivation reduces this to the single time-dependent variable of interest, E. This is done using a conservation of mass equation for E.
First the differential equation for E is written:

59
EP can be expressed as a function of the total amount of response material, EP(TOT), which remains constant over time and is defined by the following conservation of mass equation: Solving for EP and substituting into the differential equation for E gives, At saturating concentrations of agonist, the E vs t equation reduces to Eq. 5: For the slow agonist equilibration scenario, we were able to derive an analytical solution.
This involved handling a system of nonlinear differential equations.
In the specific case when 0 = 0, this simplifies to The [ ] ∞ term can be rearranged to: Substituting into the E vs t equation and rearranging gives the final equation, Eq. A.9: where,

A.6. Receptor desensitization and response degradation (Model 6)
In this model, the response is regulated by two processes; the receptor undergoes receptor desensitization, and the response is degraded. This model combines Models 2 and 3 and is represented by Scheme 6:

Scheme 6
The E vs t equation is obtained as follows. The differential equation for E is, The derivation proceeds by taking the Laplace transform for [RA] (Eq. B.4) and substituting it into the transform for E, giving, after some rearranging, The E vs t equation is now obtained by taking the inverse Laplace transform, giving Eq. A.10:

Eq. A.10
At saturating concentrations of agonist, the E vs t equation reduces to Eq. 6: For the slow agonist equilibration scenario, the Laplace transform for [RA] is Eq. B.6.
Substituting into the transform for E and rearranging gives, Taking the inverse transform gives the E vs t equation, Eq. A.11:

A.7. Receptor desensitization and response degradation with recycling (Model 7)
In this variant of the receptor desensitization and response degradation model (Appendix A.6), the response degradation process results in the reformation of the response precursor. In other words, the response recycles to the response precursor. The model is represented by Scheme 7:

Scheme 7
An analytical solution for E as a function of time could not be found for this system. The E vs t data can be simulated by numerical solution of the differential equations for E and EP. and the analytical solution for [RA]. The equations for E and EP are, The analytical equation for [RA] is Eq. B.3 for the rapid equilibration scenario and Eq. B.5 for the slow equilibration model.

A.8. Precursor depletion and response degradation (Model 8)
In this model, the precursor is depleted as it is being converted to the response and the response degrades. This is Model 4 of the original kinetic model . The model is represented by Scheme 8:

Scheme 8
The equation used in this study is a modified form of Eq. 4 in . A new, macroscopic term is introduced here called the depletion rate constant, kDEP. This term represents the rate at which the response precursor is depleted as it is being converted to the response. It is introduced as follows. The E vs t equation for the model is, as derived previously , Substituting [RA] for [ ] and introducing kτ gives, 67 kDEP is defined here as the product of the receptor concentration and the response generation rate constant:

Eq. A.12
At saturating concentrations of agonist, the E vs t equation reduces to Eq. 8: Eq. 8 For the slow agonist equilibration scenario, an analytical solution for E could not be obtained, as described previously . The E vs t data can be simulated by numerical solution of the differential equations for E and EP, and the analytical equation for [RA] (Eq. B.2). The differential equations are, We have discovered an analytical solution is obtainable for response precursor, in other words the decline of EP which results from generation of E. Owing to decoupling of [RA] from E (Appendix B.1), the analytical equation for [RA] (Eq. B.2) can be substituted into the differential equation for EP: which can be rearranged to, This ODE is related to the Gompertz growth/decay model (Kirkwood, 2015): This problem has a known solution: Applying to the initial value problem for EP gives the analytic solution for decline of the response precursor, Eq. A.13: where,

A.9. Receptor desensitization and precursor depletion (Model 9)
In this model, the receptor desensitizes and the precursor is depleted as it is being converted to the response. The model is represented by Scheme 9:

Scheme 9
70 The E vst t equation is obtained by first deriving the equation for EP and then subtracting this from EP(TOT) using the conservation of mass equation. The differential equation for EP is, Owing to decoupling of the ODEs for [RA] and E (Appendix B.1), [RA]t can be represented by a time-dependent coefficient in the differential equations and EP. This is Eq. B.3. Substituting into the differential equation for EP gives, This expression is an initial value problem of the Gompertz growth/decay model type (Kirkwood, 2015): This problem has a known solution: Applying to the expression for EP gives the analytic solution for decline of the response precursor (Eq. A.14): The E vs t equation, Eq. A.15, can now be determined from conservation of response precursor ( ( ) = + ):

Eq. A.15
For the case of a saturating concentration of agonist, is unity, and the E vs t equation reduces to Eq. 9: where, For the slow agonist equilibration scenario, we first derive an analytical equation for EP.
The governing differential equations are, The system decouples (Appendix B.1.). The resulting system for [R] and [RA] is inhomogeneous, with a forcing term: (1) and (2) are defined in Appendix B.2. We can now solve for EP by writing the differential equation for EP as a linear ODE: (2) − (1) Solving with , =0 = ( ) , we find the EP vs t equation, Eq. A.16: where, which can be rearranged to, The E vs t equation, Eq. A.17, is now obtained from the conservation of response precursor:

A.10. Receptor desensitization, precursor depletion and response degradation (Model 10)
In this model, all three regulation mechanisms are in operation. The receptor desensitizes, the response precursor depletes as it is being converted to the response, and the response degrades.
The model is represented by Scheme 10:

Scheme 10
For the rapid agonist equilibration scenario, an analytical solution for E could not be obtained. . The E vs t data can be simulated by numerical solution of the differential equations for E and EP, and the analytical solution for [RA]. The differential equations are: The analytical equation for [RA] is Eq. B.3 for the rapid equilibration scenario and Eq. B.5 for the slow equilibration model.
We were able to obtain a solution for response precursor, in other words the decline of EP which results from generation of E. These equations are the same as those for Model 9 (Appendix A.9), the receptor desensitization and precursor depletion model. (The additional step in Model 10 here compared with Model 9 is the degradation of E, which does not affect EP because there is no recycling.) For the rapid ligand equilibration scenario, the equation for EP is Eq. A.14 (and for the slow scenario, Eq. A.16 (Appendix A.9).

Appendix B: Receptor-agonist binding equations
Solving for [ ] ������ using the procedures described in (Hoare, 2017)  (2) − (1) The Laplace transform is used in the derivation for Model 2 (Appendix B. Data are from Fig. 3 of . In that study, a horizontal exponential equation was used to determine the observed rate constant (k) and the plateau (Bmax). Here, kτ is calculated by multiplying k by Bmax (Section 3.2.4). Data were extracted using a plot digitizer [WebPlotDigitizer (Rohatgi, 2018)]. and response degradation (pink). The response generation process is in green. In this process, the response precursor (EP) is converted to the response (E) by the agonist-occupied receptor (RA).

Figure legends
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 R0A that does not generate a response (because it can't couple to EP). The rate constant for desensitization is kDES. Degradation of response is represented by decay of E to D, governed by the degradation rate constant kD. This model is an extension the original kinetic operational model , extended to incorporate receptor desensitization.   (Kohout et al., 2001), extracted using a plot digitizer [WebPlotDigitizer (Rohatgi, 2018)].
Diacylglycerol was quantified using a live-cell biosensor. Diacylglycerol is degraded by diacylglycerol kinases and the AT1 receptor is susceptible to desensitization. This results in a riseand-fall exponential curve (Section 3.3.1.  Fig. 2A of (Violin et al., 2006), extracted using a plot digitizer [WebPlotDigitizer (Rohatgi, 2018)].   (Merelli et al., 1992), extracted using a plot digitizer (Graph Grabber 2.0, Quintessa Ltd., Henley-on-Thames, United Kingdom). . The E vs t profile is a Gompertz curve, which ascends more slowly than a horizontal exponential curve at later time points (Kirkwood, 2015). 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 agonistoccupied 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 operational model , extended to incorporate receptor desensitization.    Receptor desensitization: Angiotensin-stimulated inositol phosphates accumulation in mouse embryonic fibroblasts via the AT 1 angiotensin receptor (Kohout et al., 2001). Mouse fibroblasts were isolated from arrestin double knock out mice or wild type mice. In the absence of arrestin, the response is unregulated (straight line). In the presence of arrestin, the receptor becomes desensitized on application of the agonist, resulting in the horizontal exponential curve, consistent with the kinetic desensitization model (Section 3.2.1.). Data for wild-type cells were fit to the receptor desensitization model (horizontal exponential equation, Section 3.2.1), giving a k τ value of 0.10-fold.min -1 , and a k DES value of 0.042 min -1 (translating to a half-time for receptor desensitization of 17 min). Data for the arrestin knock out cells were fit to the no-regulation model (straight-line equation, Section 3.1), giving a k τ value of 0.13-fold over basal.min -1 . A maximally-stimulating concentration of angiotensin II was applied (100 nM). Data are from Fig. 3C of (Kohout et al., 2001), extracted using a plot digitizer [WebPlotDigitizer (Rohatgi, 2018)]. Receptor desensitization: GnRH-stimulated inositol phosphates accumulation via the rat GnRH 1 receptor with a C-terminal extension in HEK293 cells (Heding et al., 1998). The mammalian GnRH 1 receptor lacks a C-terminal tail and so does not undergo arrestinmediated receptor desensitization. Incorporating a C-terminal tail from a receptor that does (TRH receptor) results in arrestin-mediated receptor desensitization (Heding et al., 2000;Heding et al., 1998). Data for the C-terminally-extended receptor were fit to the receptor desensitization model (horizontal exponential equation, Section 3.2.1), giving a k τ value of 1.5%.min -1 , and a k DES value of 0.86 min -1 (translating to a half-time for receptor desensitization of 0.81 min). Data for the wild-type receptor were fit to the no-regulation model (straight-line equation, Section 3.1), giving a k τ value of 2.2%.min -1 . A maximallystimulating concentration of GnRH was applied (1 M). Data are from Fig 4 of (Heding et al., 1998), with basal response subtracted, extracted using a plot digitizer [WebPlotDigitizer (Rohatgi, 2018)].    (Davis et al., 1986). In the absence of Li + , inositol phosphates are degraded, resulting in the horizontal exponential curve, consistent with the kinetic response degradation model (Section 3.2.2.). Data for the absence of Li + were fit to the response degradation model (horizontal exponential equation, Section 3.2.2), giving a k τ value of 250 cpm IP 3 .min -1 . Data for the presence of Li + were fit to the no-regulation model (straight-line equation, Section 3.1), giving a k τ value of 580 cpm IP 3 .min -1 . A maximally-stimulating concentration of GnRH was applied (85 nM). Data are from Fig. 3 of (Davis et al., 1986), extracted using a plot digitizer [WebPlotDigitizer (Rohatgi, 2018)].  (Streb et al., 1985). In these cells, the CCK 1 receptor becomes desensitized on application of agonist (Klueppelberg et al., 1991) and IP 3 is degraded. This results in a rise-and-fall exponential curve (Section 3.3.1). Data were fit to Model 6, the receptor desensitization and response degradation model (see Section 3.3.1 for details), giving the following fitted values: k τ , 130 cpm [ 3 H]IP 3 .min -1 ; k DES , 0.14 min -1 ; k D , 0.39 min -1 . In the presence of Li + , the E vs t profile becomes a horizontal exponential curve, consistent with suppressed response degradation and with the remaining regulation process being receptor desensitization (Section 3.2.1). Data for the presence of Li + (i.e. desensitization alone) were fit to Model 2 (Section 3.2.1), giving fitted values of: k τ , 75 cpm [ 3 H]IP 3 .min -1 ; k DES , 0.20 min -1 . Data are from Fig. 10 of (Streb et al., 1985), extracted using a plot digitizer [WebPlotDigitizer (Rohatgi, 2018)]. Receptor desensitization with response degradation: Angiotensin II-stimulated diacylglycerol generation via the AT 1 angiotensin receptor in HEK-293 cells (Violin et al., 2006). Diacylglycerol was quantified using a live-cell biosensor. Diacylglycerol is degraded by diacylglycerol kinases and the AT 1 receptor is susceptible to desensitization. This results in a rise-and-fall exponential curve (Section 3.3.1.). Data for the wild-type receptor were fit to Model 6, the receptor desensitization and response degradation model (see Section 3.3.1 for details), giving the following fitted values: k τ , 2.7 DAGR ratio units.min -1 ; k DES , 0.74 min -1 ; k D , 4.1 min -1 . Receptor desensitization can be suppressed using a modified AT 1 receptor which lacks phosphorylation sites ("Desensitization deficient," open circles). The E vs t profile becomes a horizontal exponential curve, consistent with suppressed receptor desensitization and with the remaining regulation process being response degradation (Section 3.2.2). Data for the desensitization-deficient receptor (i.e. degradation alone) were fit to Model 3 (Section 3.2.2), giving fitted values of: k τ , 4.0 DAGR ratio units.min -1 ; k DES , 6.1 min -1 . Data are from Fig. 2A of (Violin et al., 2006), extracted using a plot digitizer [WebPlotDigitizer (Rohatgi, 2018)].  . The E vs t profile is a Gompertz curve, which ascends more slowly than a horizontal exponential curve at later time points (Kirkwood, 2015). A horizontal exponential curve is shown for comparison (dashed line). Data were simulated using Eq. 9 with the following parameter values: k τ , 5.0 response units.min -1 ; E P(TOT) , 10 response units; k DES , 0.20 min -1 . (c) Receptor desensitization, precursor depletion and response degradation (Model 10). Data were simulated numerically using the differential equations in Appendix A.10, with  A set to unity (representing a maximally-stimulating concentration of agonist) with the following parameter values: k τ , 20 response units.min -1 ; E P(TOT) , 40 response units; k DES , 0.20 min -1 ; k D , 0.030 min -1 . The dashed line is a fit of the simulated data to the generic rise-and-fall exponential equation (Section 3.3.4.), giving a k τ value of 17 response units.min -1 , K 1 of 0.48 min -1 and K 2 of 0.029 min -1 .

Effect of multiple agonist concentration, equilibrium agonist binding assumption
In Figure S1, response to multiple agonist concentrations is presented for each of the models, assuming agonist rapidly equilibrates with the receptor.

Model 2: Receptor desensitization
The t1/2 of the horizontal exponential curve is dependent on

[A] (decreasing as [A] increases)
whereas the plateau is the same for all [A] (Fig. S1B). As a result, EC50 decreases over time.

Model 3: Response degradation
The t1/2 of the horizontal exponential curve is the same for all [A] whereas the plateau is dependent on [A] (Fig. S1C). Consequently, EC50 is constant over time. See Fig. 10 of  for an example.

Model 4: Response degradation and recycling
Both the t1/2 and plateau are dependent on agonist concentration (Fig. S1D). EC50 decreases over time. See Table 2 of  for an example.

Model 5: Precursor depletion
The t1/2 of the horizontal exponential curve is dependent on [A] (decreasing as [A] increases) whereas the plateau is the same for all [A] (Fig. S1E). As a result, EC50 decreases over time.

Model 6: Receptor desensitization and response degradation
In this rise-and-fall exponential curve, the gradient of the rise, and the peak value, are dependent on [A] (Fig. S1F). Response declines to zero.

Model 7: Receptor desensitization and response degradation with recycling
A distorted rise-and-fall curve is evident in which the decline phase is delayed, resulting in a bulge on the right-hand, downward part of the curve. The gradient of the rise, and the peak value, are dependent on [A] (Fig. S1G). Response declines to zero.

Model 8: Precursor depletion and response degradation
In this rise-and-fall exponential curve, the gradient of the rise, and the peak value, are dependent on [A] (Fig. S1H). Response declines to zero. Note at later time points the complex relationship between response and [A]. Numerous examples are provided by cytoplasmic Ca 2+ mobilization response data: Fig. 1 of (Princen et al., 2003), Fig. 4A of (Malysz et al., 2009), Fig. 5 of (Kassack et al., 2002), and Fig. 1 of (Milligan and Rees, 1999).

Model 9: Receptor desensitization and precursor depletion
The response is described by a Gompertz curve, which resembles the horizontal exponential curve but it ascends more slowly at the later time points (see Fig. 9B). The initial gradient is dependent on [A] whereas the plateau is the same for all [A] (Fig. S1I).

Model 10: Receptor desensitization, precursor depletion and response degradation
The response closely approaches a rise-and-fall exponential curve (see Section 3.4.3 for details).
The gradient of the rise, and the peak value, are dependent on [A] (Fig. S1J). Response declines to zero.

Effect of multiple agonist concentration, slow agonist equilibration assumption
For these models, the curve shapes resemble the corresponding shape for the rapid equilibration but there is a clear lag in the rise phase of the curves (compare the corresponding panels of Figs. S1 and S2). The lag is dependent on the agonist concentration. The higher the concentration, the less evident the lag, a result of the increasing rate of agonist-receptor association due to mass action.

Depletion of response precursor models
In some of the more complicated models, it was possible to obtain an analytical solution for EP even though it was not possible to obtain the solution for E. In some experimental paradigms it is possible to measure the decline of the precursor, for example the decline of phosphatidylinositol 4,5-bisphosphate (Tewson et al., 2013). Data for these models are simulated in