Introduction

In pharmacology, efficacy serves as a measure of how much response each ligand-receptor complex can produce upon its formation1. A sound method to assess agonist efficacy is essential for research and drug discovery. The operational model of agonism (OMA) introduces the term operational efficacy2. OMA has become a golden standard in the evaluation of agonism and subsequently also of signalling bias3,4,5. The OMA describes the response of the system as a function of ligand concentration using three parameters: (1) The equilibrium dissociation constant of agonist (KA) to the receptor initiating functional response; (2) The maximal possible response of the system (EMAX); (3) The operational efficacy of agonist (τ). As we will show, in the OMA, KA represents the affinity of the agonist for the receptor. Therefore, KA is specific to a given combination of ligand and receptor. The maximal possible response of the system EMAX is specific to the system. The operational efficacy (τ) is a measure of the response to a given agonist at a given system, ranging from 0 to infinity, and is specific to a combination of ligand and system.

OMA has numerous limitations and possible pitfalls. All three parameters (EMAX, KA, and τ) are inter-dependent, thus, one of them has to be predetermined before fitting OMA to data or the global fit of multiple curves at various receptor densities has to be performed6. Besides that, there are additional methodological and conceptual issues7,8,9. Here we focus on the exponentiation of operational efficacy. In practice, positive cooperativity or positive feedback leads to steep and negative cooperativity or negative feedback leads to flat concentration–response curves that the OMA (Eq. (2)) does not fit. Therefore, an equation intended for the description of non-hyperbolic concentration curves (Eq. (7)) was introduced by Black et al.10. Since then, Eq. (7) has been commonly used. The presented simple mathematical analysis of Eq. (7) shows that the slope factor affects the relationship between observed maximal response to agonist (E’MAX) and operational efficacy (τ) and the relationship between the concentration of agonist for half-maximal response (EC50) and its equilibrium dissociation constant (KA). Some combinations of operational efficacy and slope factor in the Black & Leff equation (Eq. (7)) lead to EC50 greater than KA which is biochemically counterintuitive. In the system with receptor reserve maximum effect is reached before the occupation of all receptors is reached. Therefore, half of the effect (EC50) is reached before reaching the occupation of half of all receptors (KA). In a system without receptor reserve, full occupation of receptors is reached before reaching system EMAX. Therefore, EC50 occurs at KA. Black & Leff model, however, allows for EC50 > KA in the systems with low efficacy and high slope factor.

The Hill equation was originally formulated to describe the binding of oxygen molecules to haemoglobin11,12. The Hill equation was then incorporated into the first model of receptor function, the so-called occupation theory13. In contrast to the Black et al. slope factor, the slope factor of the Hill equation (Eq. (10)) does not affect the centre or asymptotes of the hyperbola describing concentration–response curves. Therefore, the Hill coefficient does not affect relationships between E’MAX and τ nor between EC50 and KA. This makes Hill equation-based OMA (Eq. (11)) more practical in many ways. Fitting the Black & Leff equation (Eq. (7)) to the theoretical data revealed several drawbacks, like under- or over-estimation of parameters or high levels of uncertainty of parameter estimates. Moreover, in contrast to Hill equation-based OMA (Eq. (11)), in some cases, fitting the Black & Leff equation to the experimental data resulted in the wrong ranking of agonist efficacies τ and wrong estimates of agonist KA. However, it should be noted that the Black & Leff equation should not be used firsthandily as slopes of individual response curves vary among agonists.

The general concept of the operational model of agonism

For demonstrative purposes, we will derive OMA from scratch for receptor-effector systems. This will become handy for the analysis of systems with low expression of receptors. The OMA equation for ligand-gated ion channels is the same, although, it is based on different sets of equations. In general, OMA consists of two functions. One function describes the binding of an agonist to a receptor as the dependence of the concentration of agonist-receptor complexes [RA] on the concentration of an agonist [A]. The second function describes the dependence of functional response (E) on the concentration of agonist-receptor complexes [RA]. OMA expresses the dependence of response E on the concentration of [A].

Rectangular hyperbolic OMA

Definition of OMA

In the simplest case, when both binding and response functions are described by rectangular hyperbola (Supplementary Information Figure S1), the resulting function is also a rectangular hyperbola. For example, in a bi-molecular reaction, the dependence of ligand binding to the receptor [RA] is described by Eq. (1) where [A] is the concentration of ligand and KA is its equilibrium dissociation constant that represents the concentration of ligand at which half of the total number of receptors, RT, binds the ligand and the other half of the receptors is free.

$$\left[ {RA} \right] = \frac{{R_{T} \left[ A \right]}}{{\left[ A \right] + K_{A} }}$$
(1)

If the bound ligand is an agonist, it activates the receptor and produces functional response E. Response as a function of agonist binding (agonist-receptor complexes [RA]) is given by Eq. (2).

$$E = \frac{{E_{MAX} \left[ {RA} \right]}}{{\left[ {RA} \right] + K_{E} }}$$
(2)

where EMAX is the maximum possible response of the system and KE is the value of [RA] that elicits a half-maximal effect. Various agonists produce a functional response of different strengths. The OMA was postulated to introduce the “transducer ratio” τ that is given by Eq. (3).

$$\tau = \frac{{R_{T} }}{{K_{E} }}$$
(3)

The substitution of Eq. (2) with Eq. (1) and Eq. (3) gives Eq. (4).

$$E = E_{MAX} \frac{\left[ A \right]\tau }{{\left[ A \right]\left( {\tau + 1} \right) + K_{A} }}$$
(4)

Analysis of OMA

Equation (4) is the equation of OMA2. It has three parameters: The equilibrium dissociation constant of agonist (KA) at the effect-eliciting (active) state of the receptor6 that is specific to a combination of ligand and receptor. The maximal possible response of the system (EMAX) is specific to the system. And the “transducer ratio” (τ) that is specific to a combination of ligand and system. Equation (4) is a rectangular hyperbola with the horizontal asymptote, the observed maximal response to agonist A (E’MAX), given by Eq. (5).

