To produce a response, a ligand must first bind to the receptor. Binding parametrization is typically achieved via the widely used equilibrium dissociation constant (K_d ), which is measured in units of concentration (usually molarity, M). In the simplest one-state model, where all occupied receptors are active and that forms the basis of the widely used one-parameter Clark equation ( Figure 1F ), K _d is expressed as a function of the equilibrium concentrations of the species involved as:

It is important to remember that K is not just an arbitrary model parameter: it is an entity measurable in equilibrium or kinetic binding assays, and it is related to the Gibbs free energy of binding (ΔG) via the well-known thermodynamic equation:

Here, T is the absolute temperature and R the universal gas constant, R = k _B N _A = 8.314 J/K·mol. Accordingly, to bind with 1 nM affinity (K _d = 10^–9 M) at physiological temperature (T = 310 K), a ligand requires a free energy of binding of ΔG ⁰ = –RT·lnK _d = 53.4 kJ/mol (12.8 kcal/mol).

The present model uses a similar definition for K _d. It differentiates between active and inactive states (denoted by an asterisk; R vs. R* and LR vs LR*), but contrary to most two-state models that distinguish between binding affinities for the active and inactive states (e.g., K _d and K _d/α; Figure 2 ), it considers ligand-bound and ligand-free states as ensembles of conformations characterizable by the same K _d. Hence, K _d here represents an average binding constant for the ensemble of active and inactive receptor forms that the ligand effectively sees. With this assumption, the definition of K _d will rely on the total concentration of occupied and unoccupied receptors, which in the case of the full two-state model ( Figure 1A ) becomes

This is a macroscopic equilibrium constant measurable in equilibrium binding assays that assess total binding to the receptor. In case of the minimal two-state model, where there is no constitutive activity ( Figure 1C , [R*] = 0), it reduces to the previously used form (Buchwald, 2017):

Binding, characterized by K _d, is not just snapping of a ligand into a rigid site (“key-in-the-lock”), which corresponds to that of the active conformation (inactive for an antagonist), but it also involves change in energy and possibly adaptation of the site to accommodate the ligand (“hand-in-the-glove”). The corresponding conformational shift (toward active states) leads to change in the downstream activation (induced fit) ( Figure 3 ). A key advantage of this definition is that K _d represents ensemble averages for the bound and unbound receptor forms and becomes independent from the value of the efficacy parameter (Buchwald, 2017). Another one is that it avoids the loop-related thermodynamic constrains (i.e., need for a same α for both K _d and K_ε in Figure 2 ) as will be discussed in the section Mechanism of Receptor Binding and Activation.

The need for efficacy parametrization, besides that for binding, arose from the recognition of partial agonism, i.e., the existence of compounds that can occupy all receptors without achieving full activation. Well-known partial agonist examples include prenalterol (versus the full agonist adrenaline at β-adrenoceptors), pilocarpine (versus acetylcholine at muscarinic receptors), and impromidine (versus histamine at H₂ receptors) (see Supplementary Figure S1 for representative chemical structures). Here, efficacy parametrization is done with an ε parameter that represents the fraction of ligand-bound receptors that are active (Buchwald, 2017):

Hence, efficacy ε as defined here is a unitless parameter (but not an equilibrium constant such as K_ε ) and is an intrinsic efficacy in the sense that it is a characteristic of the ligand L (for a given receptor R). It ranges from 0 for an antagonist that can occupy all receptors but produces no effect to unity (1) for a full agonist that converts all occupied receptors to active ones. The activity here is the one measured immediately post-receptor; downstream signal amplification can complicate things as response from a partial agonist can be amplified to a maximum final response (see discussion later). A set of representative response curves obtained for a fixed binding affinity (K _d) but different efficacies ε are shown in Figure 4A .

Note that the efficacy ε as defined here is different from the efficacy e as defined by Stephenson, which was introduced to measure the ability of a drug to produce a response in a tissue. In Stephenson’s definition, the stimulus was S = e[LR*]/[R_tot], and e could have values from 0 up to infinity (Stephenson, 1956). In fact, ε is more similar to the efficacy as defined by Furchgott (ϵ _F = e/[R_tot]) representing “the capacity of a ligand to initiate stimulus from one receptor” (here, as measured right after the receptor before additional signaling steps). It is also similar to the intrinsic activity α _A introduced by Ariëns (1954), which is the ratio of the maximum response produced by the partial agonist to that produced by the full agonist and is ultimately quantifying “effect per unit of pharmacon–receptor complex” (as long as there is no amplification).

With these assumptions, the fractional response (effect) f _resp, which under the assumptions of no amplification ( Figure 1D ) equals the fraction of activated receptors f _act, can be expressed as a function of ligand concentration only as (Buchwald, 2017):

This is essentially equivalent with a Clark equation that includes a scalable maximum (E _max model for partial agonism; Figure 1D ). It is also equivalent with the model proposed by Ariëns (e.g., Equation 12 in Ruffolo, 1982) incorporating an intrinsic activity term (α _A). Equation 10 links the definition of the maximum achievable effect (response) for a partial agonist, f _resp,max = ε, to a two-state model ( Figure 1D ). With this definition, ε corresponds to the fraction of maximum activation a partial agonist can achieve as compared to the full agonist—as measured right after the receptor (to avoid possible confounding effects created by downstream amplification; see later). For full agonists, ε = 1 and the above equation correspond directly to the Clark equation for response (Clark, 1926; Clark, 1933) ( Figure 1F ), which is mathematically equivalent with the Hill–Langmuir equation for ligand binding (Hill, 1909) and a special case of the versatile Hill equation (Hill, 1910) often used in pharmacological and other applications (Goutelle et al., 2008; Gesztelyi et al., 2012).

