The influence of the way of regression on the results obtained by the receptorial responsiveness method (RRM), a procedure to estimate a change in the concentration of a pharmacological agonist near the receptor

The receptorial responsiveness method (RRM) enables the estimation of a change in concentration of an (even degradable) agonist, near its receptor, via curve fitting to (at least) two concentration-effect (E/c) curves of a stable agonist. One curve should be generated before this change, and the other afterwards, in the same system. It follows that RRM yields a surrogate parameter (“c_x”) as the concentration of the stable agonist being equieffective with the change in concentration of the other agonist. However, regression can be conducted several ways, which can affect the accuracy, precision and ease-of-use. This study utilized data of previous ex vivo investigations. Known concentrations of stable agonists were estimated with RRM by performing individual (local) or global fitting, this latter with one or two model(s), using a logarithmic (logc_x) or a nonlogarithmic (c_x) parameter (the latter in a complex or in a simplified equation), with ordinary least-squares or robust regression, and with an “all-at-once” or “pairwise” fitting manner. We found that the simplified model containing logc_x was superior to all alternative models. The most complicated individual regression was the most accurate, followed closely by the moderately complicated two-model global regression and then by the easy-to-perform one-model global regression. The two-model global fitting was the most precise, followed by the individual fitting (closely) and by the one-model global fitting (from afar). Pairwise fitting (two E/c curves at once) improved the estimation. Thus, the two-model global fitting, performed pairwise, and the individual fitting are recommended for RRM, using the simplified model containing logc_x.

1 Introduction

Measuring the concentration of small molecules in biological samples is generally considered to be an issue for analytical chemistry that can now be performed to a high standard (Valcárcel, 2000; Belanger et al., 2020). However, concentration determination of small molecules with a short half-life and/or strong compartmentalization tendency in living, functioning, moreover moving tissues is still a challenge (Karsai et al., 2006; Ramakers et al., 2008; Zhang et al., 2015; Markiewicz et al., 2022). If certain requirements are met, the receptorial responsiveness method (RRM) can help to solve (or bridge) this problem. RRM can quantify an acute (and occasionally even chronic) increase in the concentration of a pharmacological agonist in the microenvironment of its receptors (Gesztelyi et al., 2004; Karsai et al., 2006). (Using this method “in reverse” (i.e., by reversing the temporality), it is of course also possible to determine a decrease in concentration.)

The principle of RRM may seem to be complicated at first glance but is in fact straightforward (the main challenge to understand it stems from its interdisciplinary nature). First, it is worth considering a thought experiment. Let us administer the same dose of an agonist twice to a biological system, and then let us compute the effect of both doses. The correct way is to define the same initial state for both effects (namely, the state before the administration of the first dose) and to take also the first dose into account when calculating the effect of the second one (because, due to its magnitude relative to the second dose, the first dose is not negligible). However, if the effect of the second dose is calculated without taking the first dose into account, and the condition after the first dose is taken as the initial state, the effect computed for the second dose will be smaller than that for the first dose (although the two doses equal).

The explanation for this smaller effect is that the agonist molecules administered with the first dose, by binding to a fraction of the receptors and then by eliciting an effect on them, consume part of the response capacity of the system. For this reason, the agonist molecules administered with the second dose will act in a system with submaximal response capacity. If, for example, the dose in question is high enough to exert a (practically) maximal effect after its first administration, no further effect can develop, thus, when administering it for the second time, the system will be unable to respond. Thus, the decrease in the responsiveness (observed this way) depends on the magnitude of the ignored dose, more precisely on the agonist concentration near the relevant receptors. This phenomenon can be exploited to quantify an unknown concentration of an agonist developed in the tissue compartment of the binding site of the specific receptors (Gesztelyi et al., 2004; Karsai et al., 2006; 2007; Grenczer et al., 2010a).

It should be emphasized that this phenomenon, although superficially reminiscent of it, is not based on receptor desensitization (where the effect is correctly computed but smaller than expected). In turn, the receptor desensitization interferes with the above-mentioned phenomenon, thus this latter can only be used for concentration estimation if the receptor in question desensitizes slowly relative to the time window of the determination (in our experiments: about 10–40 min, depending on the circumstances). The A₁ adenosine receptor belongs to the especially slowly desensitizing receptors (Willems et al., 2006; Klaasse et al., 2008; Mundell and Kelly, 2011).

Of course, the effect of an agonist depends not only on its concentration in the vicinity of its receptors and on the genuine sensitivity of these receptors but also on numerous other factors (e.g., receptor density, receptor reserve, receptor isoforms, receptor trafficking, postreceptorial signaling). Moreover, the tissue concentration of the given agonist is also influenced by several factors (such as stability of the agonist including its resistance to tissue enzymes, carriers and possibly reactive endogenous compounds) (for a review, see: Kenakin and Pharmacology primer, 2022). However, when comparing two states of a system (in our case, to assess a decrease in the responsiveness of the system), the non-varying system characteristics can be cancelled so as not to disturb the determination of a single cause of the difference (which, in our case, is the absence or presence of a disregarded agonist concentration). This is a kind of application of the pharmacological null method (or null analysis), which is a useful procedure to negate cell- and tissue-dependent characteristics irrelevant for the given investigation. Nevertheless, the cost of this advantage is that measurements based on the null method require the permanence of the biological system throughout the investigation (Kenakin and Pharmacology primer, 2022).

The phenomenon described above can also be observed if the agonist being present first and being disregarded is different from that added subsequently, but they both evoke the same response (so the same response capacity is consumed by them). To do so, it is enough if they are agonists on the same receptor, or even if they influence the same postreceptorial signaling in the same direction (although binding to different receptors). So, the essence of the phenomenon is the bias in the measurement that is caused by the neglect of the first agonist. Thus, the first and disregarded agonist can be named as “biasing” or “distorting” agonist (Gesztelyi et al., 2004; Grenczer et al., 2010a).

It should be noted that the terms “bias” and “biasing” mentioned above are not related to the emerging concept of biased agonism. According to that concept, a biased (also called functionally selective) agonist is the one that does not activate every signaling pathway linked to the given receptor, but it does stimulate at least one of them. The phenomenon of biased agonism has extensively been observed for G protein-coupled receptors (Jarpe et al., 1998; Andresen, 2011; Maguire, 2016; Kenakin, 2024), including the A₁ adenosine receptor (Franco et al., 2021; McNeill et al., 2021; IJzerman et al., 2022). Thus, to avoid misunderstanding, the agonist first present and disregarded will be called as distorting (and not biasing) agonist, and its concentration will be referred to as distorting concentration throughout this paper. Again, the term “distortion” used herein means an inevitable modification of the measured data and the consequent transformation of the curves generated from these data, under well-circumscribed conditions (see above).

It is worth noting here that NECA, CPA and CHA, three high-efficacy A₁ adenosine receptor agonists responsible for data evaluated in the present investigation (for more detail, see legends of Figures 1–3), activate identical postreceptorial signaling pathways. However, their activation patterns slightly differ, indicating some level of functional selectivity for these agonists (Verzijl and Ijzerman, 2011). Despite their differences, NECA, CPA and CHA all have relatively long half-life in the living tissues (Pavan and IJzerman, 1998; Fredholm et al., 2001), ensuring the permanence of the experimental system that is needed for RRM (Gesztelyi et al., 2004; Karsai et al., 2006).

Figure 1

Figure 1. The implementation of RRM with individual regression combined with ordinary least-squares fitting, using two equations: the Eq. 8 that contained the common logarithm of the concentration to be determined (logc_x; left panels), and the Eq. 7 containing this concentration itself (c_x; right panels), as a parameter. The measurement was performed on “Intact”–“Distorted” concentration-effect (E/c) curve pairs obtained from a previous ex vivo study (Szabo et al., 2019a). The x-axis shows the common logarithm of the molar concentration of one of three synthetic A₁ adenosine receptor full agonists (CPA, NECA, CHA) that was administered for the construction of the E/c curves. The y-axis shows the effect (expressed as the percentage decrease in the initial contractile force of isolated, paced guinea pig left atria). The “Intact” E/c curves (filled symbols) were generated in a conventional way, whereas the “Distorted” E/c curves (open symbols) were constructed in the presence of a single, disregarded extra concentration of the agonist used for the E/c curve. The symbols indicate the responses to the agonists averaged within the groups (±SEM). The continuous lines show the best-fit curves of the fitted Eqs 7, 8 (encompassing c_x or logc_x, respectively). The thick dotted lines indicate the 95% confidence bands, while the thin dotted lines denote the 95% prediction bands. RRM: receptorial responsiveness method; CPA: N⁶-cyclopentyladenosine; NECA: 5′-(N-ethylcarboxamido) adenosine; CHA: N⁶-cyclohexyladenosine.

Figure 2