$$E_{MAX}^{\prime } = E_{MAX} \frac{\tau }{\tau + 1}$$
(5)

A more efficacious agonist (having a high value of parameter τ) elicits higher E’MAX than less efficacious agonists (having a low value of parameter τ). Thus, τ is actually operational efficacy. The relationship between parameter τ and E’MAX is hyperbolic meaning that two highly efficacious agonists (e.g., τ values 10 and 20) differ in E’MAX values less than two weak agonists (e.g., τ values 0.1 and 0.2).

In Eq. (4), the concentration of agonist A for half-maximal response (EC50), is given by Eq. (6).

$${\text{EC}}_{50} = { }\frac{{{\text{K}}_{{\text{A}}} }}{{{\uptau } + 1}}$$
(6)

According to Eq. (6), for τ > 0, the EC50 value is always lower than the KA value. The KA to EC50 ratio is greater for efficacious agonists than for weak agonists. Similarly to E’MAX, the relationship between parameter τ and EC50 is hyperbolic. In contrast to E’MAX values, the ratio KA to EC50 ratio is more profound for two highly efficacious agonists (e.g., τ values 10 and 20) than for two weak agonists (e.g., τ values 0.1 and 0.2).

Limitations of OMA

The OMA has several weak points. The major drawback of OMA is the lack of physical basis of the agonist equilibrium dissociation constant KA. Equation (2) assumes that agonist binding [RA] denotes effect-producing active complexes. The agonist binding to the receptor in an inactive conformation is not observed in the response (KE → ∞; τ = 0). In the radioligand binding experiments, agonists bind to all receptor conformations including the inactive ones. For various reasons, the receptors in an active conformation may be scarce or absent from radioligand binding experiments. Then it may be impossible to determine the KA value in the radioligand binding experiments. All three parameters of OMA (EMAX, KA and τ) are interdependent6. To fit Eq. (4) to individual concentration–response curves, one of the parameters must be fixed. E.g., the maximal response of the system EMAX is determined by comparing the functional response to a given agonist in a system with a reduced population of receptors by irreversible alkylation14 or cell lines with varying receptor expression levels6. Alternatively, global fitting of multiple concentration–response curves with all parameters free can be employed. However, due to a high number of degrees of freedom global fitting is less robust than per partes methods6. Another limitation of the OMA is that the shape of the functional response is a rectangular hyperbola.

Non-hyperbolic OMA

Definition of non-hyperbolic OMA

In practice, concentration–response curves steeper or flatter than the ones described by Eq. (4) are observed. In such cases, Eq. (4) does not fit experimental data. As stated by the authors, Eq. (7) was devised for non-hyperbolic dependence of functional response on the concentration of agonist10. Eq. (7) was derived in the same way as Eq. (4) from Eq. (1) and Eq. (2) while [RA] and KE in Eq. (2) were exponentiated to factor n.

$${\text{E}} = { }\frac{{\left[ {\text{A}} \right]^{{\text{n}}} {\uptau }^{{\text{n}}} {\text{E}}_{{{\text{MAX}}}} }}{{\left[ {\text{A}} \right]^{{\text{n}}} {\uptau }^{{\text{n}}} + { }\left( {\left[ {\text{A}} \right] + {\text{ K}}_{{\text{A}}} } \right)^{{\text{n}}} }}$$
(7)

Analysis of non-hyperbolic OMA

Introduced power factor n changes the slope and shape of the functional-response curve (Supplementary information Figure S3). Nevertheless, Eq. (7) as a mathematical function has rectangular asymptotes: The horizontal asymptote (x →  ± ∞) E’MAX is given by Eq. (8) and the vertical asymptote (y →  ± ∞) is equal to -KA/(τ + 1) (Supplementary Information Eq. (S2) and Eq. (S3)). From Eq. (7), the EC50 value is given by Eq. (9).

$$E^{\prime}_{MAX} = E_{MAX} \frac{{\tau^{n} }}{{1 + \tau^{n} }}$$
(8)
$${\text{EC}}_{50} = {\text{K}}_{{\text{A}}} { }\frac{1}{{\sqrt[{\text{n}}]{{2 + {\uptau }^{{\text{n}}} }}{ } - { }1}}$$
(9)

Evidently, the introduced slope factor n affects both the observed maximal response E’MAX and the half-efficient concentration of agonist EC50 ( Fig. 1A and C). The influence of the slope factor on E’MAX is bidirectional (Supplementary Information Table S1, Figure S2). For operational efficacies τ > 1, an increase in the value of the slope factor increases E’MAX. (Figs. 1A and 2A blue lines). For operational efficacies τ < 1, an increase in slope factor decreases E’MAX (Figs. 1C and 2A yellow lines). The effect of the slope factor on E’MAX is the most apparent for low values of operational efficacy τ, making the estimation of model parameters of weak partial agonists impractical. Imagine strong agonist τ = 10 and weak agonist τ = 0.1. For n = 1: Strong agonist E’MAX is 90% and weak agonist E’MAX is 10% of system EMAX. For n = 2: Strong agonist E’MAX is 99% (one-tenth more) and weak agonist E’MAX is just 1% (ten times less).

Figure 1
figure 1

Theoretical concentration–response curves. Theoretical curves of concentration responses under equilibrium according to Eq. (7), left and Eq. (11), right. Simulation parameters: EMAX = 1; τ = 3 (top) or τ = 0.3 (bottom); KA = 10−6 M. Values of slope factors are listed in the legend.

Full size image