Since the introduction of the concept in the late 1980s (Costa and Herz, 1989), it is now well recognized that certain G-protein–coupled receptors (GPCRs) can be active even in the absence of an agonist (have constitutive signaling activity) and that some ligands can act as inverse agonists (i.e., reduce the activity of the ligand-free receptor) (Bond and Ijzerman, 2006). To incorporate such constitutively active receptors into the formalism of the present model ( Figure 1A ), a baseline receptor efficacy (ε _R0) is introduced along lines similar to the definition of ε (Equation 9), i.e., the fraction of unbound receptors that are active:

However, contrary to ε, which is a ligand characteristic, ε _R0 is a receptor and not a ligand characteristic. For such receptors, ligands known as inverse agonists reduce the signaling output below that of the basal state. In the present formalism, they will have efficacies that are smaller than this baseline receptor efficacy, ε _inv.agon.<ε _R0. Partial agonists remain those compounds that generate a response, albeit a submaximal one even at concentrations that saturate all receptor sites (ε _R0<ε<1), and so-called neutral antagonists do not discriminate between the two conformers (ε _R0=ε) and thus can block agonist, partial agonist, and inverse agonist activities. Some representative response curves that can be obtained with this parametrization are shown in Figure 4B , and further discussion is included later.

Some kind of signal amplification is needed to account for the existence of “spare receptors” or “receptor reserves”: cases where almost maximal response can be achieved by occupying only a (small) fraction of all receptors. The concept has been introduced in the mid-1950s (Nickerson, 1956; Stephenson, 1956) with a method for its quantification described by Furchgott about a decade later (Furchgott, 1966; Furchgott and Bursztyn, 1967). Well-known extreme cases include, for example, the response of the human calcitonin receptor type 2 to calcitonin (where 20% occupancy already produces close to 100% response) (Chen et al., 1997) and guinea pig ileal response to histamine (where 2% occupancy already produces close to 100% response) (Kenakin and Cook, 1976; Adham et al., 1993; Kenakin, 2018b). Another example is the stimulation of β-adrenergic receptors in the heart by epinephrine where half-maximal increase of muscle contractility already occurs at 1% to 3% receptor occupancy in rats and at 10% to 20% occupancy in humans, and the effects saturate well before the receptors do (Brown et al., 1992). Detailed quantifications, including for different ligand series, have been done for several cases, typically GPCRs, such as the muscarinic acetylcholine (Furchgott and Bursztyn, 1967; Harden et al., 1986; Eglen and Whiting, 1987), opioid (Chavkin and Goldstein, 1984; Adams et al., 1990; Fox and Hentges, 2017), dopamine (Meller et al., 1987), 5-hydroxytryptamine (5-HT) (Meller et al., 1990), A₁-adenosine (Dennis et al., 1992; Morey et al., 1998), and cannabinoid (Gifford et al., 1999) receptor systems.

Most commonly, such cases are due to strong signal amplification via a downstream cascade of saturable functions provided by second messengers or other systems (Strickland and Loeb, 1981; Koshland et al., 1982; Ferrell, 1996). Activation of adenylate cyclase and production of cyclic AMP (cAMP) as second messenger followed by further downstream cascading is a well-known example. For example, the mechanism by which epinephrine (acting on a myocyte) or glucagon (acting on a hepatocyte) sets off a cascade of phosphorylation leading to the production of glucose involves downstream amplification of the initial signal by several, possibly up to eight orders of magnitude (Nelson and Cox, 2012). As a result, epinephrine concentrations as low as 10^–10 M in the blood can stimulate liver glycogenolysis and release of physiologically adequate levels of glucose. Such a small epinephrine stimulus can generate intracellular cAMP concentrations of 10^–6 M (already a ∼10,000-fold gain) followed by three more catalytic steps leading to release of glucose via another approximately 10,000-fold amplification (Lodish et al., 2003). In many ways, the resulting cascade mechanisms of signal amplification in biological systems resemble those used in electronic circuits (Grubelnik et al., 2009). Very strong amplification is not always desirable as amplification of noise-level signals has to be avoided (Bialek and Setayeshgar, 2005).

In the present model, pathway-specific signal amplification is incorporated explicitly using a separate gain parameter γ. This is an important novel component that has not been explicitly included as such in previous models like the operational (Black and Leff) or the minimal two-state (del Castillo–Katz) model. Amplification is built in via use of a nonlinear hyperbolic-type response function (Buchwald, 2017). Hyperbolic relationships between fractional occupancy and fractional response have been found in several cases (see examples later), justifying the use of such functions. Use of hyperbolic functions also has the advantage that even if there is a cascade of such sequential saturable signal amplification functions with the output of one serving as input for the next one (as is often the case in biological systems), the overall response can be represented by a single hyperbolic-type response function (Trzeciakowski, 1999a). The nonlinear response (amplification) function used here is a classic hyperbolic function, just as it is in the operational model. However, to be able to reach its asymptotic limit, not [LR*] is used as its input, because [LR*] cannot be higher than [R_tot] and, hence, cannot cover the asymptotic part of the response function that is beyond the [R_tot] value, but its odds-ratio type (De Muth, 2014) transform: Λ = p/(1–p) (see Buchwald, 2017 for further details). In other words, the input for the hyperbolic function here is not [LR*], but Λ = [LR*]/([R_tot] − [LR*]) = f _act/(1 − f _act), with f _act = [LR*]/[R_tot] representing the fraction of active receptors. For full agonists, f _act corresponds to the fraction occupied f _occup; f _act = f _occup if ε = 1. This extends the range of the input from 0–1, which is the range for f _act (i.e., [LR*]/[R_tot]), to 0–∞. ¹ Hence, the response function linking E to [LR*] will be:

