A single unified model for fitting simple to complex receptor response data

The fitting of complex receptor-response data where fractional response and occupancy do not match is challenging. They encompass important cases including (a) the presence of “receptor reserve” and/or partial agonism, (b) multiple responses assessed at different vantage points along a pathway, (c) responses that are different along diverging downstream pathways (biased agonism), and (d) constitutive activity. For these, simple models such as the well-known Clark or Hill equations cannot be used. Those that can, such as the operational (Black&Leff) model, do not provide a unified approach, have multiple nonintuitive parameters that are challenging to fit in well-defined manner, have difficulties incorporating binding data, and cannot be reduced or connected to simpler forms. We have recently introduced a quantitative receptor model (SABRE) that includes parameters for Signal Amplification (γ), Binding affinity (Kd), Receptor activation Efficacy (ε), and constitutive activity (εR0). It provides a single equation to fit complex cases within a full two-state framework with the possibility of incorporating receptor occupancy data (i.e., experimental Kds). Simpler cases can be fit by using consecutively reduced forms obtained by constraining parameters to specific values, e.g., εR0 = 0: no constitutive activity, γ = 1: no amplification (Emax-type fitting), and ε = 1: no partial agonism (Clark equation). Here, a Hill-type extension is introduced (n ≠ 1), and simulated and experimental receptor-response data from simple to increasingly complex cases are fitted within the unified framework of SABRE with differently constrained parameters.