An increase in the value of the slope factor increases the EC50 value (Fig. 2B). Again, the effect of the slope factor on the EC50 value is more eminent at low values of operational efficacy τ (yellow lines). Paradoxically, any combination of operational efficacy τ and slope factor fulfilling the inequality in Fig. 2C (blue area) results in EC50 values greater than KA (e.g., Fig. 1C, yellow lines). For example, EC50 > KA applies if τ = 0.5 and n > 1.6, or if τ = 1 and n > 1.6, or when τ = 1.5 and n > 2.15. The upper asymptote of inequality is 2. Thus, the possibility of EC50 > KA applies to τ < 2 in these cases, due to the high slope factor n the response function (Eq. (2)) has a “lag” before the steep growth where the binding function (Eq. (1)) grows faster. Therefore, in this range EC50 > KA.

Figure 2
figure 2

Analysis of Black & Leff equation (Eq. (7)). (A) Dependency of observed E’MAX to system EMAX ratio (ordinate) on slope factor n (abscissa) and operational efficacy τ (legend). (B) Dependency of EC50 to KA ratio (ordinate) on slope factor n (abscissa) and operational efficacy τ (legend). (C) Inequality plot of slope factor n (abscissa) and operational efficacy τ (ordinate) yielding half-efficient concentration EC50 greater than equilibrium dissociation constant KA.

Full size image

The operational efficacy τ may be also considered as a measure of “receptor reserve”. In a system with a relatively small capacity of a functional response output, the strong agonist reaches its maximal response before reaching full receptor occupancy. Thus, the agonist EC50 value is lower than its affinity for the receptor. The smaller the occupancy fraction needed for the full response to a given agonist the greater is difference between agonist EC50 and KA values. According to OMA (Eq. (2)), the relation between EC50 and KA is described by Eq. (6). The greater value of operational efficacy τ, the smaller EC50 value and the greater the difference from KA. Thus, the value of operational efficacy τ is a measure of the receptor reserve of a given agonist in a given system. In a system with a large capacity of functional output, agonists do not have a receptor reserve and must reach full receptor occupancy to elicit a full signal. In such a system, even strong agonists produce small effects attributed to partial agonists. Thus, the parameter τ is specific to a combination of ligand and system.

Nevertheless, for agonists that elicit at least some response in a given system, the parameter τ must be greater than 0. Then according to Eq. (6) of the operational model of agonism, the EC50 value must be smaller than the KA value. In principle, the EC50 value greater than the KA can be achieved only by some parallel mechanism that increases the apparent K’A provided that a ratio of K’A to KA is greater than EC50 to KA. For example, such a mechanism may be negative allosteric modulation of agonist binding or non-competitive inhibition of functional response.

Limitations of the non-hyperbolic OMA

Besides all limitations of the hyperbolic OMA, the non-hyperbolic version of OMA has additional drawbacks. The most important is the lack of mechanistic background for factor n. Exponentiation of agonist concentration [A] to power factor n results in non-hyperbolic functional-response curves. Importantly, as shown above, exponentiation of operational efficacy τ to power factor n breaks the logical relationship between observed maximal response E’MAX and operational efficacy τ. That, as it will be shown later, impedes the correct estimation of τ. Further, exponentiation of τ may result in KA values smaller than EC50 (Fig. 2C, blue area).

OMA with Hill coefficient

Definition of OMA with Hill coefficient

The Hill coefficient may serve as an alternative slope factor in the OMA. Hill equation incorporates the Hill coefficient as a slope factor to rectangular hyperbola12. The major advantage of the Hill coefficient as a slope factor is that it allows for a change in the eccentricity (vertices) of the hyperbola-like curves without changing the centre (EC50) and asymptotes (E’MAX) (Supplementary Information Figure S3). The Hill equation published in 191011 can be formulated as15 Eq. (10).

$$E = \frac{{\left[ A \right]^{{n_{H} }} E_{MAX}^{\prime } }}{{\left[ A \right]^{{n_{H} }} + EC_{50}^{{n_{H} }} }}$$
(10)

where nH is the Hill coefficient. Substitution of E’MAX by Eq. (5) and EC50 by Eq. (6) gives Eq. (11). For derivation from first principles see for example Roche et al.16

$$E = \frac{{\left[ A \right]^{{n_{H} }} \frac{\tau }{\tau + 1}E_{MAX} }}{{\left[ A \right]^{{n_{H} }} + \left( {\frac{{K_{A} }}{\tau + 1}} \right)^{{n_{H} }} }}$$
(11)

As expected, the Hill coefficient does not influence the maximal observed response E’MAX or half-efficient concentration of agonist EC50 (Fig. 1B,D). Eq. (11) was suggested as suitable for the analysis of functional responses displaying symmetrical response curves16.

Implications of OMA with Hill coefficient

Analysis of the OMA with slope factor by Black et al. (Eq. (7)) has shown that the slope factor n has a bidirectional effect on the relationship between the parameters E’MAX and τ. and that the slope factor n affects the relationship between the parameters EC50 and KA. In contrast, in Eq. (11) neither the value of E’MAX nor EC50 is affected by the Hill coefficient (Fig. 1B,D). The parameters E’MAX and EC’50 can be readily obtained by fitting Eq. (10) to the single concentration–response data.

Limitations of OMA with Hill coefficient

The major criticism of the Hill equation is its parsimonious character. It is relatively simple and its parameters are easy to estimate. However, as a model, it is just an approximation. In an experiment, the slope of the concentration–response curve different from unity may be a result of the parallel signalling mechanism providing feedback or allosteric cooperativity. In the case of positive cooperativity, it results in steep concentration–response curves. In the case of negative cooperativity, it results in flat concentration–response curves.

OMA of allosteric systems