After some transformations and introduction of γ=[R_tot]/K_γ , this leads to the final general form for the present three-parameter model ( Figure 1C ) (Buchwald, 2017):

Here, the just introduced γ parameter represent a unitless amplification (gain) factor. Since it is a gain, for all practical purposes, it has a value larger than unity, γ≥1. A set of illustrative response curves for a fixed value of K _d and different values of ε and γ are shown in Figure 4A . For a given ligand acting at a specific receptor, affinity (K _d) and efficacy (ε) should be the same (as long as receptor-binding and activation are not influenced by the environment, and signaling bias, which will be discussed later, is not considered), but transduction (signal amplification) could be pathway and tissue-dependent. Hence, γ, and, in fact, both of its components, R_tot and K_γ , can be tissue-dependent.

A slightly rearranged form of this equation provides a better understanding of the interplay between its parameters:

From here, it is clear that half-maximal activity (EC₅₀) is observed at K _obs = K _d/(εγ+1–ε), and maximum (fractional) effect achievable by a given ligand is f _resp,max = εγ/(εγ+1–ε). Hence, for a full agonist at the receptor (ε = 1), K _obs = K _d/γ and γ is a straightforward multiplication factor causing a left-shift of the sigmoid response function by γ units on a semi-log scale. Thus, for such an agonist, signal amplification causes no change in the shape of the response on semi-log scale, just a left-shift by γ increasing the apparent potency γ-fold (as illustrated by the blue lines in Figure 4A ). For an agonist that produces only partial activation at the receptor (ε<1), downstream amplification can increase the maximum response, and sufficiently strong amplification can transform such a partial agonist into an apparent full or close to full agonist. With a large enough γ, the maximum response, f _resp,max = εγ/(εγ + 1 − ε), can approach unity even for small εs. However, for such partial agonists, the left shifts, (εγ + 1 − ε), are smaller than for a full agonist, (1·γ + 1 − 1 = γ), so that the change in the apparent EC₅₀ (K _obs) is less (red vs. blue lines in Figure 4A ).

Finally, to include the effect of amplification for cases with constitutive activity (completely general case represented by model A in Figure 1 ), we will calculate the fraction of activated receptors f _act, which is proportional with the effect right after the receptor, and use it as input for the amplification function. With the definitions of ε _R0, ε, and K _d, the receptor concentrations can be eliminated, and f _act can be expressed as a function of [L] as:

As before, this will serve as input via Λ=f _act/(1–f _act) for the present amplification function (Equation 12), resulting in the final form of the full four-parameter SABRE model:

It can be seen from here that basal response ([L] = 0) is ε _R0 γ/(ε _R0 γ + 1 − ε _R0), whereas the maximum effect for a ligand L ([L]→∞) remains f _resp,max = εγ/(εγ+1–ε) (see Supplementary Figure S2 for an illustration of the effects of different ε _R0 and γ on the response calculated with this equation). Rearranging this in a manner like that done for Equation 14 but also separating the basal response leads to

Hence, the half-point of the transition is at

Consistent with the true generalized nature of the present model, if there is no constitutive activity, ε _R0 can be set to 0, and this recovers the simplified forms of the model that have been introduced earlier (Buchwald, 2017) (models C–F, Figure 1 ). Conversely, if there is no amplification ( Figure 1B ), γ = 1 and the total effect becomes:

As illustrated in Figure 4B for some representative parameter values, this can describe concentration-response functions for full (ε = 1), partial (ε>ε _R0), and inverse agonists (ε<ε _R0), as well as neutral antagonists (ε=ε _R0).

The combination of all these effects can result in quite complex concentration-response curves that may make fitting with well-defined parameters difficult (see below). Nevertheless, the present four-parameter model can account for cases where the mixture of signal amplification and partial agonism causes complex responses, and several specific examples will be discussed (see also Buchwald, 2017 for the case of competitive partial agonism). It can also account for the observation of different responses created by the same agonists if observed along different pathways (biased agonism; e.g., Figures 13 and 14 ) or at different vantage points on the same downstream pathway involving multiple amplification steps (e.g., Figure 10 ).

Before closing this section on parametrization, a few general modeling related considerations have to be highlighted. The present model can fit complex data such as those illustrated in the section Fit of Complex Fractional Response Versus Fractional Occupancy Data; however, use of its fully parameterized version (Equation 2 or even Equation 4) only makes sense if sufficient data are available, and if occupancy (binding) and response (effect) can be assessed independently and under sufficiently similar conditions. Since adequate model fitting requires the availability of 5–10 (well distributed) data points for each adjustable parameter (Knofczynski and Mundfrom, 2008; Austin and Steyerberg, 2015; Kenakin, 2018b), reliable fitting of the full model can only be accomplished if a sufficiently large number of data points are available. Along these lines, it has to be mentioned that the present model uses one more parameter than the operational model, e.g., three (K _d, ε, γ) vs. two (K _d, τ) for the case of no constitutive activity or four (K _d, ε, γ, ε _R0) vs. three [e.g., K _d, ε, χ (Slack and Hall, 2012)] for cases with constitutive activity. This is relevant because obtaining well-defined parameter values could require more data points, and rigorous model selection criteria advocate the use of the simplest model that can still provide adequate fit (George, 2000; Myung and Pitt, 2004; Buchwald, 2005; Buchwald, 2007). However, the need for one extra parameter is more than compensated for by, on one hand, the intuitive nature of the present parameters (due to separation of efficacy in receptor activation from gain in signal amplification), and, on the other, the ability to use simplified forms with reduced number of parameters. Contrary to the operational models, with the present one, simplified forms can be recovered for special cases of its parameters, and these can and should be used on their own when adequate or when there is not enough data to support full parametrization. Note that Hill type extensions that involve an additional n _H parameter to account for more or less abrupt response functions are not discussed here.