Figure 2. The application of RRM with one-model global regression combined with ordinary least-squares fitting, using two equations: the Eq. 8 containing the common logarithm of the concentration to be determined (logc_x; left panels), and the Eq. 7 including this concentration itself (c_x; right panels), as a parameter. This procedure was performed on “Intact”–“Distorted” concentration-effect (E/c) curve pairs obtained from an earlier ex vivo study (Szabo et al., 2019a). The x-axis shows the common logarithm of the molar concentration of one of three synthetic A₁ adenosine receptor full agonists (CPA, NECA, CHA) that was administered for generating the E/c curves. The y-axis indicates the effect (the percentage decrease in the initial contractile force of isolated, paced guinea pig left atria). The “Intact” E/c curves (filled symbols) were constructed conventionally, whereas the “Distorted” E/c curves (open symbols) were generated in the presence of a single, disregarded surplus concentration of the agonist used for the E/c curve. The symbols indicate the responses to the agonists averaged within the groups (±SEM). The continuous lines show the best-fit curves of the fitted Eqs 7, 8. The thick dotted lines indicate the 95% confidence bands, while the thin dotted lines denote the 95% prediction bands (if any). RRM: receptorial responsiveness method; CPA: N⁶-cyclopentyladenosine; NECA: 5′-(N-ethylcarboxamido)adenosine; CHA: N⁶-cyclohexyladenosine.

Figure 3

Figure 3. The implementation of RRM with two-model global regression combined with ordinary least-squares fitting, using two equations (as alternatives in addition to the Eq. 2 representing the Hill model): the Eq. 8 that contained the common logarithm of the concentration to be determined (logc_x; left panels), and the Eq. 7 containing this concentration itself (c_x; right panels), as a parameter. This procedure was performed on “Intact”–“Distorted” concentration-effect (E/c) curve pairs obtained from a previous ex vivo study (Szabo et al., 2019a). The x-axis shows the common logarithm of the molar concentration of one of three synthetic A₁ adenosine receptor full agonists (CPA, NECA, CHA) that was administered for constructing the E/c curves. The y-axis indicates the effect (the percentage decrease in the initial contractile force of isolated, paced guinea pig left atria). The “Intact” E/c curves (filled symbols) were generated in a conventional way, while the “Distorted” E/c curves (open symbols) were constructed in the presence of a single, disregarded extra concentration of the agonist used for the E/c curve. The symbols indicate the responses to the agonists averaged within the groups (±SEM). The continuous lines show the best-fit curves of the Eq 2 and the Eqs 7, 8. The Eq. 2 was fitted to the “Intact” E/c curves, whereas the Eq. 8 (left panels) or 7 (right panels) was fitted to the “Distorted” ones. The thick dotted lines indicate the 95% confidence bands, while the thin dotted lines represent the 95% prediction bands. RRM: receptorial responsiveness method; CPA: N⁶-cyclopentyladenosine; NECA: 5′-(N-ethylcarboxamido)adenosine; CHA: N⁶-cyclohexyladenosine

A complete concentration-effect (E/c) curve (ranging from ineffective to saturative concentrations) encompasses the information needed for this concentration estimation with high redundancy, and thereby offers the possibility of a more reliable estimation than what a single concentration - effect pair could provide. Consequently, the best way to assess a change (“distortion”) in a biological response is to analyze E/c curves. Two types of E/c curves are required: one holding information on the real (intact) relationship between agonist concentration and effect in the given system (serving as a benchmark), and another one informing about the extent of the distortion (and thereby, in our case, about the magnitude of the distorting concentration). These two E/c curves can be referred to as “intact” and “distorted” that contain “intact” and “distorted” effects, respectively. The intact effect and the corresponding distorted effect are assigned to the same agonist concentration, the only difference between them is the presence and neglect of a distorting concentration in the case of the distorted effect (Gesztelyi et al., 2004; Grenczer et al., 2010a). (It is worth anticipating that if more than one distorting concentration is to be estimated in the same system, then more than two E/c curves are to be assessed with RRM. This situation can technically be handled in two different ways, see further on).

The essence of RRM is to fit its model to a distorted E/c curve, utilizing the relevant pieces of information implied in the corresponding intact E/c curve. The model of RRM is a fusion of two equations: 1) the basic equation of RRM (see later) that describes the relationship between an intact effect and its distorted counterpart, and 2) a quantitative receptor function model (that is generally the Hill equation: Gesztelyi et al., 2012a). The information in the intact E/c curve can be made accessible two main ways. On one hand, the given receptor function model can be fitted to the intact E/c curve, and then the identical parameters in the RRM’s model can be constrained at the obtained values. After this “individualization”, the RRM’s model can simply be fitted to the corresponding distorted E/c curve (Gesztelyi et al., 2004; Grenczer et al., 2010a; 2010b). In this paper (similarly to the recent ones: Szabo et al., 2019a; Viczjan et al., 2021; 2022), this procedure will be referred to as individual regression. (Alternatively, this type of fitting can be called local regression, as opposed to global regression: Anjos et al., 2020; Spatial Data Science, 2024). On the other hand, the intact and distorted E/c curves can be fitted at once, via global regression. Global regression can also be performed in two ways: fitting the intact and distorted E/c curves to one model (which is the model of RRM: Szabo et al., 2019a), or fitting these E/c curves to two models (the selected receptor function model and the compatible RRM’s model). This latter maneuver should be done so that the receptor function model fits only the intact E/c curve and the RRM’s model fits exclusively the distorted E/c curve, with the common parameters shared between the two models (Curve Fitting Guide, 2024). Hereinafter, these two procedures will be referred to as one-model global regression and two-model global regression, respectively.

Finally, RRM provides a best-fit value as a measure for the distorting concentration. This best-fit value can be a concentration (c_x) or its common logarithm (logc_x), depending on which parameter (c_x or logc_x) is included in the RRM’s model used. So, the estimate yielded by RRM is the best-fit value, if it is c_x, or the antilog of the best-fit value, if this latter is logc_x. The nature of this estimate depends on the fact whether the distorting agonist and the agonist used for the E/c curves are the same or not. If so, the estimate directly applies to the distorting concentration. If not, then the estimate is the concentration of the agonist used for the E/c curves that is equieffective with the distorting concentration (Gesztelyi et al., 2004; Grenczer et al., 2010a). Although this latter estimate is just a surrogate, it has the advantage that a stable agonist can be used to determine the distorting concentration of an agonist that may be degradable and/or exhibit strong compartmentalization. Importantly, since the estimate provided by RRM is derived from E/c curves, it applies to the (average) concentration of the distorting agonist near the receptors involved in the production of the effect measured (Gesztelyi et al., 2004; Karsai et al., 2006; 2007; Viczjan et al., 2022).

In a previous study (Szabo et al., 2019a), estimates of RRM were found to be more accurate and precise when obtained from individual regression rather than one-model global fitting. This may be surprising in light of that global regression is thought to be more powerful and reliable than the individual one, especially in the case of limited data and/or data with much scatter (Knutson et al., 1983; Motulsky and Christopoulos, 2004; Herman and Lee, 2012). The worse than expected performance of the one-model global fitting was then attributed to the lack of accurate estimates for the intact E/c curves, because, in that study (Szabo et al., 2019a), the model of RRM contained logc_x as a parameter (instead of c_x). Hence, RRM could not yield zero, the correct estimate for the intact E/c curves, as the logarithm of zero is not defined.

Earlier computer programs calculated only symmetrical confidence intervals (CIs). Therefore, quantities typically following a log-normal distribution (e.g., concentrations) were recommended to be used as common logarithms in the equations to fit, in order to get correct CIs (Motulsky and Christopoulos, 2004). If asymmetry of CIs can be handled (not uncommon for modern statistical software systems), the use of concentrations as a logarithm in the regression models seems to be less important.

Consequently, in the present study, our first goal was to explore whether the model of RRM containing c_x is superior to that including logc_x (during one-model global regression). As an alternative approach to the issue of logc_x as a fitted parameter, we also investigated whether the two-model global regression, another way to avoid fitting logc_x to an intact E/c curve, can provide better results than the one-model global regression and perhaps than the individual fitting. Furthermore, we aimed to find out whether reducing the complexity of the model of RRM improves the results of RRM. In this context, the influence of dividing by c_x in the RRM’s model (to simplify its expression) had to be considered as well. In addition, regarding their impact on the outcome of RRM, the distribution of the scatter of E/c curve data and the fitting to E/c curve families (where more than one distorted E/c curve belongs to one intact E/c curve) were also addressed.

2 Materials and methods

2.1 Data analyzed