The simplest scenario leading to variation in the slope of concentration–response curves is allosteric interaction between agonist and allosteric modulator or two molecules of agonist, e.g., in a ligand-gated channel or dimeric GPCR17,18. Positive cooperativity among agonist molecules results in steep functional-response curves and negative cooperativity results in flat functional-response curves. As the mode of cooperativity is the property of a given agonist, the slope of the functional-response curve may vary among agonists. Slope factor n in Eq. (7) is deemed the property of the system and thus the same for all agonists. In the attempt to keep slope factor n constant in allosteric systems, the second slope factor m was introduced to Eq. (7) resulting in19 Eq. (12):

$$E = \frac{{\left[ A \right]^{nm} \tau^{n} E_{MAX} }}{{\left[ A \right]^{nm} \tau^{n} + \left( {\left[ A \right]^{m} + K_{A}^{m} } \right)^{n} }}$$
(12)

where factor n links agonist concentration to the slope of the functional-response curve and factor m links agonist concentration to the slope of the receptor-binding curve. As the slope factors are interdependent one of them must be predetermined. It is only possible to predetermine the slope of the binding curve m. When fitting Eq. (12) to the data, functional-response slope factor n should become the same for all agonists. Observed maximal response E’MAX is still given by Eq. (8) but half-efficient concentration EC50 is given by Eq. (13).

$$EC_{50} = K_{A} \sqrt[m]{{\frac{1}{{\sqrt[n]{{2 + \tau^{n} }} - 1}}}}$$
(13)

Similarly to Eq. (9), the effect of the binding slope factor m on the EC50 to KA ratio depends on operational efficacy τ and functional response slope factor n. The effect of binding slope factor m depends on the radicand value. The m-th root is greater than the radicand when m > 1 and radicand > 1 or m < 1 and radicand < 1. The m-th root is smaller than the radicand when m < 1 and radicand > 1 or m > 1 and radicand < 1. Moreover, functional responses of allosteric systems may be not only steep or flat but also biphasic including bell-shaped18. Neither Eq. (7) nor Eq. (12) can accommodate such shapes and equations adequate to the mode of action are needed.

OMA of non-competitive inhibition

As shown in Fig. 2, OMA with slope factor n allows for EC50 values higher than KA. The simplest mode of interaction that increases observed EC50 above KA is non-competitive auto-inhibition20,21,22,23. Under non-competitive auto-inhibition, RA non-competitively blocks functional response by binding to effector G (Fig. 3). This type of inhibition is characterized by RA binding to a spatially distinct (allosteric) site resulting in a decreased response of effector G. In non-competitive inhibition, RA binds to both sites independently, exerting neutral binding cooperativity (absence of allosteric interaction). Non-competitive auto-inhibition results in a concentration-dependent increase in EC50 and a decrease in E’MAX. Functional response is given by Eq. (14).

$$E = E_{MAX} \frac{{\left[ {RA} \right]}}{{K_{E} + \left[ {RA} \right]}} \frac{{\left[ {RA} \right]}}{{K_{I} + \left[ {RA} \right]}}$$
(14)
Figure 3
figure 3

Non-competitive auto-inhibition of functional response. Dots, functional response to an agonist (KI = 1, EMAX = 1, RT = 1, KA = 10−6 M) in the model of non-competitive autoinhibition according to Eq. (15). Values of operational efficacies τ are indicated in the legend. Full lines, left, Black & Leff equation (Eq. (7)), right, Hill equation (Eq. (10)) fitted to the data. Parameter estimates are in Table 1.

Full size image

For KI > 0, after the substitution of Eq. (1) for [RA], Eq. (14) becomes Eq. (15) (Supplementary information Eq. (S31)).

$$E = E_{MAX} \frac{\tau \left[ A \right]}{{K_{A} + \left[ A \right]\left( {1 + \tau } \right)}} \frac{\sigma \left[ A \right]}{{K_{A} + \left[ A \right]\left( {1 + \sigma } \right)}}$$
(15)

where σ = RT/KI. The apparent maximal response E’MAX is given by Eq. (16) (Supplementary Information Eq. (S32)).

$$E^{\prime}_{MAX} = E_{MAX} \frac{\tau }{1 + \tau } \frac{\sigma }{1 + \sigma }$$
(16)

Thus, apparent operational efficacy τ’ is given by Eq. (17).

$$\tau^{\prime } = \frac{\sigma \tau }{{\sigma + \tau + 1}}$$
(17)

The EC50 value for the model of non-competitive auto-inhibition is given by Eq. (18) (Supplementary information Eq. (S36)).

$$EC_{50} = K_{A} \frac{{\sqrt {\sigma^{2} + 6\sigma \tau + 8\sigma + \tau^{2} + 8\tau + 8} + \sigma + \tau + 2}}{{2\left( {1 + \sigma } \right)\left( {1 + \tau } \right)}}$$
(18)

Non-competitive auto-inhibition decreases apparent maximal response E’MAX, increases observed half-efficient concentration EC50 and results in steep curves (Fig. 3). The resulting concentration–response curve is asymmetric with a typical slope factor of about 1.2 regardless of values of KI and KE. Both models fit well. Fitting Hill equation gives correct estimates of apparent maximal response E’MAX and thus correct estimates of apparent operational efficacy τ’ (Table 1). In contrast, for KI ≥ 1, fitting the Black & Leff equation results in underestimated values of τ and pKA. Values of τn well approximate τ’ values.

Table 1 Results of fitting Black & Leff equations to the model of non-competitive auto-inhibition. Black & Leff (Eq. (7)) and Hill (Eq. (10)) equations were fitted to model data Eq. (15) with EMAX fixed to 1 or E’MAX confined <  = 1, respectively.
Full size table