Transforming the classic form of the operational model (Equation 20) in a manner similar to that done for the present model (Equation 14 vs. 13) leads to

Hence, for the operational model, the maximum fractional response achievable by a ligand is f _resp,max = τ/(τ+1) and half-maximal activity occurs at EC₅₀ = K _obs = K _d/(1+τ). Because for full (or close to full) agonists, the maximum (fractional) response, τ/(τ+1), needs to be close to 1, τ needs to have large values, and those are difficult to obtain in a well-defined manner, as the τ/(τ + 1) fraction is no longer sensitive to changes in τ when approaching unity. Hence, fitting by nonlinear regression can result in large and ill-defined τ values, making the linked K _d values also badly defined. Since the calculated K _d is (τ + 1)-fold different from the observed K _obs, the operational model (Equation 20) can end up not just with ill-defined, but also unrealistic K _d values that are essentially meaningless from a binding perspective. Hence, this model is difficult to fit for full or close to full agonists, and results can be cumbersome to interpret (see, e.g., Table 1 and Figure A1 in Buchwald, 2017 for specific illustrations). In agreement with this, mathematical identifiability analysis and simulation for the operational model has shown that when only functional assay data are available, only the transduction coefficient (τ/K _d) and not τ can be estimated precisely (Zhu et al., 2018). There are indeed applications employing the operational model either only to determine unresolved τ/K _d ratios (Kenakin et al., 2012) or with experimental K _d values to constrain the regression (Rajagopal et al., 2011). Because in the present model, ε is restricted to the 0 to 1 range and K _d is independent of ε, fitting by nonlinear regression does not face these issues, and all parameters can be fitted in well-defined manner for partial and full agonists as long as there are sufficient data points. Several specific examples are included in the next chapter for illustration.

As highlighted by Equation 21, the K _d parameter of the operational model obtained from data fitting is different from the apparent (observed) K _obs as its value also depends on τ: K _obs = K _d/(τ+1). Although this allows the concentration–response curve to shift from the concentration–binding curve, which is why such complex models are needed in the first place, it also makes K _d an empirical parameter not directly related to binding. As discussed above, for full or close to full agonists, τ needs to be large [e.g., τ>10 is needed for E _max = τ/(τ+1)>0.9], and K _obs values will be shifted considerably compared to K _d—often to unrealistically high values that are clearly far from the actual binding affinity of the ligand. Hence, it has become accepted to use K _d as an empirical parameter not necessarily related to receptor binding, and in some implementations, such as in GraphPad Prism, full agonists are not fitted at all with this model. It has been pointed out that for the operational model, changes in binding (K _d) and in conformation (τ) become indistinguishable for very efficacious agonists making interpretations difficult and cumbersome (Colquhoun, 1998).

Contrary to this, all parameters of the present model are straightforward, intuitive, and clearly related to their corresponding processes. Affinity (binding) is characterized by K _d, the binding constant, and values measured in equilibrium assays can be used directly in the model. This K _d is uncoupled from the post-binding ability to activate the receptor (efficacy) as well as the strength of the post-activation amplification. Efficacy (ability to activate the receptor) is characterized by ε, a ligand-specific unitless parameter ranging from 0 (for a ligand that keeps all receptors inactive) to 1 (for an agonist that shifts all receptors into active state). Nonlinear transduction, due to post-receptor amplification (gain), is characterized by γ, a unitless parameter ranging from 1 (no amplification) to infinity.

The present model also overcomes an essential hypothesis-related problem of the operational (Black and Leff) and del Castillo–Katz models, namely, that although the final equations are mathematically identical (Equation 20), they arrive at it from two conceptually different approaches that are both incomplete. On one hand, the minimal two-state (del Castillo–Katz) model ( Figure 2 without an R* state), which is now a generally accepted approach to describe switching of the receptor between active and inactive states, does not formally incorporate signal amplification (nonlinear transduction), which, as discussed, is known to exist. On the other, the operational (Black and Leff) model allows nonlinear response, but it is a single-state model that does not formally incorporate the possible existence of active and inactive ligand-bound receptor states (i.e., the possibility that not all occupied receptors are active), which is also known to exist. It makes up for this by limiting its output function to a maximum of τ/(τ+1). Consequently, these models, in fact, merge together two different effects in their τ parameter: the “intrinsic efficacy” of the (bound) ligand to activate the receptor, which can lead to partial activation (partial agonism) even with a linear response function, and the “efficacy” of the post-activation amplification downstream from the receptor, which can create fractional response in excess of the fractional occupancy (“receptor reserve”) and can be tissue- or organ-specific. These issues are overcome by the present model, by the introduction of separate efficacy ε and gain γ parameters that are also decoupled form the binding affinity, K _d.

Finally, another important advantage is that, contrary to previous models, the present one is a true generalized model: simplified forms can be recovered as special cases of its parameters, e.g., ε _R0 = 0 for no constitutive activity, γ = 1 for no amplification (no receptor reserve), and ε = 1 for full agonism only ( Figure 1 ). When adequate and/or when there is not enough data to support full parametrization, these simplified forms can and should be used on their own. The operational model, and consequently all of its extensions, cannot be reduced back to simpler forms, such as the E _max model for partial agonism or the Clark equation, as there are no τ values for which Equation 20 converts back to any of them (see Buchwald, 2017 for details). Therefore, one cannot transition back to simpler forms despite being in general desirable for complex models to be able to recover simpler ones for special cases of their parameters. In contrast, the present general model can be reduced back to a whole series of consecutively nested simpler forms ( Figure 1 ). For example, ε _R0 = 0 indicating no constitutive activity, reduces the general four-parameter model ( Figure 1A ) to the three-parameter model used before ( Figure 1C ). Further, if there is no post-receptor amplification, γ = 1, this model ( Figure 1C ) collapses back to an E _max model for partial agonists with efficiency ε ( Figure 1D ):