Investigations of the present study were carried out on cumulative, graded E/c curves, originally published in two earlier studies from our lab (Gesztelyi et al., 2004; Szabo et al., 2019a). The E/c curves were constructed with A₁ adenosine receptor agonists, by measuring the contractile force of isolated, paced guinea pig left atria. These E/c curves were E/c curve pairs or families (one intact E/c curve with one or more than one related distorted E/c curve(s), respectively). These E/c curves were reevaluated with RRM implemented by combining several ways of regression, in order to compare the usefulness of the different possibilities. The intact - distorted E/c curve pairs (obtained from Szabo et al., 2019a) were reexamined by combining the following options: 1) the RRM’s model containing logc_x vs. c_x; 2) the originally developed (Gesztelyi et al., 2004), more complex RRM’s model (containing c_x) vs. an RRM’s model simplified with a division by c_x (containing c_x or logc_x); 3) individual regression vs. one-model global regression vs. two-model global regression; and 4) ordinary least-squares (simply: ordinary) fitting vs. robust fitting. To reevaluate the E/c curve families (obtained from Gesztelyi et al., 2004), the following options were combined: 1) individual fitting vs. one-model global fitting vs. two-model global fitting; 2) ordinary regression vs. robust regression; and 3) fitting all related E/c curves at once (i.e., all-at-once manner) vs. fitting the intact E/c curve together with only one related distorted E/c curve at once (i.e., pairwise manner).

The intact E/c curves were generated with CPA, NECA or CHA, stable, synthetic A₁ adenosine receptor agonists with long half-life (for more detail, see legends of Figures 1–3), in the absence of any previously added adenosine receptor agonist. The distorted E/c curves were constructed with the same agonists, but in the presence of a single, known, previously administered concentration of the agonist used for the E/c curve. The distorted E/c curves were started when the effect of the surplus agonist concentration had fully developed. The surplus concentration (and its effect) was disregarded during the evaluation to make it a distorting concentration.

All E/c curves, mentioned and presented in this study, were produced from individual E/c curves (i.e., one atrium–one curve) through averaging them within the experimental groups by the curve fitting software. The averaging process preserved the individuality of the effect values (as “replicate Y values”) (Curve Fitting Guide, 2024).

Each intact - distorted E/c curve pair (Szabo et al., 2019a) held information about one distorting concentration near to the EC₅₀ of the given agonist, while each intact - distorted E/c curve family (Gesztelyi et al., 2004) informed about three distorting concentrations spanning two orders of magnitude including the EC₅₀ of the agonist used.

2.2 Regression manners dealing with models, variable parameters and sharing

To characterize the intact E/c relationship, the Hill model was chosen (Gesztelyi et al., 2012a):

$E=E_{\mathit{max}}\bullet \frac{c^n}{c^n+{{EC}_{50}}^n}=\frac{E_{\mathit{max}}}{1+{\left(\frac{{EC}_{50}}{c}\right)}^n}(1)$

To fit the intact E/c curves, instead of the classical form of the Hill equation (the left-hand expression of the Eq 1), a simplified version (the right-hand expression of the Eq 1) was used, in a form where all concentrations (c and EC₅₀) were expressed as a logarithm (10^logc and 10^logEC50, respectively):

$E=\frac{E_{\mathit{max}}}{1+10^{n\cdot \left(\mathit{log}E{C}_{50}-\mathit{log}c\right)}}(2)$

where: c: the concentration of the agonist that was administered during the construction of the E/c curve; E: the effect of c (considered to be intact because of the lack of any distorting factor); E_max: the maximal effect (achievable with the given agonist in the system in question); EC₅₀: the agonist concentration producing half-maximal effect; n: the Hill coefficient (slope factor). During curve fitting, in every case, logc was the independent variable and logEC₅₀ served as a parameter (see Appendix A).

To determine the distorting concentration, RRM was applied. The basic equation of RRM is as follows (Gesztelyi et al., 2004; Grenczer et al., 2010a):

$E^{\prime }=100-\frac{100\bullet \left(100-E\right)}{100-E_x}(3)$

where: E’: the distorted effect that was calculated as if it had exclusively been the effect of c, regardless of the presence of c_x (c_x, the distorting concentration, is always attributed to the same agonist as c by RRM, and in this study, c_x and c belonged indeed to the same agonist); E: the intact counterpart of E′, i.e., the effect that properly reflects the co-action of c and c_x; E_x: the effect of c_x alone (that is also intact, similarly to E).

The model of RRM was derived from the fusion of Eqs 1, 3 (a procedure similar to the combination of the Hill model and the Schild equation: Waud et al., 1978; Motulsky and Christopoulos, 2004; Gesztelyi et al., 2012b). The model of RRM can be expressed in several, algebraically equivalent forms, five of which, relevant for the present investigation, are as follows:

$E^{\prime }=100-\frac{100\cdot \left(100-E_{\mathit{max}}\bullet \frac{{\left(c_x+c\right)}^n}{{\left(c_x+c\right)}^n+{{EC}_{50}}^n}\right)}{100-E_{\mathit{max}}\bullet \frac{{c_x}^n}{{c_x}^n+{{EC}_{50}}^n}}(4)$

$E^{\prime }=100-\frac{100\cdot \left(100-\frac{E_{\mathit{max}}}{1+{\left(\frac{{EC}_{50}}{c_x+c}\right)}^n}\right)}{100-\frac{E_{\mathit{max}}}{1+{\left(\frac{{EC}_{50}}{c_x}\right)}^n}}(5)$

$E^{\prime }=100-\frac{100\cdot \left(100-E_{\mathit{max}}\bullet \frac{{\left(c_x+{10}^{logc}\right)}^n}{{\left(c_x+{10}^{logc}\right)}^n+{10}^{n\bullet \mathit{log}{EC}_{50}}}\right)}{100-E_{\mathit{max}}\bullet \frac{{c_x}^n}{{c_x}^n+{10}^{n\bullet \mathit{log}{EC}_{50}}}}(6)$

$E^{\prime }=100-\frac{100\cdot \left(100-\frac{E_{\mathit{max}}}{1+10^{n\cdot \left(\mathit{log}E{C}_{50}-\mathit{log}\left(c_x+10^{\mathit{log}c}\right)\right)}}\right)}{100-\frac{E_{\mathit{max}}}{1+10^{n\cdot \left(\mathit{log}E{C}_{50}-\mathit{log}\left(c_x\right)\right)}}}(7)$

$E^{\prime }=100-\frac{100\cdot \left(100-\frac{E_{\mathit{max}}}{1+10^{n\cdot \left(\mathit{log}E{C}_{50}-\mathit{log}\left(10^{\mathit{log}{c}_x}+10^{\mathit{log}c}\right)\right)}}\right)}{100-\frac{E_{\mathit{max}}}{1+10^{n\cdot \left(\mathit{log}E{C}_{50}-\mathit{log}{c}_x\right)}}}(8)$

where: c: the concentration of the agonist administered during the construction of the E/c curve; c_x: the distorting concentration; E’: the distorted effect; E_max, logEC₅₀ and n: the parameters from the Hill equation (see Eqs 1, 2).

The Eq 6 stems from the Eq 4, by expressing c and EC₅₀ in logarithmic form (as 10^logc and 10^logEC50). In turn, Eqs 7–8 are derived from the Eq 5, the version of the Eq 4 simplified by a division by c_x, also by expressing c and EC₅₀ in logarithmic form. While Eqs 6, 7 have retained c_x from the previous equations, the Eq 8 contains logc_x as a parameter (Table 1).

Table 1

Table 1. The five algebraically equivalent forms of the model of RRM presented in this study (Eqs 4–8), and three of them used for curve fitting (Eqs 6–8), highlighted in bold. The comparison of Eqs 6, 7 provides information about the influence of the complexity of the fitted equation as well as about the effect of the division by c_x (as c_x should be allowed to be zero, theoretically it should not be used for division). In turn, the comparison of Eqs 7, 8 informs about the impact of the fitted form of the distorting concentration (c_x or logc_x). c (and logc): the concentration of the agonist administered during the construction of the E/c curve (and its logarithm), the independent variable in the model of RRM; EC₅₀ (and logEC₅₀): the half-maximal effective concentration of the agonist administered during the construction of the E/c curve (and its logarithm), a parameter in the model of RRM; c_x (and logc_x): the concentration of the distorting agonist (and its logarithm), the obligate variable parameter in the model of RRM.

Thus, Eqs 6, 7 gave the sought estimate itself, whereas the Eq 8 provided the estimate after taking the antilog of its best-fit value. In Eqs 6-8, E_max, logEC₅₀ and n were parameters in addition to c_x or logc_x, while logc and E’ were the independent variable and the dependent variable, respectively (see Appendix A).

It should be emphasized that, during its derivation, the Eq 6 has not undergone a division by c_x, and it contains c_x as a parameter. Consequently, the Eq 6 allows, in all ways, c_x to be zero. In turn, the Eq 7, although it incorporates c_x as a parameter, stems from a division by c_x. Thus, it can be regarded theoretically questionable whether the Eq 7 allows c_x to be zero. Nevertheless, the Eq 7 is simpler than the Eq 6. Eventually, the Eq 8 comes from a division by c_x and includes logc_x, so it does not permit c_x to be zero at all (Table 1).

To fit the intact and distorted E/c curves, the Eqs 6-8 were used in three ways: individually, globally with one model, and globally with two models.