Fitting Eq. (15) with fixed system EMAX to the model of functional response of non-competitive inhibition yields correct parameter estimates that are associated with the low level of uncertainty only when correct initial estimates of τ and σ are given (Supplementary Information Figure S7 and S8). In the case of KI = 5 (Supplementary Information Figure S9), estimates of operational efficacy τ and inhibition factor σ are swapped pointing to the symmetry of Eq. (15). This symmetry makes calculation of τ and σ impossible as any τ and σ combination resulting in an appropriate apparent efficacy τ’ (Eq. (17)) corresponds well to a given functional-response data (Supplementary Information Figure S8). Fitting the Black & Leff equation (Eq. (7)) yields wrong estimates of KA and underestimated values of τ. Importantly, the extent of underestimation varies. The τ of 0.2 was underestimated by 17%. The τ of 5 was underestimated sixfold. In contrast, the calculation of apparent operational efficacies from the fitting of the Hill equation (Eq. (10)) is very close.

Signalling feedback

Signalling feedback occurs when outputs of a system are routed back as inputs24. Negative feedback is a very common auto-regulatory (auto-inhibitory) mechanism in nature, for example at G-protein coupled receptors25. Positive feedback also occurs in biology to propagate signals that would be otherwise dampened by other mechanisms, e.g. neuronal action potential. In the receptor-effector systems, an increase in output signal [RAG] proportionally either decreases (negative feedback) or increases (positive feedback) input [RA]26. Functional response in employing feedback mechanisms is then given by Eq. (19). For derivation see Supplementary Information Eqs. (S25) to (S37).

$$E = E_{MAX} \frac{{\tau \left[ A \right]\left( {\tau \left[ A \right] + \left[ A \right] + K_{A} } \right)}}{{\left[ A \right]^{2} \left( {\delta \tau + \tau^{2} + \tau + 1} \right) + \left[ A \right]K_{A} \left( {\delta \tau + \tau + 2} \right) + K_{A}^{2} }}$$
(19)

where δ is a feedback factor. Values greater than 1 denote negative feedback. Values smaller than 1 denote positive feedback. The maximal response E’MAX to an agonist with operational efficacy τ is related to the maximal response of the system EMAX according to Eq. (20).

$$E{^\prime}_{MAX} = E_{MAX} \frac{{\tau^{2} + \tau }}{{\delta \tau + \tau^{2} + \tau + 1}}$$
(20)

Negative feedback decreases the observed maximal response to an agonist E’MAX while positive feedback increases it. And apparent operational efficacy τ’ of a given agonist can be calculated according to Eq. (21).

$$\tau^{\prime } = \frac{{\tau^{2} + \tau }}{\delta \tau + 1}$$
(21)

Furthermore, negative feedback decreases the KA to EC50 ratio while positive feedback increases it (Supplementary Information Eq. (S54)). For derivations see Supplementary Information Eq. (S37) through Eq. (S54). In general, negative feedback results in flat functional-response curves as with the increase in signal output proportionally more of an agonist is needed for the same increase of the signal. Conversely, positive feedback results in steep functional-response curves as signal output proportionally increases signal input (Supplementary Information Figure S9).

Importantly, in a system with constant negative feedback, like in Fig. 4, the steepness of the curve depends on operational efficacy τ. Functional-response curves to agonists with high operational efficacy are flatter than the ones of agonists with low operational efficacy (Table 2).

Figure 4
figure 4

Functional response with signal feedback. Dots, functional response to an agonist (EMAX = 1, RT = 1, KA = 10−6 M) in the model of the system with constant negative feedback (δ = 5) according to Eq. (19). Values of operational efficacies τ are indicated in the legend. Full lines, left, Black & Leff equation (Eq. (7)), right, Hill equation (Eq. (10)) fitted to the data. Parameter estimates are in Table 2.

Full size image
Table 2 Results of fitting Black & Leff equations to the model of the system employing signalling feedback. Black & Leff (Eq. (7)) and Hill (Eq. (10)) equations were fitted to model data with EMAX fixed to 1 or E’MAX confined <  = 1, respectively.
Full size table

Fitting Eq. (19) with fixed system EMAX to the model employing feedback mechanisms yields correct parameter estimates that are associated with the low level of uncertainty when correct initial estimates of τ and δ are given (Supplementary Information Figure S9 and S10). Fitting the Black & Leff equation (Eq. (7)) yields wrong estimates of KA and underestimated values of τ. Again, the extent of underestimation varies. The τ of 0.2 was underestimated by 67%. The τ of 5 was underestimated fourfold. However, values of τn well approximate τ’ values. In contrast, calculations of apparent operational efficacies from the fitting of the Hill equation (Eq. (10)) are very close.

Systems with a similar expression of receptor and effector ([RT]≈[GT])

The Eq. (2) and consequently Eq. (4) are valid only when either [RA] or [G] is constant. That requires [RT] >  > [GT] or [RT] <  < [GT]. Usually, systems exhibit receptor reserve, indicating [RT] >  > [GT]. In systems with low receptor expression, [RT] and [GT] may be similar. In a such system [RAG] as a function of [RA] is given by Eq. (22) (Supplementary Information Eq. (S62)).

$$\left[ {RAG} \right] = \frac{1}{2} \left( {K_{E} + \left[ {RA} \right] + \left[ {G_{T} } \right] - \sqrt {K_{E}^{2} + 2K_{E} \left( {\left[ {RA} \right] + \left[ {G_{T} } \right]} \right) + (\left[ {RA} \right] - \left[ {G_{T} } \right])^{2} } } \right)$$
(22)

where [RA] denotes the concentration of all receptor agonist complexes, free RA plus RA in complex with G (RAG) (Supplementary Information Eq. (S56)). [RAG] as a function of [A] (Supplementary Information Eq. (S64)) has only approximate solutions. Therefore, [RA] values as a function of [A] were calculated as [RA] according to Eq. (1) and used in Eq. (22) to model functional responses of the system with low receptor-expression level (Fig. 5). The resulting curves are steep and asymmetric. In the low receptor-expression system, the response reaches the maximum at lower agonist concentration due to receptor depletion that cuts off the signal early which results in an apparent increased steepness and, in turn, in a curve asymmetry. The greater the operational efficacy is, the steeper and more asymmetric functional-response curves are (Table 3). However, the observed operational efficacy τ’ is equal to the modelled operational efficacy.