Finally, if only full agonists are considered (ε = 1), all occupied receptors are active ( Figure 1F ), and the corresponding equation collapses back to the well-known Clark equation, which forms the basis of the entire quantitative receptor theory:

Consequently, one can start fitting with the simplest most adequate form corresponding to constrained parameters, and these constrains can be removed sequentially as needed based on the complexity of the data being fitted and the number of data points available. For the simplest case, one can start with ε _R0 = 0 (no constitutive activity), γ = 1 (no amplification—fractional response overlaps with fractional occupancy), and ε = 1 (full agonism only) corresponding to a Clark model, and then remove these constrains as needed. On the other hand, the general model with all of its parameters released can fit complex cases where the fractional occupancy and the fractional response data are measured separately and do not correspond. In many cases, such data cannot be fit within the formalism of the operational model (e.g., to accommodate measured K _d values); a number of specific examples will be discussed below for illustration.

With the present model, one can directly link the fractional response, f _resp = E _{/ E max}, to the affinity-determined fractional occupancy, f _occup, which here is

To do so (for the case of no constitutive activity), [L] in the expression of f _resp (Equation 13) is replaced with its expression as a function of f _occup, from Equation 26 above, [L]=f _occup K _d/(1–f _occup), leading to:

Because of its two parameters (ε, γ), the present model allows quite flexible profiles. To avoid overparametrization-related problems, fit for a given response-system should be done using a common gain (γ) parameter across all ligands, and this should allow a good estimate of the overall amplification. At full occupancy (f _occup = 1), the maximum fractional response is f _resp,max = εγ/(εγ − ε +1), which can be less than unity for low-efficacy compounds. ³

These data involve concentration-dependent responses in a series of imidazoline type α-adrenoceptor agonists (e.g., phenylephrine, oxymetazoline, naphazoline, clonidine, and tolazoline) (Ruffolo et al., 1979) and are often used as a textbook example to illustrate mismatch between fractional receptor occupancy and response (Rang et al., 2015). Contractions of isolated rat aorta were measured as response, and receptor binding (K _d) was assessed separately by two different methods. Response alone can be fitted well by a standard E _max model (Equation 10; Supplementary Table S1.B) or the operational model (Equation 20) giving fits of essentially identical quality (Buchwald, 2017). However, fit with the operational model results in K values that are different from the measured K _ds and are quite meaningless for the full agonist phenylephrine (τ = 10^6.0; logK = −1.53 vs. logK _d = −6.46) (Buchwald, 2017). While these models can be used for empirical fitting of the response data, they cannot connect the response to the independently determined binding data (K _d). This can be done, however, with the present model (Equation 4) without significant loss in the quality of fit while also reducing the overall number of adjustable parameters (n = 6, 1 γ + 5 εs, versus the previous n = 10, 5 compounds, each with a K _d and E _max or K _A and τ) ( Figure 9A ; Table S1.C). Decreasing the number of parameters is important as model simplification is always an essential consideration (George, 2000; Myung and Pitt, 2004; Buchwald, 2005; Buchwald, 2007). Hence, the present model can do more than just an empirical fit of the response data and can account for complex cases where f _resp can either exceed or lag behind f _occup depending on the ligand ( Figure 9B ). Fit here predicts a reasonable 11.9-fold amplification for this system (γ = 11.88 ± 2.02) and intrinsic efficacies ranging from 1.0 for phenylephrine to 0.009 for tolazoline (Table S1).

Direct fit of the f _resp vs. f _occup data ( Figure 9B ) with the newly derived Equation 27 here constrained to a single amplification γ parameter for all compounds also results in very good fit and essentially identical ε and γ parameters (Table S1). The value of the correlation coefficient, r ² = 0.995, indicates that the model accounts for 99.5% of the variability in the quite complex response vs. occupancy data. Note that for oxymetazoline the fractional response exceeds the fractional occupancy at low occupancy (f _occup<50%), but lags behind it at higher ones.

A second set of data used for illustration involves the M₃ muscarinic receptor and a set of seven agonists including acetylcholine, carbachol, methacholine, oxotremorine, pilocarpine, and others (Sykes et al., 2009). This provides an even more complex test as two different responses were measured at consecutive vantage points after receptor activation: the stimulation of GTP binding to Gα subunits and the subsequent increase in intracellular calcium levels. Occupancy estimates are from pK values obtained from equilibrium competition experiments with N-methyl-[³H]scopolamine, but dynamic measurements were also performed to assess association and dissociation rate constants (k _on, k _off), and they resulted in somewhat shifted, but very similar pK values (r ² = 0.99) (Sykes et al., 2009). If all data are placed on the same graph, it is already apparent that the amplification for the two subsequent responses quantified here is quite different: while it has to be close to unity for GTP binding (as K _ds and K _obs≈EC₅₀s roughly overlap), it has to be around four orders of magnitude for Ca increase (as the corresponding EC₅₀s are left-shifted by about four log units) ( Figure 10A ). Fit of each of the two concentration-dependent responses with the present model accounts for 98% of the variability in the data ( Figure 10B ). However, because amplification is very strong in the Ca increase assay and all agonists achieve essentially maximum activation, well-defined efficacies cannot be determined from these data alone. Therefore, a unified fit was performed on both responses using only a single set of ε efficacies (one each for the seven compounds) plus two γ gain parameters (one for each response) as adjustable parameters (Table S2). This confirmed that the two amplifications are indeed very different as the γ values obtained from this overall fit were 2.1 and 10,089 for the GTP and Ca readout, respectively. The overall fit of these two response data sets with the same set of only 9 adjustable parameters was again quite good accounting for 95% of the variability in the data (r ² = 0.95; Figure 10A ) with betanachol, acetylcholine, and maybe pilocarpine fitting less well. Obtained relative overall efficacies are summarized in Table S2. In agreement with the original observation, which looked at intrinsic activities and log τ values obtained from the operational model (Sykes et al., 2009), efficacies (ε) obtained here also correlate very well with the pk _off values (r ² = 0.86) while having no correlation with the pK _d affinities (r ² = 0.01). In conclusion, the present model has the potential to connect activity data assessed at different vantage points and involving different amplification steps along the same downstream pathway to affinity and intrinsic efficacy data for complex series of partial agonists.