During the individual regression, one of Eqs 6-8 was fitted to one E/c curve at once (either distorted or intact), in a manner that the Hill parameters (E_max, logEC₅₀ and n) were constrained to constant values. These values were provided by a previous fitting of the Eq 2 to the corresponding (or the same) intact E/c curve. So, during individual regression, of the four parameters in Eqs 6-8, only one was variable (logc_x or c_x).

During global regression, all parameters of Eqs 6-8 were variable (E_max, logEC₅₀, n and logc_x or c_x). When using one model, the related intact and distorted E/c curves were simultaneously fitted to one of Eqs 6-8 with the Hill parameters shared between (or among) the curves. When using two models, the related intact and distorted E/c curves were also simultaneously fitted but in a way that the Eq 2 was fitted to the intact E/c curve, while one of Eqs 6-8 was fitted to the related distorted E/c curve(s), with Hill parameters shared between the equations as well as between (or among) the curves.

When more than one distorted E/c curve belonged to one intact E/c curve, an additional choice was available for the global regression (irrespective of the number of models): curve fitting at once to all related E/c curves (all-at-once manner), and curve fitting simultaneously to only two related E/c curves, one of which was compulsorily the intact E/c curve (pairwise manner). Both options were carried out.

During the individual and one-model global ways of regression, c_x or logc_x values were determined for the intact E/c curves as well, as an additional control (since the expected value for the estimate was zero). The two-model global fitting provided no estimates for the intact E/c curves because they were fitted only to the Hill model (Eq 2).

2.3 Regression manners addressing data scatter

The above-mentioned ways of regression were further combined with ordinary or robust fitting, addressing the distribution of the scatter of data points around the best-fit curve. Upon robust fitting, only accuracy (but not precision) could be compared, because robust regression prevented obtaining 95% CIs and 95% confidence and prediction bands (Curve Fitting Guide, 2024).

2.4 Data processing and presentation

Curve plotting and fitting was performed with GraphPad Prism 9.5.1 for Windows (GraphPad Software Inc., La Jolla, CA, USA). Some operations were made with Microsoft Excel for Microsoft 365 (Microsoft Co., Redmond, WA, USA).

Accuracy of the regression was judged by the distance of the estimate from the corresponding known distorting concentration.

Precision of the regression was characterized by the width of the 95% CI of the best-fit value (c_x or logc_x) addressing the distorting concentration. In addition, precision of the curve fitting and precision of the E/c curve data were characterized by the distance of the 95% confidence and prediction bands, respectively, from the corresponding best-fit curve.

Further information supplied by a 95% CI was the position of the best-fit value within it (Curve Fitting Guide, 2024). If the best-fit value was well centered (i.e., the 95% CI was symmetrical or close to it), the parameterizations of the model were regarded appropriate.

When setting the way in which the software checks how well the experimental data define the model, both available options were used, sc. “identifying ambiguous fits” and “identifying unstable parameters” (Curve Fitting Guide, 2024). The reason for this was that more estimates of the distorting concentration could be obtained with “identifying ambiguous fits” (an option dealing with the extent of parameter intertwining) than with “identifying unstable parameters” (the default option). When RRM was implemented with “identifying unstable parameters”, some estimates (labelled as “ambiguous” by the former option) were missing, while the rest had practically the same values and 95% CIs.

For almost all cases, the default option for computing 95% CIs was “asymmetrical” (Curve Fitting Guide, 2024). Where not (during defining the models to fit), the “asymmetrical” option was chosen. For every setting not addressed above, the default option was used.

3 Results

3.1 The influence of the form of the fitted equation (regarding its complexity and the expression of the distorting concentration)

Upon individual fitting and two-model global regression for the distorted E/c curves, all equations provided practically the same estimates (not affected by the ordinary or robust way of fitting). Upon individual regression for the intact E/c curves, the estimates, while being somewhat different for the different equations, showed small, close-to-zero or zero values (as expected). For the sake of simplification, only data provided by Eqs 7, 8 have been presented throughout this paper (cf. Table 2, 3).

Table 2

Table 2. The logc_x best-fit values (provided by the Eq 8), with 95% confidence intervals (95% CI) and antilog values (c_x, converted to nmol/L, as estimates, highlighted in bold), obtained with the receptorial responsiveness method (RRM) from data of six groups (see column headers) of a previous ex vivo study (Szabo et al., 2019a), using three fitting ways combined with two another fitting options (see row headers). na: not applicable; nM: nmol/L; CPA: N⁶-cyclopentyladenosine; NECA: 5′-(N-ethylcarboxamido)adenosine; CHA: N⁶-cyclohexyladenosine

Table 3

Table 3. The c_x best-fit values (provided by the Eq. 7) converted to nmol/L (as estimates, highlighted in bold), with 95% confidence intervals (95% CI), obtained with the receptorial responsiveness method (RRM) from data of six groups (see column headers) of a previous ex vivo study (Szabo et al., 2019a), using three fitting ways combined with two another fitting options (see row headers). na: not applicable; nM: nmol/L; CPA: N⁶-cyclopentyladenosine; NECA: 5′-(N-ethylcarboxamido)adenosine; CHA: N⁶-cyclohexyladenosine

Upon one-model global regression, however, the fitting of the Eq 6 failed completely, it provided no or obviously wrong estimates (accompanied by error messages from the curve fitting software). The estimates yielded by the Eqs 7, 8 differed, although to a moderate extent (regardless the use of ordinary or robust fitting). The estimates for the intact E/c curves provided by the Eq 8 were clearly better (i.e., closer to zero) than those yielded by the Eq 7 (cf. Table 2, 3).

According to the 95% CIs as well as 95% confidence and prediction bands (being accessible only with ordinary fitting), there was no considerable difference in the precision between the models containing c_x and logc_x (if the fitting was possible) (cf. Table 2, 3, furthermore the adjacent left and right panels of Figures 1–3). As a surprising exception, during one-model global fitting to NECA E/c curves, confidence and prediction bands were only displayable (i.e., only then were they not undisplayably poor) when using the Eq 7 (Figure 2).

Summarizing, in general, the Eq 8 was better than the Eq 7, while the Eq 7 was better than the Eq 6. So, three conclusions could be drawn: 1) the use of c_x (instead of logc_x) in the model of RRM did not improve either accuracy or precision of the estimation; 2) the simpler model of RRM worked more reliably than the more complex, although algebraically equivalent model; 3) avoiding division by c_x was of no benefit for the estimation (even in the case of the intact E/c curves).

3.2 The influence of the number and use of models

Accuracy was the best in the case of individual fitting, which was followed by the other options with only moderate differences. When using the model of RRM with logc_x, accuracy was not affected at all by the fact whether one- or two-model global fitting was chosen. Accuracy was slightly influenced by whether ordinary or robust regression was carried out, and it was moderately influenced by whether c_x or logc_x was in the fitted model (this latter observation was especially true in the case of one-model global fitting). All estimates for the intact E/c curves showed small values, especially when using logc_x to fit (Table 2, 3).

According to the 95% CIs, precision was the best upon individual fitting, but only when the distorted E/c curves were assessed. Precision was acceptable when using two-model global fitting (always), whereas it was unquantifiably poor in the case of one-model global regression (always) and individual fitting (when assessing intact E/c curves) (Table 2, 3). For the one-model global regression, the cause of the ambiguous fit was that c_x and logc_x in the Eqs 7, 8, respectively, were highly intertwined (i.e., showed strong correlation) with the other parameters. The use of c_x (in the Eq 7) substantially increased the correlation even between parameters other than c_x. Moreover, the use of the Eq 6 (also containing c_x) prevented the successful curve fitting during one-model global regression.

The results provided by the 95% confidence and prediction bands were supportive of those yielded by the 95% CIs, with one exception. Namely, the 95% confidence bands for the intact E/c curves, obtained with individual fitting, seem to be surprisingly narrow, when considering the very wide 95% CIs for the corresponding best-fit values (Table 2, 3; Figure 1). This phenomenon was especially conspicuous when using the model with logc_x (Table 2; Figure 1, left panels). This might be due to that the best-fit values involved (both c_x and logc_x) were very small (as expected in the case of intact E/c curves), so even a relatively wide range around these c_x and logc_x values could appear to be narrow. Therefore, these narrow confidence bands might be misleading, while the much wider 95% prediction bands, holding information about the uncertainty of the E/c curve data in addition to the uncertainty of the curve fitting, could be considered more illustrative (Figure 1).

Overall, the two-model global regression could be regarded as the most precise, but the individual regression was not far behind, especially if considering that the estimates for the intact E/c curves served only as an additional control (Table 2, 3; Figures 1–3).

Regarding the convenience of use, the global way of regression was undoubtedly superior to the individual one. Specifically, the one-model global fitting was the most convenient (because it did not require creating a multiline model in the software; see Appendix A). In terms of convenience, it had no relevance whether c_x or logc_x was fitted and in which equation, furthermore whether ordinary or robust regression was performed.

In addition, every best-fit value, for which 95% CI could be computed, was well centered in it, indicating the proper parametrizations of the fitted models.