Figure 5
figure 5

Functional response of the system with a similar expression of receptor and effector ([RT]≈[GT]). Dots, functional-response data modelled in two steps. First, binding was calculated according to Eq. (1). Then resulting [RA] was used in Eq. (22). ET = 10−6 M, RT = 10−5 M, KA = 10−7 M. Values of operational efficacies τ are indicated in the legend. Full lines, left, Black & Leff equation (Eq. (7)), right, Hill equation (Eq. (10)) fitted to the data. Parameter estimates are in Table 3.

Full size image
Table 3 Results of fitting Black & Leff equations to the model of the system with a similar expression of receptor and effector ([RT]≈[GT]). Black & Leff (Eq. (7)) and Hill (Eq. (10)) equations were fitted to model data with EMAX fixed to 1 or E’MAX confined <  = 1, respectively.
Full size table

Fitting the Black & Leff equation (Eq. (7)) to the model system with a low receptor-expression level yields wrong estimates of KA for extremely high operational efficacy (τ = 1000). Values of τ are far off for all modelled efficacies. While τ of 1000 is overestimated more than sixfold, the other τ values are underestimated up to threefold. Values of τn approximate model τ values of 33.3, 10 and 3.33. However, the τn value is overestimated 20,000 and 2.5 times for model values of 1000 and 100, respectively. Except for the τ value of 1000, operational efficacies calculated according to the Hill equation (Eq. (10)) are less than 10% off the model values.

The case study

How fitting the Black & Leff equation to experimental data can affect estimates of the operational efficacy and subsequent analysis is demonstrated in the following example of measurement of the GTPγS binding as a functional response of M2 receptor to muscarinic agonists carbachol, iperoxo, N-desmethyl clozapine (NDMC) and oxotremorine at five subtypes of inhibitory G-proteins: Gi1, Gi2, Gi3, GoA and GoB expressed in Sf9 cells. According to the saturation of [3H]NMS binding, Sf9 cell membranes contained about 7 pmol of binding sites per mg of membrane protein. Co-expression of individual subtypes of inhibitory G-proteins did not affect the expression level of M2 receptors. According to competition with [3H]NMS binding, tested agonists had the same affinity for M2 receptors in all co-expression systems (Supplementary Information Table S2).

When analysing the functional response, first, system EMAX values were determined using the procedure described earlier6. They ranged from 84 to 88% of the maximum binding capacity of G-proteins. After subtraction of basal binding, responses were expressed as a fraction of EMAX. The system EMAX was set equal to 1.

Signalling profiles varied among subtypes of G-proteins. For example, carbachol and oxotremorine reached similar maximal responses E’MAX at Gi1, GoA and GoB. At Gi2, the E’MAX of oxotremorine was substantially greater than the E’MAX of carbachol. In contrast, at Gi3, the E’MAX of oxotremorine was substantially lower than the E’MAX of carbachol (Fig. 6). Also, the steepness of functional responses to agonists varied among agonists as well as among subtypes of G-proteins. The functional response to carbachol was normal (nH ≈ 1), except for a flat response (nH = 0.68) at GoB (Supplementary information Table S2). The functional response to NDMC was steep (nH > 1.2) at Gi1 and Gi3, normal at Gi2 and flat (nH < 0.85) at GoA and GoB. As the exact reason for variation in the steepness of functional-response curves is unknown, the Hill equation (Eq. (11)) was fitted to the experimental data for comparison with the fitting of the Black & Leff equation (Eq. (7)). Results of fitting are summarized in Supplementary Information Table S2).

Figure 6
figure 6

Functional response of muscarinic M2 receptors. GTPγS binding to Gi1 (upper left), Gi2 (upper middle), Gi3 (upper right) GoA (lower left) and GoB (lower right) G-proteins upon stimulation of M2 muscarinic receptors by carbachol (red), iperoxo (green), NDMC (blue) or oxotremorine (yellow) at the concentration indicated on the x-axis is expressed as the fraction of maximal response of the system (EMAX). Data are means ± SD from 3 to 5 independent experiments performed in quadruplicates. Solid curves, Black & Leff equation (Eq. (7)). Dotted curves, Hill equation (Eq, (11)). Parameter estimates are in Supplementary Information Table S2.

Full size image

The analysis of estimates of equilibrium dissociation constants KA obtained from fitting Black & Leff (Eq. (7)) and Hill (Eq. (11)) equations shows that in comparison to Hill fits, estimates from Black & Leff fits are associated with greater variability (among individual fits) and uncertainty (individual fits) (Fig. 7). In the majority of cases, KA estimates according to the Black & Leff equation were substantially lower than according to the Hill equation. Except for NDMC, KA estimates for both models are higher than KI from competition with [3H]NMS. Notably, agonist KI values are the same for all subtypes of G-protein α-subunit while KA values vary among them (Supplementary Information Table S2 and S3).

Figure 7
figure 7

Analysis of pKA estimates. Estimates of equilibrium dissociation constant KA from fitting Black & Leff (Eq. (7)) (blue) or Hill (Eq. (11)) (yellow) equation are expressed as negative logarithms (pKA). Data are means ± SD from fits of 5 independent experiments performed in quadruplicates. *, different from Black & Leff (p < 0.05) according to ANOVA and Tukey HSD post-test.

Full size image