Within the framework of the present model, there is a specific connection between the stimulus input, represented by the concentration of active receptors [LR*], and response as shown by Equation 12. The loss of total receptors available for ligation due to irreversible inhibition, [R_tot]’=q[R_tot], leads to a corresponding loss in the [LR*] related stimulus. From the definition of the efficacy ε (Equation 9), the concentration of active receptors is

After inactivation leaves only a q fraction of R_tot, this becomes

Assuming that inactivation does not affect the post-receptor signal amplification function (Equation 12), which remains the same γ-dependent function, and only the input stimulus is altered, the fractional response after inactivation will be

After some transformations, this leads to:

By comparing this to the expression of response before inactivation (Equation 4), it can be seen that within the framework of the present model, a q-fold decrease in R_tot translates into an apparent q-fold reduction of efficacy, ε’=qε. Hence, the model can be used with this formalism of “fractional efficacy” ε’ to fit partial inactivation data obtained in Furchgott type experiments; two examples will be discussed below.

A notable problem with the receptor reserve concept is that it was originally defined as the fraction of receptors not required to achieve maximal response (for a full agonists) (Neubig et al., 2003). However, because the transduction functions linking response to occupancy are always of asymptotic nature (i.e., they asymptotically approach a limited maximum response as ligand concentration increases), in most cases, virtually all receptors are needed for maximum response, but often, only a relatively small fraction is needed for an almost maximal response. Hence, the “reserve” or, in other words, the discrepancy between the measured f _resp and f _occup varies strongly depending on where (i.e., at what response level) it is assessed. As the half-maximal point is often used as a reference point, a frequently used way to express how much “reserve” or “spare” receptors are is to estimate the percent of receptors occupied that already produces half-maximal response (e.g., (Furchgott and Bursztyn, 1967; Kenakin and Cook, 1976; Meller et al., 1987; Adham et al., 1993; Chen et al., 1997; Morey et al., 1998; Kenakin, 2018b)). For the present model, this can be obtained from the reverse of Equation 27 linking f _resp and f _occup:

Hence, in a system with a gain parameter γ, a ligand with efficacy ε will produce half-maximal effect (f _resp = 0.5) at

For a compound that produces full activation at the receptor (ε = 1), half-maximal effect is produced at f _occup = 1/(1 + γ). This means that in a system with 10-fold amplification (γ = 10), half-maximal response is produced at a receptor occupancy of 9.1%, while in a system with 100-fold amplification, it is produced already at 0.99% occupancy. For agonists with lower efficacies (ε < 1), these values will be different even in cases where the amplification is strong enough. For ligands with low enough efficacies, half-maximal response might not even be achieved at all, as γ + 1 > 1/ε is needed to have f _occup<1 in Equation 34. This interplay between efficacy and amplification can explain why the discrepancy between f _resp and f _occup can be different for different compounds even if tissue response (amplification) remains the same—a problem that confounded these studies and the notion of receptor reserve for a long time (Kenakin, 1986).

A first illustration is provided with a detailed data set where amplification could be estimated with two different approaches. Data are for a series of compounds acting at the dopamine receptor, which can be inhibited by the irreversible antagonist N-ethoxycarbonyl-2-ethoxy-1,2-dihydroquinoline (EEDQ) (Meller et al., 1987). In vivo dose response curves were generated in rats for the dopamine agonist reversal of γ-butyrolactone-induced striatal L-DOPA (L-3,4-dihydroxyphenylalanine) accumulation for four compounds: N-propylnorapomorphine (NPA), EMD 23,448, and (+/–)3-(3-hydroxyphenyl)-N-n-propylpiperidine [3-PPP(+) and 3-PPP(–)]. Data from one set of experiments following the Furchgott method indicates that response to NPA follows a typical hyperbolic response and is dose-dependently inhibited by EEDQ ( Figure 11A ). Fit with the present model (Equation 33, using a single K _d for NPA and a single γ for the response pathway) is very good (r ²=0.99) and suggests an about 43-fold signal amplification (γ=42.6±65.1) with EEDQ caused fractional inhibitions (q) estimated from the ε values that are in good agreement with those obtained in the original paper (Meller et al., 1987) using the classic Furchgott approach (Equation 29) (Table S3). The common K _d of the fit (19.0 μg/kg) also agrees well with that obtained via regression from the Furchgott approach for the highest EEDQ dose (24.4 μg/kg).