3.3 The influence of the number of distorted E/c curves fitted at once

It is worth noting in advance that, upon global regression, the assessment of an E/c curve family offers an extra choice for the implementation of RRM: in addition to all-at-once fitting, it is possible to fit in a pairwise manner (that increases the impact of the intact E/c curve on the results). For simplicity, only Eqs 2, 8 were used to evaluate the E/c curve families.

In terms of accuracy, all modes of regression provided similarly good estimates for the greater two of the three distorting concentrations (within a given E/c curve family). However, the estimates for the smallest distorting concentration were acceptable only when individual regression or global regression with pairwise fitting was performed. The estimates for the intact E/c curves (where relevant) were close to zero (Table 4, 5).

Table 4

Table 4. The logc_x best-fit values (yielded by the Eq 8), with 95% confidence intervals (95% CI) and antilog values (c_x, converted to nmol/L, as estimates, highlighted in bold), obtained with the receptorial responsiveness method (RRM) from data of four groups (see column headers) of a previous ex vivo study (Gesztelyi et al., 2004), using three fitting ways combined with two another fitting manners (see row headers). In the case of global regression, curve fitting was carried out by pairing each Distorted group with the related Intact one (“pairwise technique”). na: not applicable; nM: nmol/L; CPA: N⁶-cyclopentyladenosine

Table 5

Table 5. The logc_x best-fit values (provided by the Eq 8), with 95% confidence intervals (95% CI) and antilog values (c_x, converted to nmol/L, as estimates, highlighted in bold), obtained with the receptorial responsiveness method (RRM) from data of four groups (see column headers) of a previous ex vivo study (Gesztelyi et al., 2004), using the two manners of global regression combined with two another fitting options (see row headers). Curve fitting was carried out simultaneously to all groups (“all-at-once technique”). na: not applicable; nM: nmol/L; CPA: N⁶-cyclopentyladenosine

Precision achieved with global regression carried out with pairwise technique was similar in all respects to that described in the previous subsection for the Eq 8 (Table 4, 5; Figure 4). Upon global regression performed in all-at-once manner, however, precision was unquantifiably poor (irrespective of the number of models). Thus, the two-model global regression, implemented with pairwise fitting, was considered the most precise, but it was closely followed by the individual regression.

Figure 4

Figure 4. The implementation of RRM with one-model global regression (panel above), individual regression (left panel below) and two-model global regression (right panel below) combined with ordinary least-squares fitting, using exclusively the common logarithm of the concentration to be determined (logc_x) as a parameter. The procedure was performed on concentration-effect (E/c) curve families, each consisting of one “Intact” and three “Distorted” E/c curves, obtained from an earlier ex vivo study (Gesztelyi et al., 2004). Both ways of global regression were performed in pairwise manner (i.e., only two E/c curves at once, one of which was always the “Intact” E/c curve). The x-axis shows the common logarithm of the molar concentration of CPA, a synthetic A₁ adenosine receptor full agonist, administered for constructing the E/c curves. The y-axis indicates the effect (the percentage decrease in the initial contractile force of isolated, paced guinea pig left atria). The “Intact” E/c curves (black symbols) were generated conventionally, while the “Distorted” E/c curves (symbols filled with a color from dark gray to white) were constructed in the presence of a single, disregarded surplus CPA concentration (which latter increased as the color of the symbol lightened). The symbols indicate the responses to CPA averaged within the groups (±SEM). The continuous lines show the best-fit curves of the fitted Eqs 2, 8 (chosen according to the way of fitting and to the type of the E/c curve). The thick dotted lines indicate the 95% confidence bands, while the thin dotted lines denote the 95% prediction bands (if any). RRM: receptorial responsiveness method; CPA: N⁶-cyclopentyladenosine.

Regarding manageability, the one-model global regression performed with all-at-once fitting was the most convenient to use, it was followed by the one-model global regression with pairwise fitting, then the two-model global regression with all-at-once fitting, and then the two-model global regression with pairwise fitting. Eventually, the most complicated way to perform RRM was the individual regression in this case as well.

Every best-fit value was well centered in its 95% CI (if any), thus the parametrizations of the models were found to be appropriate here as well.

The main results of the study have been summarized in Table 6.

Table 6

Table 6. Comparison of the results of the receptorial responsiveness method (RRM) performed with individual or global fitting and all-at-once or pairwise techniques. These conclusions have been drawn for the model containing logc_x (generally found to be better than models including c_x), but roughly the same is true for the models with c_x.

4 Discussion

In the present study, different implementations of RRM were tested by estimating known concentrations of three stable adenosine receptor agonists, NECA, CPA and CHA, in the vicinity of the A₁ adenosine receptors located in the guinea pig atrial myocardium.

RRM is a simple (i.e., non-multiple) nonlinear regression-based procedure with a unique model containing two variables (c or its logarithm logc, and E′) plus at least one variable parameter (c_x or its logarithm logc_x). The c is the concentration of an agonist that is administered to generate an E/c curve, while E’ is the effect (related to the E/c curve) that is partly (but not completely) evoked by c in a biological system. As for the obligate variable parameter (c_x or logc_x), its role is, in the broadest sense, to quantify something, which (before the generation of the E/c curve) has decreased the responsiveness of the given biological system. This quantification is made by RRM with the concentration of the agonist used for the E/c curve that is capable of producing the same reduction in the responsiveness as the original evoking factor. In a simple case, the “something” to be quantified is a single, constant concentration of an agonist (called “distorting agonist”, the quantification of which is thus the goal of RRM). In turn, in the simplest case, the agonist used for the E/c curve and the distorting agonist are the same, a case when RRM directly estimates the concentration of the distorting agonist (as c_x or as logc_x) (Gesztelyi et al., 2004; Grenczer et al., 2010a; 2010b).

Although RRM deals with an inherently linear issue (see the Eq. 3), the relationship between a concentration (or dose) and a biological effect is typically nonlinear. Thus, we arranged the model of RRM (intended for curve fitting) to be nonlinear via combining the basic equation of RRM with a nonlinear receptor function model (being generally the Hill equation). This way, the model of RRM became suitable to evaluate raw (or minimally transformed) E/c data (Gesztelyi et al., 2004; Grenczer et al., 2010a). This yields more reliable results because the more the data to be fitted are transformed, the greater the risk that the biological variability and measurement errors in the raw data will significantly bias the results (Motulsky and Christopoulos, 2004).

Beyond the simple concentration estimation (such as the determination of the interstitial accumulation of the endogenous adenosine during nucleoside transport blockade: Karsai et al., 2006; 2007), RRM also meets other challenges. RRM is suitable to correct E/c curves of adenosine receptor agonists for the distortion caused by a change in the endogenous adenosine level. Accordingly, RRM helped to uncover the slight positive inotropic effect of EHNA that resulted probably in part from phosphodiesterase type 2 (PDE2) inhibition and partly from an unclear mechanism characteristic of potent adenosine deaminase inhibitors (Gesztelyi et al., 2003; Kemeny-Beke et al., 2007; Pak et al., 2015). The use of RRM also enabled the discovery of a hitherto unknown property of FSCPX, a well-established irreversible A₁ adenosine receptor antagonist, namely, the ability to decrease the interstitial adenosine concentration (via an unknown mechanism), in a series of in silico (Zsuga et al., 2017; Szabo et al., 2019b) and ex vivo investigations (Erdei et al., 2018; Viczjan et al., 2021). Furthermore, RRM contributed to the assessment of the A₁ adenosine receptor reserve for the direct negative inotropic effect of adenosine, an agonist difficult to quantify ex vivo (and in vivo) due to its short half-life (Kiss et al., 2013; Pak et al., 2014; Zsuga et al., 2017). In addition, RRM was used to decide which one of two possible actions of cannabidiol, a widely used phytocannabinoid, can be responsible for the enhancement of the adenosinergic activity (this was found to be the nucleoside transport blockade rather than direct A₁ adenosine receptor agonism) (Viczjan et al., 2022).

The basic equation of RRM (Eq 3) describes the link between a distorted effect (measured in a system with decreased responsiveness) and the corresponding intact effect (determined in the same system in its intact state). Encompassing this relationship, the model of RRM (Eqs 4–8) connects E’, the distorted effect, with c, the known concentration of an agonist (or with logc, its logarithm) and with the cause of the decreased responsiveness that is expressed as c_x, a distorting concentration (or logc_x, its logarithm). Importantly, for simplicity, c_x is attributed to the same agonist as c. The model of RRM incorporates a quantitative model of the receptor function as well and thereby uses its parameters (herein and mostly it is the Hill model: Eq 2). It should be noted that, from a theoretical point of view, the model of RRM can be considered oversimplified when applied to data resulted from the co-action of two different agonists. This circumstance might be responsible for some of the shortcomings of this method (Grenczer et al., 2010a; Grenczer et al., 2010b).