Similarly to estimates of KA, in comparison to Hill fits, estimates of operational efficacies τ from Black & Leff fits are associated with greater variability and uncertainty (Fig. 8). Overall, τ values tend to be underestimated for flat curves and overestimated for steep curves. Following analysis of the Black & Leff equation, this phenomenon is apparent for low operational efficacies, especially for τ < 1. In the extreme case of steep response to NDMC (n = 1.5) and flat response to oxotremorine (n = 0.7) at Gi3, Black & Leff estimates of τ are the same for both agonists despite oxotremorine elicited E’MAX three-times higher than NDMC. In the less extreme case, at GoB G-protein, oxotremorine was about 60% more efficacious than carbachol but estimated τ values by Black & Leff were the same for both agonists although functional-response curves of both agonists are flat (nH = 0.6 and 0.5, respectively). Thus, even a small change in the steepness of the functional-response curve has a profound effect on the estimation of the τ value.

Figure 8
figure 8

Analysis of τ estimates. Estimates of operational efficacy τ from fitting Black & Leff (Eq. (7)) (blue) or Hill (Eq. (11)) (yellow) equation are expressed as decadic logarithms (log τ). Data are means ± SD from fits of 5 independent experiments performed in quadruplicates. *, different from Black & Leff (p < 0.05) according to ANOVA and Tukey HSD post-test.

Full size image

Discussion

The operational model of agonism (OMA)2 is widely used in the evaluation of agonism. The OMA characterizes a functional response to an agonist by the equilibrium dissociation constant of the agonist (KA), the maximal possible response of the system (EMAX) and the operational efficacy of the agonist (τ) (Eq. (4)). To fit non-hyperbolic functional responses slope factor n was introduced to the OMA (Eq. (7))10. Analysis of the Black & Leff equation (Eq. (7)) has shown that the slope factor n has a bidirectional effect on the relationship between the parameters E’MAX and τ (Fig. 1A versus C) and also affects the relationship between the parameters EC50 and KA.

Fitting Black & Leff equation (Eq. (7)) to the models of non-competitive auto-inhibition (Fig. 3, Table 1), signalling feedback (Fig. 4, Table 2) and system with a similar expression of receptor and effector ([RT]≈[GT]) (Fig. 5, Table 3) resulted in inaccurate estimates of τ and KA values. In the presented examples, the degree of over- or under-estimation of τ is not the same for all its values but depends on the value of τ, distorting relations among estimates of τ vales. In contrast, fitting the Hill equation (Eq. (10)) to the model data gave more or less accurate estimates of apparent operational efficacies τ’ from which operational efficacies can be calculated (according to Eq. (17), and Eq. (21)), provided that mechanism of action is known.

Biased agonists stabilize specific conformations of the receptor leading to non-uniform modulation of individual signalling pathways27. To measure an agonist bias, the parameters τ and KA must be determined and log(τ/KA) values of tested and reference agonists compared at two signalling pathways3. It is evident from the analysis of the Black & Leff equation, that as far as the EC50 value is dependent on parameters n and τ (Fig. 1A,C), log(τ/KA) values cannot be compared to judge possible signalling bias unless the parameter n is equal to 1 for both tested and reference agonist.

Despite the dire effects of slope factor n, the Black & Leff equation is widely accepted4,28,29,30,31,32,33,34. It even entered textbooks5. Very little concern on factor n has risen. For example, Kenakin et al.3 analysed in detail the effects of slope factor n on EC50 and τ to KA ratio but did not deal with the bi-directional effect of n on τ nor proposed an alternative approach to avoid potential pitfalls. To force a proper shape on functional-response curves whilst keeping slope factor n constant for all ligands, Gregory et al.19 introduced the second slope factor into their OMA and operational model of allosterically modulated agonism (OMAMA) analysis making equations even more complex. So far, the greatest criticism of OMA was voiced by Roche et al.16, noting that to accommodate the shape of theoretical curves Black & Leff equation tends to overestimate equilibrium dissociation constant KA and operational efficacy τ and thus be misleading. They advocated for different expressions of operational models including OMA modified by the Hill coefficient in the case of symmetric concentration–response curves. Theoretical models show that in cases of non-competitive auto-inhibition and signal feedback where the slope of functional response does not vary, Black & Leff equation can be used to estimate apparent operational efficacy τ’ (Tables 1 and 2). However, estimates of operational efficacy do not approximate τ of the models.

The fitting of the Hill equation (Eq. (10)) to the functional response is straightforward and easier than fitting the Black & Leff equation. As shown in Figs. 3, 4, 5, 6, the Hill equation fits well with various functional-response curves, often better than the Black & Leff equation. In contrast to the Black & Leff equation, the Hill equation gives correct estimates of maximal response to agonist E’MAX and its half-efficient concentration EC50 as documented in Tables. 1, 2, 3 and Figs. 7 and 8. In the case of the Hill equation, neither the value of E’MAX nor the value of EC50 is affected by the Hill coefficient (Fig. 1B,D). Therefore, biased signalling may be inferred from the comparison of the ratio of intrinsic activity (E’MAX/EC50) of tested agonist to the intrinsic activity of reference agonist at two signalling pathways as in the case of the Hill equation the E’MAX/EC50 ratio is equivalent to τ/KA ratio35,36. Further, if needed, apparent operational efficacy τ’ can be calculated from known E’MAX values and known maximal response of the system EMAX. From the relationship between τ’ and EC506,17, the mechanism of functional response can be inferred by comparison to explicit models.

The case study of the functional response of individual subtypes of inhibitory G-proteins to activation of the M2 muscarinic receptor (Fig. 6) demonstrated the pitfalls of exponentiation of operational efficacy τ: At Gi3 G-protein, oxotremorine reached 3-times higher maximal response E’MAX than NDMC. However, estimates of operational efficacy τ by the Black & Leff equation were the same for both agonists (Fig. 8), rendering it unsuitable for such data. An analysis of the Black & Leff equation implies that the source of the discrepancy lies in the profound effects of slope factor n on the estimation of τ value. The τ estimates by the Hill equation reflected differences in E’MAX values, making it suitable for such a scenario. The root of the problem is variation in the slope of the functional response among tested agonists, indicating that the slope is not the property of the system required by the Black & Leff model. Rather agonists exert various degrees of cooperativity. As muscarinic receptors possess only one orthosteric binding site the data indicate that muscarinic receptors may function as oligomers18. Although cooperativity does not automatically mean oligomerization37, several lines of evidence indicate that muscarinic receptors may indeed oligomerize38,39,40,41.