and its Hill-type extension cannot be used as they assume responses proportional with occupancy; hence, do not allow separation between fractional response, f resp = E /Emax , and occupancy, f occup = [LR occup ]/[LR max ] (Fig. 1A,B). They do not include parametrization for efficacy, only for occupancy via K d , the classic equilibrium dissociation constant characterizing receptor binding, which is defined in terms of the concentrations of the species involved (see Fig. 2, bottom row) and hence measured in units of concentrations: In most of these more complex cases, one can use the more cumbersome operational (Black & Leff) model (Eq. 4) 17 or its variations for fitting response data because it has an additional efficacy type parameter (τ): However, this model in its strict two-parameter version cannot connect response to occupancy because it cannot incorporate experimental binding affinities since its K D parameter is a fitted one that is not related to the experimental dissociation constant, K d . The discrepancy is especially serious for full or close to full agonists, where τ needs to have large values 16 . Related to this, an important limitation of the operational model is that most of the time, it is difficult to fit in a non-ambiguous manner, i.e., to obtain well-defined parameters. If just functional data are available (i.e., a single concentration-response curve), only the so-called transduction Top row: receptor response (red) and occupancy (blue) as a function of ligand concentration on typical semilog scales (f resp and f occup as a function of log C = log [L]). For the simplest cases (A-B), they overlap, and the red line is not visible. Bottom row: corresponding response versus occupancy curves (f resp as a function of f occup ). Note that in some of the more complex cases, the fractional response can be both ahead and behind the fractional occupancy even for the same compound (red and blue arrows indicating deviations from the unity line in E and F). While these are simulated data, experimental data illustrating such cases are also available (see, e.g., references [11][12][13][14] and discussions in 15,16   www.nature.com/scientificreports/ coefficient (τ/K D ) and not K D and τ independently can be estimated precisely due to identifiability issues during regression 18,19 .
To overcome these limitations, we have recently introduced a quantitative receptor model (SABRE) that can fit all these cases of increasing complexity with a unified equation using parameters for Signal Amplification (γ), Binding affinity (K d ), and Receptor activation Efficacy (ε) (plus constitutive activity, ε R0 , if needed) 15,16 . Its most general form that also includes a Hill coefficient (n), which will be derived here, is: As explicit incorporation of a signal amplification parameter is a main novelty, the model was designated SABRE as highlighted by the capitalization above. Not only is SABRE using parameters that are more intuitive and easier to interpret than those of the operational model or the del Castillo-Katz minimal two-state model, but, contrary to those, it can fit both simple and complex cases with the same equation (Eq. 5) with differently constrained parameters. It can be collapsed into consecutive simplified forms by fixing its parameters as special values (Fig. 2), and these can and should be used on their own when adequate. By introducing independent parametrization for the (post-receptor) signal amplification, SABRE allows a more clear conceptualization of receptor signaling as separate processes of binding, activation, and signal transduction (amplification) that can now be characterized and quantified via their own distinct parameters: K d , ε, and γ, respectively (Fig. 2).
Here, receptor response data of increasing complexity were fitted with SABRE using different levels of parameter constraining to illustrate the advantage of a unified model that allows nested simplifications (Fig. 2). Main details of model parametrization (including its Hill-type extension) are summarized in the Method section below followed by illustrative fittings of simulated and experimental response versus occupancy data in the Results section. These range from the simplest case, E max -type response only data, to complex multiple responses and biased agonism examples.

Methods
Model concepts. The present model maintains the main assumptions of the two-state receptor theory, e.g., that ligand-bound (occupied) and active receptor states do not fully correspond but introduces a slightly different and more intuitive parametrization-a detailed discussion is included in Ref. 16 . Briefly, binding of the ligand is assumed to alter the likelihood of activation: receptors can be active or inactive in both their ligand-free and ligand-bound forms (Fig. 2, top right); however, the corresponding probabilities (i.e., times spent in the respective conformations) can be quite different. Hence, ligand-free (R) and ligand-bound (LR) states are considered as an equilibrium ensemble of active ( * ) and inactive conformations present: R ⇌ R * and LR ⇌ LR * , respectively (not excluding the possibility that multiple, possibly overlapping active states might exist). In general, ligand-free receptors are overwhelmingly in their inactive conformation, R. In cases where there is no constitutive activity, they are entirely so. Binding of an agonist, which is governed by the affinity parameter K d , shifts the equilibrium toward the active state. The ability of a bound ligand to do so is characterized by an (intrinsic) efficacy parameter, ε. For receptors with constitutive activity, a basal receptor efficacy, ε R0 , is used to account for baseline activation even in absence of a ligand. The signal (response) generated by the active receptor (R * or LR * ) can be amplified downstream, and this is characterized by a pathway-specific gain parameter γ. Hence, the model uses four parameters: K d , the equilibrium dissociation constant characterizing binding affinity; ε, the efficacy characterizing the ability of bound ligand to activate the receptor (0 ≤ ε ≤ 1); ε R0 , a basal receptor efficacy characterizing constitutive activity (0 ≤ ε R0 ≤ 1); and γ, a gain (amplification) parameter characterizing the nonlinearity of (post-activation) signal transduction (1 ≤ γ). To accommodate responses that are more (or less) abrupt than those corresponding to a straightforward law of mass action, an additional a Hill-coefficient parameter, n, is also introduced here. Corresponding equations, schematics, and all successive nested simplifications leading to special cases are summarized in Fig. 2.
Parametrization. Occupancy (binding) parametrization, is achieved via a K d parameter similar to that of the original simple definition (Eq. 3), but with some modifications. SABRE differentiates between active and inactive receptor states (denoted by an asterisk, i.e., R * vs R and LR * vs LR), but considers ligand-bound and ligand-free states as an ensemble of conformations, so that K d represents an average binding constant for these ensembles of active and inactive forms that the ligand effectively sees. Hence, K d is defined in terms of the overall concentrations of occupied and unoccupied (ligand-free) receptors ( Fig. 2A) Accordingly, SABRE does not distinguish between binding affinities for the active and inactive states (as done by other two-state models, e.g., K d and K d /α; 16 ). It uses a single, ensemble-averaged K d that is a macroscopic equilibrium constant and, hence, experimentally measurable in equilibrium binding assays that assess total binding to the receptor. While this deviates from the assumptions of other models, the validity of such a single, experimentally measurable binding constant is nicely supported by previous works quantifying binding via both static (K d ) and kinetic (k off /k on ) methods in parallel with multiple downstream effects that found the experimental binding constants derived by these two methods to be in very close agreement (e.g., for the M 3 muscarinic receptor 12 and for the μ-opioid receptor 20 ).
Scientific RepoRtS | (2020) 10:13386 | https://doi.org/10.1038/s41598-020-70220-w www.nature.com/scientificreports/ Efficacy parametrization here is achieved via an ε parameter that represents the fraction of ligand-bound receptors that are active 15,16 : Hence, ε is a unitless parameter in the 0 to 1 range, and represents an intrinsic efficacy measured immediately post-receptor. Note, however, that it is not an equilibrium constant such as K ε = [LR * ]/[LR]). Constitutively active receptors are incorporated into the formalism of SABRE via a baseline receptor efficacy (ε R0 ) that is defined along similar lines, but for ligand-free receptors (i.e., the fraction of ligand-free receptors that are active): While ε is a ligand characteristic, ε R0 is a receptor characteristic. With these definitions, inverse agonists that reduce the signaling output below that of the basal state will have ε < ε R0 . Full agonists have ε = 1; partial agonists that generate a response but cannot reach the maximum one even at concentrations that saturate all receptor sites have ε R0 < ε < 1, and neutral antagonists have ε = ε R0 .
Finally, the present model explicitly incorporates pathway-specific signal transduction (amplification) via a separate gain parameter γ. Signal amplification has to be built into receptor models to account for cases where almost maximal responses can be achieved at relatively small fractional occupancies, i.e., cases that were traditionally designated as having "receptor reserve" or "spare receptors". In SABRE, the fraction of active receptors is linked to the fractional response, f resp = E /Emax , via a hyperbolic function. Such functions provide convenient ways to incorporate signal amplification cascades and have been shown to exist in response vs occupancy data. However, in SABRE not f act itself, but its odds-ratio type transform serves as input (see Ref. 15

for details):
With this definition, γ represents a unitless amplification (gain) factor (γ ≥ 1). After some transformations, this results in the final form of the full four-parameter model linking fractional response, f resp = E/E max , to ligand concentration [L]: For cases with no constitutive activity (ε R0 = 0; no active unbound receptor, R * ), this reduces to the threeparameter minimal two-state model previously introduced 15 (third row of Fig. 2) and the simplified equation: A detailed derivation of these equations, including its generalization for a Hill-type extension, is included in Supporting Information, Appendix 1. A better interpretation of these parameters can be gleaned from a slightly rearranged form of this three-parameter equation, which corresponds to a case with Hill coefficient slope of unity (n = 1): From here, it is clear that half-maximal response (EC 50 ) is observed at K obs = K d /(εγ-ε + 1) and the maximum (fractional) effect achievable by a given ligand is f resp,max = εγ/(εγ-ε + 1). For full agonists at the receptor (ε = 1), K obs = K d /γ; therefore, the gain γ is a straightforward multiplication factor causing a left shift of the sigmoid response function by γ units on a semi-log scale. By explicitly separating the parametrization of pathway-specific amplification (γ) from that of ligand-specific receptor activation (ε), SABRE more clearly outlines than other models the intrinsic efficacy concept, which proved to be somewhat elusive in pharmacology 21,22 . Hill-type extension (cooperative binding). As a final step, a Hill-type parametrization can also be introduced via a Hill slope, n, to allow either more (n > 1) or less (n < 1) abrupt concentration-dependent responses: Scientific RepoRtS | (2020) 10:13386 | https://doi.org/10.1038/s41598-020-70220-w www.nature.com/scientificreports/ This function, originally introduced by Hill based on empirical considerations 23 , provides a versatile mathematical function often used in pharmacological 24,25 and other applications 26 . It shows analogy with the logistic function, one of the most widely used sigmoid functional forms, as it is equivalent with a logarithmic logistic function, y = f(x) = R max /(1 + βe -nlnx ) 27 . In fact, the classic sigmoid shapes obtained in typical semi-log scale graphs (such as those shown in Fig. 1) are those of the logistic function. The IUPHAR recommendation is to use "Hill equation" for the relationship between ligand concentration and effect (such as Eq. 15 here), whereas "Hill-Langmuir equation" should be used for relationship between ligand concentration and occupancy 27 . Hill slopes different from unity are typically indications of interacting binding sites with positive (n > 1) or negative (n < 1) cooperativity 28 . The Clark equation (Eq. 1), as well as the analogous Michaelis-Menten equation, represent a special case (n = 1) of the Hill equation-a nice example of how more complex models can be collapsed into simplified forms for special cases of their parameters.
Following a detailed derivation (see Supporting Information, Appendix 1), the Hill-type extension of the general form of the present SABRE model that includes constitutive activity is a straightforward generalization of Eq. 12 introducing the Hill-coefficients as exponents of the concentration terms: In case of no constitutive activity, this simplifies to the Hill-type extension of Eq. 13: Consecutive simplifications. An important feature of the present model is that, contrary to previous quantitative receptor models such as those based on the operational model, its general form (Eq. 5) can be reduced to consecutively nested simplified forms for special cases of its parameters (Fig. 2), and these can be used on their own when adequate. Hence, the very same model can be used to fit data of various complexity levels using different sets of constrained parameters. First, by setting n = 1, it reproduces the simpler form of SABRE as introduced before 16 , just as setting n = 1 in the Hill equation (Eq. 2 or 15) reproduces the Clark equation (Eq. 1) as its simpler form.
The way the Hill equation can be reduced to the Clark equation by constraining one of its parameters to a special value (n = 1) has been a long-known example in quantitative pharmacology. With SABRE, one can do the same and more, as one can continue with multiple consecutive simplifications of its parameters. Cases with no constitutive activity (no R * form) can be obtained by constraining ε R0 to 0 leading to the three-parameter minimal two-state SABRE model that was the first form introduced 15 : Further, if there is no amplification (or if it cannot be reliably evaluated due to lack of independently assessed occupancy/binding data), γ should be constrained to 1, and the model collapses further into the simple twoparameter E max model of partial agonism (Fig. 2): Finally, if there is no partial agonism (i.e., all occupied receptors are active), ε can be constrained to 1, and this leads to the simple one-parameter Clark equation (Fig. 2, bottom row): Regarding these consecutively nested simplified forms, it is important to always use the simplest one (i.e., the one with the largest number of fixed parameters) that still provides adequate fit to limit the number of adjustable parameters [29][30][31][32] . Adequate fitting requires the availability of 5-10 (well-distributed) data points for each adjustable parameter 5,33,34 ; hence, reliable fitting of the full model can only be accomplished if sufficiently large number of data points are available. Because of its separate amplification parameter, SABRE uses one more parameter than the operational model, e.g., three (K d , ε, γ) versus two (K d , τ) for the case of no constitutive activity (Eq. 19 vs. 4) or four (K d , ε, γ, ε R0 ; Eq. 18) versus three (e.g., K d , ε, χ; Eq. 22) for cases with constitutive activity-if comparing with the Slack and Hall version of the extended operational model shown below [35][36][37] : (2020) 10:13386 | https://doi.org/10.1038/s41598-020-70220-w www.nature.com/scientificreports/ However, the need for one extra parameter is more than compensated for by (1) the more intuitive nature of the present parameters (due to their more clear connection to receptor binding, activation efficacy, and signal amplification), and (2) the ability to use simplified forms with constrained parameters so that fewer need fitting. Contrary to the operational model and its variations, SABRE can be reduced to simplified forms for special cases of its parameters, and these can and should be used on their own when adequate or when there are not enough data to support full parametrization. Implementation and data fitting. All data are normalized and have no baseline (i.e., in the 0-100% range). Data fittings were performed in GraphPad Prism (GraphPad, La Jolla, CA, USA, RRID:SCR_002798). Fitting with the present model were done with a custom implementation corresponding to the general Eq. 5, which is available for download (Supporting Information) and with parameters constrained as indicated for each case. Simulated data were generated with the same model in Prism using the "Simulate XY data" algorithm with 5% random error). Experimental data used for illustrations of model fit are reproduced from previous publications as indicated in the corresponding figures.

Results and discussion
In systems with signal amplification, ligands of different efficacies (e.g., full, partial, and possibly inverse agonist) can produce complex concentration-response functions leading to complicated connections between fractional response (f resp ) and occupancy (f occup ) ( Fig. 1) that can be fit only by multi-parameter models. So far, no single quantitative receptor model could fit simple as well as complex cases within a unified framework. Here, it is shown how this can be done with the present model, which allows consecutive nested simplifications, using a set of examples of increasing complexity. Simulated as well as experimental data were fitted with the same general equation (Eq. 5), but with different levels of parameter constraining. To avoid excessive parametrization, all cases assume n = 1, i.e., no Hill-type extension (except for one specific illustration in Fig. 4), as well as ε R0 = 0 (no constitutive activity).

Response only data (E max model with possible Hill-type extension).
As the simplest first case, let us consider the fitting of plain response data to determine EC 50 (K obs ) values for agonists of different potencies in a given assay system. EC 50 is the half-maximal (or median) effective concentration (or dose, ED 50 ), which is the concentration (dose) of an agonist that produces 50% of the maximum possible activity of that agonist 27 . Such EC 50 determinations are routinely done in pharmacology using either the one-parameter Clark equation (if only full agonists are present) or the two-parameter E max equation (if partial agonists are also present). These models are widely used, familiar to those performing such nonlinear regression-based fit, and implemented in a variety of software packages, e.g., "log(agonist) vs. normalized response" and "log(agonist) vs. response (three parameters)" models in GraphPad Prism. The same fit can be obtained with SABRE (Eq. 5) by constraining all but one (K d ) or two (K d , ε) of its parameters. An illustrative example using simulated data for three different hypothetical agonists is shown in Fig. 3. Fit with SABRE reproduced perfectly both the log K d values (which here correspond to log EC 50 s; − 6, − 7, and − 8) and maximum fractional response (efficacy) values (100%, 90%, and 60%) used to generate the data. To have realistic data with some scatter, a 5% random error was allowed. Hence by restricting its parameters, SABRE can be reduced to reproduce simpler forms that can be used on their own, something that cannot be done with models based on the operational model.
If response data are either more or less abrupt than predicted by the standard law of mass action (unity Hill slope, n = 1), the Hill coefficient of Eq. 5 also needs to be released resulting in three-parameter fit (K d , ε, n). An illustration of this is shown in Fig. 4 with data generated as in Fig. 3 but with Hill slopes of 2.0 and 0.66. In most cases, n should be constrained to the same value for all compounds assessed in the same assay; here, this was not done to illustrate the effect of both n > 1 and n < 1.
Connecting response to independently measured occupancy data. Next, as a more complex case, let us consider fitting of the same type of response data (EC 50 ), but with integration of occupancy (K d ) data obtained from a different, independent assay in the same system. As an example, the same simulated response data as above was used ( Fig. 3; log EC 50 s of − 6, − 7, and − 8), but with an additional set of K d s (logK d s of − 5.2, − 6.6, and − 6.7). Fit with SABRE can be considered consistent if adequate fit can be obtained with all K d s restricted to their experimental value and a single amplification parameter γ characterizing the system (pathway) as a whole. The corresponding fit for the present data is shown in Fig. 5A. Very good fit could be obtained with only four adjustable parameters (3ε + 1γ): > 98% of variability could be accounted for with a well-defined gain parameter: r 2 = 0.984, γ = 21.1 ± 4.7. While this was done here with simulated data customized to illustrate such fit, examples of experimental data that can be fit with SABRE using a unified amplification parameter are also available. For example, contraction of isolated rat aorta data obtained for imidazoline type α-adrenoceptor agonists 11 , which are used as textbook illustration for possible mismatch between f resp and f occup 38 , could be fitted very well with SABRE with the assumption of a single gain (amplification) parameter for data from five compounds (γ = 11.9 ± 2.0; r 2 = 0.996; see Fig. 9 and Table S1 in Ref. 16 ).
Further, within the framework of SABRE one can not only fit such concentration-dependency data (Fig. 5A), but also directly link fractional response, f resp , to fractional occupancy, f occup 16 (Fig. 5B). By using f occup to replace [L] (see Supporting Information, Appendix 1), Scientific RepoRtS | (2020) 10:13386 | https://doi.org/10.1038/s41598-020-70220-w www.nature.com/scientificreports/ one gets the corresponding function that can be used for direct fitting: Typically, binding data are used to generate K d estimates (usually via a Cheng-Prusoff correction 39 ). If these K d s are then used to calculate occupancy (f occup ) at the [L] values were response data (f resp ) are available, fit with the above equation provides the same parameters as obtained from direct fit of the concentration-response data. Illustration for the present data is provided in Fig. 5B. While models based on the two-parameter operational model (Eq. 4) can be used to fit the response data (with some known problems fitting full and close to full agonists), they have difficulties connecting the response to independently determined binding data (K d ) as done here with SABRE. One exception is a three-parameter "special edition" extension of the operational model "with given K d values" introduced by Rajagopal and Onaran for bias quantification that uses experimental K d values as its K D 19,40 . In this model, K D s are not fitted, but replaced with experimental ones to constrain the regression. To achieve this, an additional scaling factor α (r max ) needs to be introduced in Eq. 4 to allow a scalable E max : The three parameters of this model, however, can be linked directly to those of the present model (see Supporting Information, Appendix 2) so that parameters obtained from fitting SABRE can be used to derive those of this model (K d , ε, γ → K D , τ, α): Figure 3. Fit of response only (EC 50 ) data (full and partial agonists) with the present SABRE model. Simulated data (symbols) for three different compounds were generated in GraphPad Prism ("Simulate XY data" algorithm with 5% random error) and then fitted with SABRE (lines) using the constrained version as shown (Eq. 5; ε R0 = 0, γ = 1, n = 1). This is equivalent with fit with the E max model, e.g., the widely used "log(agonist) vs. response (three parameters)" model, and results in exactly the same parameters (assuming a 0 baseline).

Scientific RepoRtS |
(2020) 10:13386 | https://doi.org/10.1038/s41598-020-70220-w www.nature.com/scientificreports/ Conversely, the parameters of this "special edition" operational model can be transformed into those of SABRE using (Appendix 2): This highlights again an important advantage of the present model, namely that it clearly separates the ability of the ligand to activate the receptor (ε) from the pathway-specific (post-receptor) signal amplification (γ), while these two are intermixed in the τ "transducer ratio" ("coupling efficiency") parameter of the operational model and its different variants [i.e., τ = ε(γ-1) here; Eq. 26]. Hence, from the perspective of the present framework, the τ parameter of the operational model mixes together ligand-and pathway-specific effects separated into ε and γ here (see Appendix 2). Note also that this "special edition" operational model (Eq. 25) is somewhat unusual as it formally allows an E max larger than 100% by having α = γ/(γ-1) > 1.
Responses assessed at different vantage points. Next, let us consider the case of multiple responses measured at different downstream vantage points along the same pathway. This is of particular relevance for G-protein coupled receptors (GPCRs), known to involve signaling cascades with multiple second messengers. A set of such data generated within the framework of the present assumptions is shown in Fig. 6; they were chosen so as to illustrate a possible nonintuitive case (shift in the order of potencies) resulting from the intermingling of the effects caused by different efficacies at the receptor on one hand and increasing downstream signal amplification on the other. Simulated responses were generated for three hypothetical compounds of different affinities (log K d of − 6, − 7, and − 5) and efficacies (ε of 0.5, 0.1, and 1.0) at three consecutive readout points of increasing amplifications (γ of 1, 20, and 500). At all three points, responses are shown as a function of both concentration (f resp vs log C) and occupancy (f resp vs f occup ) (Fig. 6). Despite the unusual nature of the data, SABRE can provide integrative fit with a unified set of parameters (n p = 9; 2 × 3 for the binding affinities and efficacies of the three www.nature.com/scientificreports/ agonists plus 3 for the amplifications of the three readout assays), not conceivable with operational model-based approaches that could only provide separate sets of K D s and τs for each individual response (i.e., n p = 18; 2 × 3 parameters of the three agonists for each of the 3 readout assays). A further illustration is provided by fitting of a set of experimental data obtained from such consecutive responses (Fig. 7). They were measured at two consecutive vantage points after M 3 muscarinic receptor activation by agonists such as carbachol, oxotremorine, oxotremorine-M, and methacholine: the stimulation of GTP binding to Gα subunits and the subsequent increase in intracellular calcium levels 12 . Receptor occupancy estimates are from log K d values obtained from equilibrium competition experiments with N-methyl-[ 3 H]scopolamine in the same work 12 . Unified fit of n d = 72 data points for both responses from four compounds using SABRE with a single set of ε efficacies (one for each compound) and two gain parameters γ (one for each response) as adjustable parameters (n p = 6) resulted in reasonable amplifications estimates (γ of 2.31 ± 0.42 and 10,243 ± 1,351 for the GTP and Ca readouts, respectively) and good overall fit (accounting for 97.8% of the variability in the data, r 2 = 0.978) (Fig. 7). Because of a very strong (~ 10,000-fold) amplification in the second (Ca response) pathway, efficacy estimates are mainly defined by the first (GTP binding) response-obtained values and standard errors are summarized in Fig. 7. Such fittings can be considered consistent only if the unified single set of efficacies can provide acceptable fit at all the response levels considered for each compound including a full agonist. Nevertheless, unlike any other previous model, SABRE has the potential to connect response data assessed at different vantage points k (E k /E k,max ) to affinities (log K d ) and intrinsic efficacies (ε) for multiple compounds as long as reasonable overall fit can be obtained with unified gain parameters (γ k ).
Biased agonism. Finally, let us consider the application of the present model for quantification of biased agonism. It has been recognized since at least the 1990s that some receptors can engage multiple downstream Multiple response data at different vantage points generated within the framework of the present SABRE model. Simulated data (symbols) for three different compounds were generated with the parameter values as shown (three compounds with log K d s of − 6, − 7, and − 5 and εs of 0.5, 0.1, and 1.0 at three different readouts with γs of 1, 20, and 500) and 5% random error; lines indicate corresponding model fits obtained with SABRE using a unified set of parameters (n p = 9). Responses at the three different readout levels are shown as a function of both concentration (f resp vs log C; middle column) and occupancy (f resp vs f occup ; rightmost column). Data were chosen so as to illustrate a possible nonintuitive shift in the order of potencies that can result from the strong downstream amplification of the response caused by a weak agonist: compound 2 (red; ε = 0.1) that generates the weakest response right after the receptor (assay 1, top) becomes the most potent one in assay 3 (bottom).
An illustration of such data generated within the framework of the present assumptions for two divergent pathways involving different signal amplifications and three hypothetical compounds of different efficacies is shown in Fig. 8. Corresponding data are shown as typical concentration-response curves (f resp vs log C), response versus occupancy ones (f resp vs f occup ), and a typical bias plot used in such cases (i.e., relative response plot of f resp1 vs f resp2 44,54 ). Data were generated for two different pathways (P 1 , P 2 ) originating from the same receptor, but having different amplifications (γ P1 = 4 and γ P2 = 20) and three hypothetical compounds having different affinities (log K d s of − 6, − 7, and − 5) and efficacies ε. Two compounds (1 and 3) had the same efficacy for both pathways (ε 1,P1 = ε 1,P2 = 0.5 and ε 3,P1 = ε 3,P2 = 1.0), while one (2) had a higher efficacy for pathway 1 (ε 2,P1 = 0.5, ε 2,P2 = 0.1) generating a biased response as most evident in the bias plot of Fig. 8 (bottom row, center).
Because of the intermix of effects related to binding, receptor activation, and signal amplification, quantifying signaling bias is difficult, and it may not be achievable in many cases 19,55 . Current bias quantification methods typically rely on calculating ΔΔlog(τ/K D ) or ΔΔlog(E max /EC 50 ) versus a selected reference compound, e.g., a logarithmic bias factor is obtained as (log(E max,P1,L /EC 50,P1,L )-log(E max,P2,L /EC 50,P2,L )) − ((log(E max,P1,Lref / EC 50,P1,Lref )-log(E max,P2,Lref /EC 50,P2,Lref )) 19,49 . This originated, in fact, from the concept of ratio of equiactive molar ratios, later termed intrinsic relative activity (RA i ) introduced by Ehlert and co-workers, which is the ratio of τ/K D fractions with the operational model and becomes that of E max /EC 50 s for a Hill slope of 1, i.e., (E max,L /EC 50,L ) / (E max,Lref /EC 50,Lref ) 56,57 . Because SABRE explicitly separates pathway amplification from ligand efficacy, it might allow a conceptually different approach, as discussed briefly before 16 . If data can be fitted with sufficiently welldefined gain parameters (γ Pk , one for each pathway P k ), then pathway-specific differences in amplification can be separated from those in ligand-specific efficacies, and these εs can serve as cleaner indicators of bias, possibly even without predefined reference agonists. If receptor occupancy (K d ) data are available, and SABRE can fit the data for each pathway adequately, one can calculate ligand efficacies for each pathway separately and compare them for indication of bias. If fitting can be done so that full or close to full agonists (ε = 1) are identified for each pathway, ligands that have efficacy ratios (ε Pk /ε Pl ) significantly different from 1 can be considered as being biased agonists. Otherwise, εγ products need to be compared. There are methods to estimate efficacies with two-state models 58 ; nevertheless, because SABRE is the first model that explicitly uses separate parameters for efficacy (ε) and amplification (γ), it can untangle these connections in a manner not possible either with direct empirical comparisons (e.g., E max /EC 50 ) or with fitting of previous models, which intermingled these two effects within the same parameter (e.g., τ, χ).
To illustrate the process of bias detection with the present model, fit was done on recent experimental data obtained at the angiotensin II type 1 receptor (AT1R): two responses (Gq-mediated inositol monophosphate  12 ). Experimental data are for carbachol (Cpd1, blue), oxotremorine-M (Cpd2, red), oxotremorine (Cpd3, green), and methacholine (Cpd4, purple). Fitting of n d = 72 total data points with SABRE (lines) was done using the experimental K d values (indicated on the x axis) and optimizing only n p = 6 adjustable parameters: two γs (one for each pathway) and four εs (one for each test compound) as shown.
If binding data are lacking (no experimental K d s), which is the case for most biased agonism studies published so far, one can still use the present approach, but with fitted K d s. This should be done by enforcing a single unified K d for each ligand to minimize the number of adjustable parameters n p . Nevertheless, even with this restriction, due to the more limited amount and nature of data, it typically becomes difficult to obtain well-defined parameters and a clear bias quantification. For example, for the data of Fig. 9, the fitting becomes ambiguous and none of the amplification or efficacy parameters can be obtained with well-defined values. The efficacy ratios for TRV026 and TRV023 remain about the same as before, but the standard errors are much wider than the values themselves (13.9 ± 15,285 and 17.6 ± 19,365). This is not surprising as the number of data points per adjustable parameters, n d /n p , here is only 4.1 (45/11) compared to 5.6 (45/8) before, which was much closer to the desired range of 5-10 needed for adequate fitting 5,33,34 . Accordingly, to improve the quantitative assessment of multiple responses and biased agonism, it is suggested to include experimental assessment of binding in the same system where response is measured whenever possible.
Finally, if needed, within the framework of SABRE one can directly fit bias plot 44,54 (i.e., relative response) data that show one response as a function of another at the same ligand concentration (e.g., Figs. 8D and 9E). The function directly connecting fractional responses f resp,P1 and f resp,P2 generated along different pathways P k at the same ligand concentrations [L] (hence, at the same f occup ) can be expressed as 16 : Figure 8. Illustration of biased agonism with response data for two different downstream pathways originating from the same receptor (A) generated with the assumption of the present SABRE model. Simulated data (symbols) for three hypothetical compounds were generated as before (Fig. 6) using the parameter values shown (CpdTst2 having two different efficacies ε for pathways P 1 and P 2 as highlighted in yellow; see text for details). Data for two pathways involving different signal amplifications (γ P1 = 4 left and γ P2 = 20 right) are shown as classic semi-log concentration-response curves (f resp vs log C; B1, B2), fractional response versus occupancy curves (f resp vs f occup ; C1, C2), and a bias plot (f resp1 vs f resp2 ; D).
Scientific RepoRtS | (2020) 10:13386 | https://doi.org/10.1038/s41598-020-70220-w www.nature.com/scientificreports/ Figure 9. Fit of experimental response data obtained from divergent pathways with the present SABRE model. Data are for two responses measured along different downstream pathways (Gq-mediated inositol monophosphate increases and β-arrestin2 endocytosis) at the angiotensin II type 1 receptor (AT1R) (data after 53 ). Experimental data are for three agonists (angiotensin II, Cpd1, blue), TRV023 (Cpd2, red), and TRV026 (Cpd3, green) as indicated. (A) Fractional response versus log concentration data shown in a single composite graph. Fitting of n d = 45 total data points with SABRE was done using the experimental K d values and optimizing only n p = 8 adjustable parameters (2 γs + 3 × 2 εs) as shown. (B) Same response versus log concentration data, but shown in two separate graphs, one for each response. (C) Fractional occupancy versus log concentration data. The concentration dependency shown was calculated from the experimental K d values with standard Hill-Langmuir equation. (D) Fractional response vs occupancy data. Same data as in A, just shown as a function of calculated fractional receptor occupancy, not log concentration. (E) Bias (relative response) plot showing one fractional response as a function of the other. Compound 1 (AngII) is a balanced (nonbiased) agonist, the curvature is due to the slightly different amplifications for the two responses. Cpd2 and 3 are significantly biased and they elicit a stronger response along pathway 1 (β-arrestin) than pathway 2 (Gq). Scientific RepoRtS | (2020) 10:13386 | https://doi.org/10.1038/s41598-020-70220-w www.nature.com/scientificreports/ Accordingly, pathways that have different amplifications (γ P1 ≠ γ P2 ) lead to relative response plots that are curvilinear even for balanced ligands L i that do not show bias (i.e., ε i,P1 = ε i,P2 ): While K d s are not needed for this type of comparison, quantitative fitting (Eq. 28) is likely to be ambiguous due to lack of sufficient data points per adjustable parameter. For example, for the case shown in Fig. 9, n d /n p is only 21/8 = 2.6 resulting in very badly defined ε ratios. Nevertheless, bias plots can still be useful as they might allow identification of differences among ligands that are less evident in concentration-response or response vs occupancy plots (e.g., Fig. 8). However, both responses might be difficult to obtain at the same ligand concentration [L] and strong curvatures might confound assessments, especially in cases where pathways have considerably different amplifications.

Conclusion
In conclusion, SABRE is a general two-state ensemble receptor model incorporating a corresponding quantitative form that can fit complex fractional response versus occupancy data. Contrary to previous models, it not only has more intuitive parameters that can be related directly to binding affinity (K d ), activation efficacy (ε), and signal amplification/gain (γ), but it can incorporate experimental K d values. Notably, it provides a unified framework to fit both complex and simple receptor data, as by constraining its parameters, the general equation (Eq. 5) can be converted into consecutively nested simplified models all the way down to the one-parameter Clark equation (Fig. 2).

Data availability
Data used for illustrations of model fit are either simulated data generated as described or reproduced from previous publications as indicated in the corresponding figures. The datasets generated and/or analyzed and a GraphPad Prism file with the implementation of the model discussed here are available from the corresponding author upon reasonable requests.