As input for RRM, at least two (types of) E/c curves are needed: an “intact” E/c curve that displays the effect of the agonist used to generate the E/c curves in the system possessing its intact responsiveness, and a “distorted” one that shows the effect of this agonist in the system with decreased response capacity. The distorted E/c curve should be fitted to the model of RRM using some information obtained from the related intact E/c curve (Gesztelyi et al., 2004; Grenczer et al., 2010a; Szabo et al., 2019a). The way this information is transferred is what distinguishes the individual (local) regression, the one-model global regression and the two-model global regression (as three options for implementing RRM). In the present investigation, one of our goals was to find out which of these options is optimal for RRM. In addition, the two-model global fitting was also used as an alternative tool to clarify the significance of logc_x in the model of RRM (in terms of fitting the intact E/c curves, see below).

A surprising property of regression (at least for non-mathematicians) is that the fitting of different but algebraically equivalent expressions can lead to somewhat different results (Young, 2017). Consequently, another goal of this study was to explore whether the simplification of the algebraic form of the RRM’s model can improve the outcome of RRM. Because the simplification used here included a division by c_x, a parameter that should otherwise be allowed to be zero, results obtained this way could reflect the consequences of both the simplification and the division by c_x, which latter operation theoretically excludes the possibility of c_x to be zero.

Not independently from the fact that different but equivalent equations can provide different results, it was recommended that a parameter following log-normal (Galton) distribution (e.g., all kinds of concentration) should be used as a common logarithm in models intended for curve fitting, in order to obtain a symmetrical CI (Motulsky and Christopoulos, 2004). In the model of RRM, the distorting concentration can also be expressed as either a logarithm (logc_x) or a numerus (c_x), options that may affect the results of the estimation. However, as the software used for curve fitting in the present study can handle asymmetrical CIs as well, it seemed reasonable to revisit, which form of the distorting concentration is better for fitting. Therefore, in the present study, three different but algebraically equivalent models of RRM (Eqs 6-8) were used to explore the influence of the following three interventions on the accuracy and precision of RRM: 1) simplification of the model of RRM; 2) division by c_x in the RRM’s model (inseparable from the simplification); and 3) changing the form of the distorting concentration (c_x or logc_x) in the RRM’s model.

The regression ways mentioned so far can be combined with two additional fitting options, i.e., ordinary and robust regression, to consider the distribution of the scatter of data points around the best-fit curve. This distribution can be closer to either Gaussian or Lorentzian distribution (the two extremes of the t distribution regarding the number of degrees of freedom), which conditions require ordinary or robust regression, respectively. It should be noted that the use of robust regression limits the usefulness of the estimation as it prevents the determination of precision (Curve Fitting Guide, 2024).

Furthermore, if more than one distorted E/c curve belongs to one intact E/c curve, and global regression is chosen, there will be two possibilities to perform RRM: fitting all related E/c curves at once or fitting them pairwise (i.e., each distorted E/c curve separately with the intact E/c curve). These two options may also lead to different results to be addressed.

In the present study, accuracy of RRM has generally been found to be acceptable for all kinds of regression used herein (Table 2–5). Nevertheless, there were three noteworthy exceptions: 1) the one-model global fitting of the Eq 6 (the non-simplified model of RRM with c_x); 2) the one-model global fitting for the intact E/c curves with the Eq 7 (the simplified model with c_x) (Table 3); and 3) the global fitting performed in all-at-once manner to determine the smallest distorting concentration (obtained from Gesztelyi et al., 2004) (Table 5). From these observations, three major conclusions can be drawn: 1) in the model of RRM, the use of c_x is rather a disadvantage than an advantage; 2) the algebraically simplified model is better for RRM (even despite the theoretical concern of dividing by c_x, a parameter that should be allowed to be zero); and 3) a small distorting concentration is a challenge for RRM, especially if it is determined by fitting an E/c curve family in all-at-once manner. These conclusions are detailed as follows.

Ad 1): Rewriting the RRM’s model to replace logc_x with c_x did not improve either accuracy or precision of the estimation. On the contrary, in some cases, the use of c_x worsened precision by increasing the correlation between some parameters, moreover, the combination of the Eq 6 (the original complex model entirely allowing c_x to be zero) with the one-model global fitting made the determination impossible. This finding has refuted our previous assumption about the poor performance of the one-model global regression implemented with a model containing logc_x (Szabo et al., 2019a). So, foibles of RRM in any case do not stem from the fit of logc_x that would hinder the correct determination of the zero value of c_x for the intact E/c curves. This conclusion has been confirmed by another finding of the present study. Namely, the two-model global regression (with a model containing either logc_x or c_x) did not improve the accuracy of RRM in comparison to either the one-model global fitting or the individual regression (Table 2, 3), although it has proven to be the most precise (Table 2, 3; Figures 1–3).

Ad 2): The problems in the estimation observed in some cases may be the results of the relative complexity of the model of RRM (even for the algebraically simplified model expressed by the Eq 8). The greater the complexity of a regression model, the greater the degree of correlation that can occur between the parameters of the model (Motulsky and Christopoulos, 2004; Curve Fitting Guide, 2024). Thus, for RRM, a model as simple as possible should be used (Table 1). In addition, in cases when the agonist used for the E/c curves and the distorting agonist are different, the causal role of the theoretical oversimplification included in the model of RRM may also be significant in terms of accuracy as well as precision (Grenczer et al., 2010a; 2010b).

Ad 3): This conclusion is consistent with our earlier in silico finding that too small and too large concentrations might be estimated with substantial inaccuracy by means of RRM (Grenczer et al., 2010a; 2010b). (A distorting concentration can be regarded small or large if its effect, namely, E_x in the Eq 3, is small or large relative to the maximal effect of the agonist used for the E/c curves, i.e., E_max in Eqs 4–8, respectively.) However, for E/c curve families, accuracy (when a small concentration was to be determined) and precision (in every case) of the global regression could be substantially improved by fitting the appropriate E/c curves pairwise (cf. Table 4, 5). Thus, despite the acceptable estimation for most cases and the good usability, we do not recommend fitting globally to more than two E/c curves at once (for RRM). In line with this, unless otherwise indicated, we discuss those results of global regression that were obtained with pairwise fitting technique.

As for the further fitting options, in terms of accuracy, the individual fitting proved to be the best, followed closely by the global regression (when fitting the model with logc_x, irrespective of the number of models). In turn, regarding precision (i.e., the reliability of the estimates), the two-model global fitting was the best, followed closely by the individual regression and, from afar, by the one-model global fitting (Table 2–5). As regards manageability, the one-model global regression was the best, followed by the moderately complicated two-model global fitting, and then by the most complicated individual regression (Table 6).

The poor precision of the one-model global regression has spectacularly been indicated by the fact that (of course upon ordinary fitting) neither confidence intervals nor confidence and prediction bands could be obtained for the results (an interesting exception was seen for the NECA E/c curves when fitting the Eq 7) (Figure 2). The source of this uncertainty may be that the parameter estimating the distorting concentration (c_x or logc_x) was greatly dependent on the other parameters when fitting the intact E/c curves. The other occasion, when similarly high level of parameter correlation occurred, was the attempt to determine the smallest distorting concentration by fitting simultaneously more than two E/c curves (affecting both types of global regression).

In favor of the oldest way to perform RRM, it should be noted that the individual regression imposes the fewest requirements for the curve fitting software used. Nevertheless, there are nowadays much software with great capabilities to choose from.

Additionally, most differences among the various implementations of RRM have been found in terms of precision and convenience of use. Accordingly, the choice between ordinary or robust regression had little impact in this study, as the robust regression prevented calculating confidence intervals for the best-fit values as well as confidence and prediction bands for the best-fit curves. Besides, the robust regression did not affect the convenience (manageability) of curve fitting at all (Table 2–5).

As a final conclusion: the use of RRM is recommended with 1) the simplified model containing logc_x; 2) two-model global fitting (pairwise, if appropriate), or possibly individual fitting; 3) ordinary least-squares regression, or possibly robust regression.

Data availability statement

The original contributions presented in the study are included in the article, further inquiries can be directed to the corresponding author.

Ethics statement

The two animal studies, sources of the data reevaluated in the present investigation, were approved by the Local Ethics Committee of Animal Research, University of Debrecen, Hungary (codes: DEMÁB 45/2001 and 25/2013/DEMÁB). The study was conducted in accordance with the local legislation and institutional requirements.

Author contributions

IO: Investigation, Writing–original draft. GV: Investigation, Writing–review and editing. TE: Investigation, Writing–review and editing. BT: Investigation, Writing–review and editing. VT: Investigation, Writing–review and editing. JZ: Writing–review and editing, Data curation. MS: Writing–review and editing, Investigation. ZS: Writing–review and editing, Funding acquisition. BJ: Writing–review and editing, Resources. RG: Writing–review and editing, Conceptualization, Formal Analysis.

Funding