Conclusions

Analysis of the Black & Leff equation has shown that (i) The slope factor n has a bidirectional effect on the relationship between the parameters E’MAX and τ. (ii) The slope factor n affects the relationship between the parameters EC50 and KA. Fitting the Black & Leff equation gives wrong estimates of τ and KA values when slope factor n varies among concentration–response curves, limiting the use of the Black & Leff equation to evaluate concentration–response curves with the same slope. Analysis of the Hill equation has shown that the Hill coefficient does not affect the relationship between the parameters E’MAX and τ nor between the parameters EC50 and KA. Fitting the Black & Leff equation to the experimental data demonstrated the drawbacks of exponentiation operational efficacy τ. In contrast, fitting the Hill equation to the experimental data gave more realistic estimates of KA and τ. Black & Leff equation may be safely used only for systems where the slope of functional response does not vary substantially to estimate apparent operational efficacy.

Methods

Models and equations

Models and equations were derived from scratch as described in Supplementary Information. For modelling the theoretical curves and fitting curves to the experimental data the Python code employing numpy, scipy and matplotlib libraries was written.

Preparation of cells and membranes

Spodoptera frugiperda cells (Sf9) (Gibco) were maintained as a suspension culture in serum-free insect cell growth medium SF900 III (Gibco) in a plastic Erlenmeyer flask in a shaking incubator at 27 °C and 135 rpm in a non-humidified environment. The cultures were maintained at a density of 1–4 × 106 cells/ml. The density of the cells was determined with a haemocytometer, and viability was assessed by the exclusion of 0.2% trypan blue (Sigma-Aldrich). Human M2 receptors and α-subunits of (Gi1, Gi2, Gi3, GoA or GoB) G-proteins were expressed in recombinant baculoviruses, which were constructed and generated according to Bac-to-Bac® Baculovirus Expression System manual (Life Technologies)42.

One hundred ml of Sf9 cell suspension at a density of 2 × 106 cells/ml were co-infected with baculoviruses encoding the M2 receptor and Go or Gi α-subunit at a multiplicity of infection MOI = 0.1: 0.1. All infections were allowed to proceed for 69 h. Infected cells were harvested by centrifugation at 500 × g for 5 min and frozen at − 80 °C.

The pellets of harvested cells were suspended in the ice-cold homogenization medium (100 mM NaCl, 20 mM Na-HEPES, 10 mM EDTA, pH = 7.4) and homogenized on ice by two 30 s strokes using a Polytron homogenizer (Ultra-Turrax; Janke & Kunkel GmbH & Co. KG, IKA-Labortechnik, Staufen, Germany) with a 30-s pause between strokes. Cell homogenates were centrifuged for 5 min at 1000 × g. The supernatant was collected and centrifuged for 30 min at 30,000 × g. Pellets were suspended in the washing medium (100 mM NaCl, 10 mM MgCl2, 20 mM Na-HEPES, pH = 7.4), left for 30 min at 4 °C, and then centrifuged again for 30 min at 30,000 × g. The resulting membrane pellets were kept at − 80 °C until assayed.

[3H]NMS binding

Membranes (approximately 10 μg of membrane proteins per sample) from Sf9 cells were incubated in 96-well plates for 3 (saturation) or 5 h (competition) at 25 °C in the incubation medium (100 mM NaCl, 20 mM Na-HEPES,10 mM MgCl2, pH = 7.4). The incubation volume for competition and saturation experiments with [3H]NMS was 400 μl or 800 μl, respectively.

In saturation experiments, eight concentrations of the [3H]NMS ranging from 94 to 1000 pM were used. Agonist binding was determined in competition experiments with 1 nM [3H]NMS. Nonspecific binding was determined in the presence of 10 μM unlabelled atropine. Incubation was terminated by filtration through Whatman GF/C glass fibre filters (Whatman) using a Brandel harvester (Brandel, USA). Filters were dried in a microwave oven (3 min, 800 W), and then solid scintillator Meltilex A was melted on filters (105 °C, 90 s) using a hot plate. The filters were cooled and counted in a Wallac Microbeta scintillation counter (Wallac, Finland).

GTPγS binding

Agonist-stimulated [35S]GTPγS binding was performed as currently described43. Briefly, it was performed in 96-well plates at 30 °C in the incubation medium described above that was supplemented with freshly prepared dithiothreitol at a final concentration of 1 mM. Suspension of membranes of Sf9 cells expressing M2 + G-protein α-subunit were preincubated with GDP and agonists for 15 min, then [35S]GTPγS was added for an additional 20 min. The final concentration of GDP and [35S]GTPγS was 20 µM and 500 pM, respectively. The maximum binding capacity of G-proteins was determined in the absence of GDP. Nonspecific binding was determined in the presence of 1 µM non-labelled GTPγS. Incubations were terminated by filtration through GF/C filtration plates (Millipore) using a Brandel cell harvester (Perkin Elmer, USA). Plates were dried in a microwave oven at 800W for 3 min and then 40µl of ROTISZINT® Eco Plus (ROTH) was added. The plates were counted in the Wallac Microbeta scintillation counter.

Analysis of experimental data

Experimental data were analysed using the two-step procedure described earlier6,17. First, the maximum system response EMAX was determined. After subtracting the value of the basal signal, functional responses were expressed as the fraction of corresponding EMAX. Then, the Hill equation (Eq. (10)) and the Black & Leff equation (Eq. (7)) were fitted to the data. In the case of the Hill equation, the maximum response to an agonist E’MAX was confined to less than 1. In the case of the Black & Leff equation, system EMAX was set equal to 1.