Next, the data obtained in this assay for all four compounds, but without EEDQ-mediated inhibition can also be fitted with the present model using the f _resp vs. f _occup method (Equation 27) as applied previously ( Figures 9B and 10B ). Fractional responses plotted as a function of f _occup (with occupancies calculated from K _ds values derived originally with the Furchgott method; Figure 7 in Meller et al., 1987) were fitted directly with the present model using Equation 27 with a single γ value. This also results in very good overall fit that accounts for as much as 95% of the overall variability in these data (r ²=0.95) ( Figure 11B ). The relative efficacies as compared to NPA as a full agonist obtained here are 0.19, 0.15, and 0.04 for EMD 23,448, 3-PPP(+), and 3-PPP(–) (Table S3), and they match very well those obtained originally by the Furchgott approach: 0.19, 0.12, and 0.05 (Meller et al., 1987). Finally, the two gain parameters obtained from two different data sets for this dopamine receptor system, γ=42.6±65.1 from the partial irreversible inhibition and γ=28.0±51.5 from the comparative activity of different agonists are consistent enough to support the applicability of the model here.

A second illustration is provided with an example that has been used to illustrate the capabilities of the τ-based operational model to explain unusual cases where full agonists become partial ones following receptor inactivation and even the order of apparent potencies changes (Kenakin, 1993; Kenakin and Christopoulos, 2011). Data are contraction of guinea-pig ileum following activation of muscarinic receptors with carbachol and oxotremorine in normal tissue and after the inactivation of muscarinic receptors by controlled alkylation with phenoxybenzamine (PHB) at two different strengths (10 μM for 10 min and 3 μM for 20 min) (Kenakin, 1993). With the present model, unified fit of all data can be achieved using a single γ for this pathway, single K _ds for carbachol and oxotremorine, and the same q values for the PBA-induced inactivation for both compounds ( Figure 12 , Table S4). This ε- and γ-based fitting can not only account for the apparent reversal of potencies, but even results in slightly better fit than that obtainable with the operational model and with more meaningful parameters (Buchwald, 2017). Fit indicates a strong amplification (≈21,000) that is then reduced more than a thousand-fold by alkylation and a 26-fold higher affinity (logK _d of –6.04 vs. –4.62), but a 56-fold less efficacy (ε of 0.018 vs. 1.0) for oxotremorine as compared to carbachol (Table S4). These K _d estimates are in very good agreement with those obtained for guinea pig ileum in a different work by the Furchgott method (–5.83 and –4.52 for carbachol and oxotremorine, respectively) (Eglen and Whiting, 1987) or measured for the M₃ receptor in another work (–5.61 and –4.09) (Sykes et al., 2009).

The present model might allow a better and more clearly parametrized approach. If it can fit the data, it allows the clear separation of pathway-specific differences in amplification from those in efficacies, which then can serve as cleaner indications of bias without a need for a reference agonist. If one assumes that the amplification for each signaling pathway is the same for all ligands (i.e., γ_k for pathway k, k = 1…l, is the same for all ligands L _i , i = 1…j) then biased agonism can be suspected for ligands that show significantly different efficacies for two (or more) pathways initiating from the same receptor (i.e., ε _im vs. ε _in for the efficacy of ligand i, L _i , for pathways m versus n). Hence, if receptor occupancy (binding affinity) data are available, one approach with the present model is to fit the data for each pathway with Equation 4 (Equation 2 if there is constitutive activity), and then compare efficacies between pathways. Bias is likely to be present for ligands with an efficacy ratio significantly different from 1. For cases where receptor occupancy data are available, the ability of the present model to connect (fractional) responses to (fractional) occupancy data (via Equation 4 or 27) provides a clear advantage, as it allows first an assessment of the adequacy of the model (i.e., whether unified gain parameters can provide adequate fit and well-defined efficacies for each response pathway) and then a detection of possible bias based on the obtained efficacy values. This is not possible with previous models, such as those based on the operational model, as they cannot incorporate experimental K _d values in their assessments. A modification of the operational model that uses experimental K _d values to constrain the regression has been proposed for bias detection by Rajagopal and co-workers (Rajagopal et al., 2011; Onaran et al., 2017). The bias quantification proposed here is in fact most similar to this one, since the εγ values of the present model will be very similar to the τ values obtained with this modified operational model, especially at large γ amplifications where εγ ≈ εγ − ε (cf. Equations 4 and 20). Under such conditions, comparing these τ ratios (Rajagopal et al., 2011; Onaran et al., 2017) is essentially the same as comparing ε ratios (as γ values are fixed for pathways). However, this operational model based approach cannot fit full agonists in systems with low amplification, where K _d and K _obs are similar (e.g., Figure 14B and D ), because full agonism can only be achieved with large τ’s so that K _obs = K _d/(τ + 1) will be shifted away from K _d. This limitation is avoided with the present parametrization.

If binding data are lacking, one can use fitted K _d ^’s enforcing a single set of values for all pathways. However, this way well-defined parameters are difficult to obtain, and simultaneous fit cannot be done for multiple pathways and ligands in commercial software such as GraphPad Prism. An alternative is to use the present model to directly fit the data in so-called bias plots (Gregory et al., 2010; Kenakin and Christopoulos, 2013) that are, in fact, relative response plots showing one response as a function of another one produced at the same ligand concentration (see, e.g., Figure 13 ). To derive the function directly connecting fractional responses f _resp1 and f _resp2 for the present model generated at the same ligand concentrations (hence, at the same fractional occupancies), f _occup will be expressed from f _resp for one of the pathways, and then used in the expression of the other. From Equation 27 (for pathway 2)

and its reverse (Equation 34) for pathway 1

one can link f _resp2 to f _resp1 directly by eliminating f _occup:

This results in various curvilinear relative response plots such as those shown in Figure 13 generated with the parameters as indicated [intended to resemble Figure 2 of the general review (Kenakin and Christopoulos, 2013)]. Note that for pathways that have different amplifications (γ ₁≠γ ₂), the relative response plots are still curvilinear even for ligands that do not show bias (i.e., balanced or “non-biased” agonists that have the same efficacy for both pathways, ε ₁=ε ₂):

Because of this nonlinearity, bias quantification is difficult especially for pathways that have considerably different amplifications. If the present model can fit the data, it allows the clear separation of pathway-specific differences in amplifications from those in ligand efficacies, which can then serve as indication of bias. Most biased agonism assessments published so far include only response measurements and no binding data; an example that includes affinity assessments too is discussed below for illustration. Due to the additional binding affinity data, such sets are better suited to judge the adequacy of the present model (i.e., whether unified gain parameters can provide adequate fit for each response pathway), and then detect actual bias based on the obtained efficacy values.

Biased signaling at the μ-opioid GPCR is of obvious interest for the possibility of achieving improved analgesia while reducing the unwanted side effects of opiate therapeutics, and oliceridine, a μ-opioid receptor (MOR) biased agonist, is one of the first products showing the clinical promise in developing biased agonists (Wadman, 2017). To illustrate the fitting of biased agonism data with the present model, we will use data from a recent detailed study with a series of ligands including morphine, endomorphin, met-enkephalin, DAMGO (L-tyrosyl-D-alanyl-glycyl-N-methyl-L-phenylalanyl-glycinol), and Pfizer standard-1 (2-(L-tyrosylamino-1-[N-acetly-L-phenylalanyl)-amino]-2-methylpropane hydrochloride) in cAMP and β-arrestin2 assays with the wild-type and two mutated μ-opioid receptors in which binding affinity estimates were also obtained using a [³H]diprenorphine competition assay (Hothersall et al., 2017).

Fit with the present model ( Figure 14A and B ) indicates that while there is essentially no amplification for the β-arrestin2 response (γ ₂ ≈ 1), there is a strong amplification for the cAMP response (γ ₁ = 62.0 ± 53.7) (Table S5). This is in agreement with the general observation that there is usually much stronger amplification in the G-protein mediated pathway than in the β-arrestin mediated one (Ehlert, 2018). Overall, quite good fits can be obtained accounting for 97% and 96% of the variability in the data for the cAMP and β-arrestin2 assays, respectively, even though the K _d-based model struggles somewhat to fit the data for DAMGO and met-enkephalin in the β-arrestin2 assay ( Figure 14B ). Part of the problem is more evident in the corresponding fractional response versus occupancy plot ( Figure 14D ), where it is clear that the f _resp vs. f _occup plot for these two compounds has an unusual upward curvature unsuitable to fit with the present model (Equation 27). Because of the large difference in the amplifications (>60-fold), the relative response (“bias”) plot is highly curved even for balanced (non-biased) agonists ( Figure 14E ). This makes the assessment of bias difficult, as even five-fold differences in the efficacies for activation of the receptor pathways (ε ₁ vs ε ₂) can be easily overshadowed by the differences in amplifications (γ ₁ vs γ ₂). Due to the amplification in cAMP, the corresponding ε values cannot be well determined and have significant uncertainties when fitting directly the relative response plot with Equation 38 ( Figure 14E ); hence, without binding information (K _ds), there is not enough data to determine bias in a well-defined manner. Including the binding data, and fitting with the K _d-based model either directly (Equation 4; Figure 14A and B ) or via the response vs occupancy data (Equation 27; Figure 14C and D ) gives better defined efficacies (Table S5). Bias calculated this way as the ratio of the two efficacies (ε _cAMP/ε _β-arrestin2) ranges from 0.281 ± 0.235 for DAMGO to 1.549 ± 0.525 for endomorphin-2. Of these, only that for DAMGO, a compound sometimes used as reference (Thompson et al., 2015), can be considered as indicating significant bias as it is the only one where the ratio ±2SD does not include 1, suggesting possible bias toward β-arrestin (Table S5). Note, however, that because of the somewhat poorer fit of the K _d-based model of the β-arrestin response for DAMGO, this conclusion should be treated with some caution. Further, in agreement with the original publication, this bias is not significant compared to Pfizer standard-1 as reference, since that compound also has a slight bias toward β-arrestin (0.783 ± 0.662).

A set of fractional receptor-G protein coupling (FRC) data assessed via a cell-free bioluminescence resonance energy transfer (BRET) assay for δ- and μ-opioid receptors (DOR, MOR) with different opioid compounds that share a peptidomimetic scaffold (Vezzi et al., 2013) was fitted with the present model assuming constitutive activity, but no amplification (ε _R0>0, γ=1; Figure 1B , Equation 19). This is included simply to illustrate the ability of the model as is to fit such data with non-zero baseline ( Figure 15 ) and to reproduce fitting with classic three-parameter (top, bottom, ED₅₀) models. In both cases, very good fits were obtained that account for more than 99% of the variability in the data (r ²=0.991) despite the wide range of ligand efficacies (Table S6). Furthermore, the fitted curves indeed fully overlapped with those obtained with a standard three-parameter log agonist vs. response model (E _max with adjustable bottom); they are both shown but are indistinguishable in Figure 15 . With the present model, ε _R0 values for each of the receptors were assumed to be shared values across the ligand series (as they are receptor characteristics) resulting in values of 0.464 ± 0.007 and 0.087 ± 0.008 for the δ- and μ-opioid receptors, respectively. Note that amplification may be present, but its magnitude cannot be determined from these data alone; hence, γ=1 was used. An illustration of the effects of different amplifications (γ) on the response generated by a receptor with constitutive activity with the present model (Equation 16) for a set of ligands of various efficacies (ε _R0, ε) is included in Figure S2.