The author(s) declare financial support was received for the research, authorship, and/or publication of this article. This research was supported by the European Union and the State of Hungary under grant numbers TKP2021-EGA-18 and TKP2020-NKA-04. (Project no. TKP2021-EGA-18 has been implemented with the support provided by the Ministry of Culture and Innovation of Hungary from the National Research, Development and Innovation Fund, financed under the TKP2021-EGA funding scheme. Project no. TKP2020-NKA-04 has been implemented with the support provided from the National Research, Development and Innovation Fund of Hungary, financed under the 2020-4.1.1-TKP2020 funding scheme.) In addition, this research was supported by the University of Debrecen Program for Scientific Publication.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

Andresen, B. T. (2011). A pharmacological primer of biased agonism. Endocr. Metab. Immune Disord. Drug Targets 11, 92–98. doi:10.2174/187153011795564179

PubMed Abstract | CrossRef Full Text | Google Scholar

Anjos, R. S. D., Nóbrega, R. S., Ferreira, H. D. S., Lacerda, A. P., and Sousa-Neves, N. (2020). Exploring local and global regression models to estimate the spatial variability of Zika and Chikungunya cases in Recife, Brazil. Rev. Soc. Bras. Med. Trop. 53, e20200027. doi:10.1590/0037-8682-0027-2020

PubMed Abstract | CrossRef Full Text | Google Scholar

Belanger, M. C., Anbaei, P., Dunn, A. F., Kinman, A. W. L., and Pompano, R. R. (2020). Spatially resolved analytical chemistry in intact, living tissues. Anal. Chem. 92, 15255–15262. doi:10.1021/acs.analchem.0c03625

PubMed Abstract | CrossRef Full Text | Google Scholar

Curve Fitting Guide (2024). Curve fitting Guide. Available at: https://www.graphpad.com/guides/prism/latest/curve-fitting/index.htm (accessed on April 21, 2024).

Google Scholar

Erdei, T., Szabo, A. M., Lampe, N., Szabo, K., Kiss, R., Zsuga, J., et al. (2018). FSCPX, a chemical widely used as an irreversible A₁ adenosine receptor antagonist, modifies the effect of NBTI, a nucleoside transport inhibitor, by reducing the interstitial adenosine level in the Guinea pig atrium. Molecules 23, E2186. doi:10.3390/molecules23092186

CrossRef Full Text | Google Scholar

Franco, R., Rivas-Santisteban, R., Reyes-Resina, I., and Navarro, G. (2021). The old and new visions of biased agonism through the prism of adenosine receptor signaling and receptor/receptor and receptor/protein interactions. Front. Pharmacol. 11, 628601. doi:10.3389/fphar.2020.628601

PubMed Abstract | CrossRef Full Text | Google Scholar

Fredholm, B. B., Ijzerman, A. P., Jacobson, K. A., Klotz, K. N., and Linden, J. (2001). International Union of Pharmacology. XXV. Nomenclature and classification of adenosine receptors. Pharmacol. Rev. 53, 527–552. PMID:11734617.

PubMed Abstract | Google Scholar

Gesztelyi, R., Kiss, Z., Zsuga, J., Pak, K., Papp, C., Galajda, Z., et al. (2012b). Thyroid hormones decrease the affinity of 8-cyclopentyl-1,3-dipropylxanthine (CPX), a competitive antagonist, for the Guinea pig atrial A(1) adenosine receptor. Gen. Physiol. Biophys. 31, 389–400. doi:10.4149/gpb_2012_043

PubMed Abstract | CrossRef Full Text | Google Scholar

Gesztelyi, R., Zsuga, J., Hajdú, P., Szabó, J. Z., Cseppento, A., and Szentmiklósi, A. J. (2003). Positive inotropic effect of the inhibition of cyclic GMP-stimulated 3',5'-cyclic nucleotide phosphodiesterase (PDE2) on Guinea pig left atria in eu- and hyperthyroidism. Gen. Physiol. Biophys. 22, 501–513. PMID: 15113122.

PubMed Abstract | Google Scholar

Gesztelyi, R., Zsuga, J., Juhász, B., Dér, P., Vecsernyés, M., and Szentmiklósi, A. J. (2004). Concentration estimation via curve fitting: quantification of negative inotropic agents by using a simple mathematical method in Guinea pig atria. Bull. Math. Biol. 66, 1439–1453. doi:10.1016/j.bulm.2004.03.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Gesztelyi, R., Zsuga, J., Kemeny-Beke, A., Varga, B., Juhász, B., and Tosaki, A. (2012a). The Hill equation and the origin of quantitative pharmacology. Arch. Hist. Exact. Sci. 66, 427–438. doi:10.1007/s00407-012-0098-5

CrossRef Full Text | Google Scholar

Grenczer, M., Pinter, A., Zsuga, J., Kemeny-Beke, A., Juhasz, B., Szodoray, P., et al. (2010a). The influence of affinity, efficacy, and slope factor on the estimates obtained by the receptorial responsiveness method (RRM): a computer simulation study. Can. J. Physiol. Pharmacol. 88, 1061–1073. doi:10.1139/y10-078

PubMed Abstract | CrossRef Full Text | Google Scholar

Grenczer, M., Zsuga, J., Majoros, L., Pinter, A., Kemeny-Beke, A., Juhasz, B., et al. (2010b). Effect of asymmetry of concentration-response curves on the results obtained by the receptorial responsiveness method (RRM): an in silico study. Can. J. Physiol. Pharmacol. 88, 1074–1083. doi:10.1139/y10-089

PubMed Abstract | CrossRef Full Text | Google Scholar

Herman, P., and Lee, J. C. (2012). “The advantage of global fitting of data involving complex linked reactions,” in Allostery. Methods in molecular biology. Editor A. Fenton (New York: Springer). doi:10.1007/978-1-61779-334-9_22

CrossRef Full Text | Google Scholar

Ijzerman, A. P., Jacobson, K. A., Müller, C. E., Cronstein, B. N., and Cunha, R. A. (2022). International union of basic and clinical pharmacology. CXII: adenosine receptors: a further update. Pharmacol. Rev. 74, 340–372. doi:10.1124/pharmrev.121.000445

PubMed Abstract | CrossRef Full Text | Google Scholar

Jarpe, M. B., Knall, C., Mitchell, F. M., Buhl, A. M., Duzic, E., and Johnson, G. L. (1998). [D-Arg1,D-Phe5,D-Trp7,9,Leu11]Substance P acts as a biased agonist toward neuropeptide and chemokine receptors. J. Biol. Chem. 273, 3097–3104. doi:10.1074/jbc.273.5.3097

PubMed Abstract | CrossRef Full Text | Google Scholar

Karsai, D., Gesztelyi, R., Zsuga, J., Jakab, A., Szendrei, L., Juhasz, B., et al. (2007). Influence of hyperthyroidism on the effect of adenosine transport blockade assessed by a novel method in Guinea pig atria. Cell Biochem. Biophys. 47, 45–52. doi:10.1385/cbb:47:1:45

PubMed Abstract | CrossRef Full Text | Google Scholar

Karsai, D., Zsuga, J., Juhasz, B., Der, P., Szentmiklosi, A. J., Tosaki, A., et al. (2006). Effect of nucleoside transport blockade on the interstitial adenosine level characterized by a novel method in Guinea pig atria. J. Cardiovasc Pharmacol. 47, 103–109. doi:10.1097/01.fjc.0000196239.51018.a0

PubMed Abstract | CrossRef Full Text | Google Scholar

Kemeny-Beke, A., Jakab, A., Zsuga, J., Vecsernyes, M., Karsai, D., Pasztor, F., et al. (2007). Adenosine deaminase inhibition enhances the inotropic response mediated by A1 adenosine receptor in hyperthyroid Guinea pig atrium. Pharmacol. Res. 56, 124–131. doi:10.1016/j.phrs.2007.04.017

PubMed Abstract | CrossRef Full Text | Google Scholar

Kenakin, T. (2024). Bias translation: the final frontier? Br. J. Pharmacol. 181, 1345–1360. doi:10.1111/bph.16335

PubMed Abstract | CrossRef Full Text | Google Scholar

Kenakin, T., and Pharmacology Primer, A. (2022). Techniques for more effective and strategic Drug discovery. Academic Press. ISBN: 9780323992893.

Google Scholar

Kiss, Z., Pak, K., Zsuga, J., Juhasz, B., Varga, B., Szentmiklosi, A. J., et al. (2013). The Guinea pig atrial A1 adenosine receptor reserve for the direct negative inotropic effect of adenosine. Gen. Physiol. Biophys. 32, 325–335. doi:10.4149/gpb_2013041

PubMed Abstract | CrossRef Full Text | Google Scholar

Klaasse, E. C., Ijzerman, A. P., de Grip, W. J., and Beukers, M. W. (2008). Internalization and desensitization of adenosine receptors. Purinergic Signal 4, 21–37. doi:10.1007/s11302-007-9086-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Knutson, J. R., Beechem, J. M., and Brand, L. (1983). Simultaneous analysis of multiple fluorescence decay curves: a global approach. Chem. Phys. Lett. 102, 501–507. doi:10.1016/0009-2614(83)87454-5

CrossRef Full Text | Google Scholar

Maguire, J. J. (2016). Evidence for biased agonists and antagonists at the endothelin receptors. Life Sci. 159, 30–33. doi:10.1016/j.lfs.2016.02.069

PubMed Abstract | CrossRef Full Text | Google Scholar

Markiewicz, R., Litowczenko, J., Gapiński, J., Woźniak, A., Jurga, S., and Patkowski, A. (2022). Nanomolar nitric oxide concentrations in living cells measured by means of fluorescence correlation spectroscopy. Molecules 27, 1010. doi:10.3390/molecules27031010

PubMed Abstract | CrossRef Full Text | Google Scholar

McNeill, S. M., Baltos, J. A., White, P. J., and May, L. T. (2021). Biased agonism at adenosine receptors. Cell Signal 82, 109954. doi:10.1016/j.cellsig.2021.109954

PubMed Abstract | CrossRef Full Text | Google Scholar

Motulsky, H. J., and Christopoulos, A. (2004). “Fitting models to biological data using linear and nonlinear regression,” in A practical guide to curve fitting (San Diego: GraphPad Software Inc.). Available at: https://cdn.graphpad.com/faq/2/file/Prism_v4_Fitting_Models_to_Biological_Data.pdf (accessed on April 8, 2023).

Google Scholar

Mundell, S., and Kelly, E. (2011). Adenosine receptor desensitization and trafficking. Biochim. Biophys. Acta 1808, 1319–1328. doi:10.1016/j.bbamem.2010.06.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Pak, K., Papp, C., Galajda, Z., Szerafin, T., Varga, B., Juhasz, B., et al. (2014). Approximation of A1 adenosine receptor reserve appertaining to the direct negative inotropic effect of adenosine in hyperthyroid Guinea pig left atria. Gen. Physiol. Biophys. 33, 177–188. doi:10.4149/gpb_2013079

PubMed Abstract | CrossRef Full Text | Google Scholar

Pak, K., Zsuga, J., Kepes, Z., Erdei, T., Varga, B., Juhasz, B., et al. (2015). The effect of adenosine deaminase inhibition on the A1 adenosinergic and M2 muscarinergic control of contractility in eu- and hyperthyroid Guinea pig atria. Naunyn Schmiedeb. Arch. Pharmacol. 388, 853–868. doi:10.1007/s00210-015-1121-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Pavan, B., and Ijzerman, A. P. (1998). Processing of adenosine receptor agonists in rat and human whole blood. Biochem. Pharmacol. 56, 1625–1632. doi:10.1016/S0006-2952(98)00270-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Ramakers, B. P., Pickkers, P., Deussen, A., Rongen, G. A., van den Broek, P., van der Hoeven, J. G., et al. (2008). Measurement of the endogenous adenosine concentration in humans in vivo: methodological considerations. Curr. Drug Metab. 9, 679–685. doi:10.2174/138920008786049249

PubMed Abstract | CrossRef Full Text | Google Scholar

Spatial Data Science (2024). Spatial data science. Available at: https://rspatial.org/analysis/6-local_regression.html (accessed on May 9, 2023).

Google Scholar

Szabo, A. M., Erdei, T., Viczjan, G., Kiss, R., Zsuga, J., Papp, C., et al. (2019b). An advanced in silico modelling of the interaction between FSCPX, an irreversible A1 adenosine receptor antagonist, and NBTI, a nucleoside transport inhibitor, in the Guinea pig atrium. Molecules 24, E2207. doi:10.3390/molecules24122207

PubMed Abstract | CrossRef Full Text | Google Scholar

Szabo, A. M., Viczjan, G., Erdei, T., Simon, I., Kiss, R., Szentmiklosi, A. J., et al. (2019a). Accuracy and precision of the receptorial responsiveness method (RRM) in the quantification of A1 adenosine receptor agonists. Int. J. Mol. Sci. 20, 6264. doi:10.3390/ijms20246264

PubMed Abstract | CrossRef Full Text | Google Scholar

Valcárcel, M. (2000). “The analytical problem,” in Principles of analytical chemistry (Berlin, Heidelberg: Springer). doi:10.1007/978-3-642-57157-2_7

CrossRef Full Text | Google Scholar

Verzijl, D., and Ijzerman, A. P. (2011). Functional selectivity of adenosine receptor ligands. Purinergic Signal. 7, 171–192. doi:10.1007/s11302-011-9232-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Viczjan, G., Erdei, T., Ovari, I., Lampe, N., Szekeres, R., Bombicz, M., et al. (2021). A body of circumstantial evidence for the irreversible ectonucleotidase inhibitory action of FSCPX, an agent known as a selective irreversible A1 adenosine receptor antagonist so far. Int. J. Mol. Sci. 22, 9831. doi:10.3390/ijms22189831

PubMed Abstract | CrossRef Full Text | Google Scholar

Viczjan, G., Szilagyi, A., Takacs, B., Ovari, I., Szekeres, R., Tarjanyi, V., et al. (2022). The effect of a long-term treatment with cannabidiol-rich hemp extract oil on the adenosinergic system of the Zucker diabetic fatty (ZDF) rat atrium. Front. Pharmacol. 13, 1043275. doi:10.3389/fphar.2022.1043275

PubMed Abstract | CrossRef Full Text | Google Scholar

Waud, D. R., Son, S. L., and Waud, B. E. (1978). Kinetic and empirical analysis of dose-response curves illustrated with a cardiac example. Life Sci. 22, 1275–1285. doi:10.1016/0024-3205(78)90096-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Willems, L., Reichelt, M. E., Molina, J. G., Sun, C. X., Chunn, J. L., Ashton, K. J., et al. (2006). Effects of adenosine deaminase and A1 receptor deficiency in normoxic and ischaemic mouse hearts. Cardiovasc Res. 71, 79–87. doi:10.1016/j.cardiores.2006.03.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Young, D. S. (2017). Handbook of regression methods. New York: Chapman and Hall/CRC. doi:10.1201/9781315154701

CrossRef Full Text | Google Scholar

Zhang, X., Reeves, D., Shi, Y., Gimi, B., Nemani, K. V., Perreard, I. M., et al. (2015). Toward localized in vivo biomarker concentration measurements. IEEE Trans. Magn. 51, 1–4. doi:10.1109/TMAG.2014.2324993

PubMed Abstract | CrossRef Full Text | Google Scholar

Zsuga, J., Erdei, T., Szabó, K., Lampe, N., Papp, C., Pinter, A., et al. (2017). Methodical challenges and a possible resolution in the assessment of receptor reserve for adenosine, an agonist with short half-life. Molecules 22, E839. doi:10.3390/molecules22050839

CrossRef Full Text | Google Scholar

Appendix

In this section, the equations have been presented as they were written into the curve fitting software. Consistently, the independent and dependent variables are X and Y here (while elsewhere they are logc and E or E’, respectively).

As Hill model, the software’s built-in equation “sigmoidal dose-response (variable slope)”, with a “bottom” parameter constrained to zero (Curve Fitting Guide, 2024), was used (see the Eq. 2):

Y = Emax/(1 + 10 ^ (n*(logEC50-X)))

Eqs 6 and 7 and 8, respectively, prepared for the individual fitting and one-model global fitting:

Y = 100–100*(100-Emax*(cx+10^X)^n/((cx+10^X)^n+10^(n*logEC50)))/(100-Emax*cx^n/(cx^n+10^(n*logEC50)))

and

Y = 100–100*(100-Emax/(1 + 10^(n*(logEC50-log (cx+10^X)))))/(100-Emax/(1 + 10^(n*(logEC50-log (cx)))))

and

Y = 100–100*(100-Emax/(1 + 10^(n*(logEC50-log (10^logcx+10^X)))))/(100-Emax/(1 + 10^(n*(logEC50-logcx))))

Multiline models involving Eq. 2 combined with Eqs 6 and 7 and 8, respectively, prepared for the two-model global fitting:

Hill = Emax/(1 + 10^(n*(logEC50-X)))

RRM = 100–100*(100-Emax*(cx+10^X)^n/((cx+10^X)^n+10^(n*logEC50)))/(100-Emax*cx^n/(cx^n+10^(n*logEC50)))

<A>Y=Hill

<∼A > Y = RRM

and

Hill = Emax/(1 + 10^(n*(logEC50-X)))

RRM = 100–100*(100-Emax/(1 + 10^(n*(logEC50-log (cx+10^X)))))/(100-Emax/(1 + 10^(n*(logEC50-log (cx)))))

<A>Y=Hill

<∼A > Y = RRM

and

Hill = Emax/(1 + 10^(n*(logEC50-X)))

RRM = 100–100*(100-Emax/(1 + 10^(n*(logEC50-log (10^logcx+10^X)))))/(100-Emax/(1 + 10^(n*(logEC50-logcx))))

<A>Y=Hill

<∼A > Y = RRM

(For more information, see “fitting different models to different data sets” in: Curve Fitting Guide, 2024.)