An intact model for quantifying functional selectivity.

A ligand that acts on a target receptor to activate particular multiple signalling pathways with activity that is distinct from other ligands is termed ligand bias. Quantification of ligand bias is based on applying the operational model to each pathway separately and subsequent calculation of the ligand bias metric (ΔΔlogR). This approach implies independence among different pathways and causes propagation of error in the calculation. Here, we propose a semi-mechanism-based model which allows for receptor selectivity across all the pathways simultaneously (termed the 'intact operational model'). The power of the intact model for detecting ligand bias was evaluated via stochastic simulation estimation studies. It was also applied to two examples: (1) opposing effects of Gi/Gs signalling of α2-adrenergic receptors and (2) simultaneous measurement of arachidonic acid release and inositol phosphate accumulation following 5-HT2C receptor activation. The intact operational model demonstrated greater power to detect ligand bias in the simulation. In the applications, it provided better precision of estimation and identified biased ligands that were missed by analysis of traditional methods. Issues identified in both examples might lead to different interpretations of the data. The intact operational model may elucidate greater understanding of the underlying mechanisms of functional selectivity.

Introduction

Different ligands can differentially activate multiple signalling pathways when coupled to a single receptor. This feature is termed functional selectivity[1]. This complex pharmacological process involves multiple biological steps, such as receptor binding and activation of a range of intracellular signalling pathways. Simplified models have been adopted to delineate this phenomenon, in which each pathway is considered independently[2,3].

In current practice, the most commonly used simplified model is the operational model, proposed in 1983 by Black and Leff[4]. In this approach, the operational model is applied separately to each pathway[5]. Since it is known that the pathways are linked biologically, this constitutes a simplification by assuming the pathways are independent. In modelling terms, this would be described as marginal and thus we denote this as ‘the marginal operational model’. When the operational model is applied in isolation, a composite parameter, R (the transduction coefficient; the ratio of the transducer ratio, τ, and functional affinity, K), is obtained for quantifying the effect of a ligand on a single pathway. Ligand bias (the relative preference of a ligand for a particular pathway) is then computed in a post hoc analysis where a ligand’s transduction coefficient is normalised to that of a reference ligand in order to accommodate system and observational bias[6]. This yields the normalised transduction coefficient ΔlogR. This value is then further normalised such that a ligand’s transduction coefficient for one pathway is related to its ΔlogR for a second pathway. This second normalisation step provides a widely reported metric for ligand bias, ΔΔlogR. Calculating this metric for a given ligand in multiple pairwise pathway comparisons will formulate a signature profile of biases (preferences) across the range of pathways of interest, and these biases are often displayed as a radar plot[7].

The use of the marginal operational model makes the assumption that each signalling pathway is independent. This precludes insights pertaining to functional selectivity from observed phenomena, such as the natural correlation of a ligand’s C50 values in different pathways. In addition, there is propagation of errors associated with the required post hoc parameter manipulation when the marginal operational model is used as the basis for determining ligand bias, due to the requirement for two successive normalisation steps.

In order to reflect the linkages between different signalling pathways, a mechanism-based three-state model has been introduced[8]. Within this model framework, equilibria are linked between different receptor conformational states and there is mutual depletion of these receptor states. This means that C50 values should be exactly the same for different signalling pathways. This strict constraint makes the three-state model less flexible and limits its applicability.

Here, we propose a simple extension of the operational model to allow its use across all pathways simultaneously and provide a direct estimate of ΔΔlogR. This model allows all pathways to be simultaneously modelled. Since we believe this provides mechanistic insight into what we consider to be an intact system, we denote this model as ‘the intact operational model’. This approach may shed light into the underlying mechanism of functional selectivity and act as the first step in the development of a full mechanistic model of functional selectivity. In addition, by allowing data from different pathways to inform the common components of the intact model, then the parameters of this model (e.g. logR) should be more precisely estimated.

The aims of this study are (1) to describe the intact operation model, (2) to evaluate its performance compared to the commonly utilised marginal operational model via theoretical assessment (power analysis) and practical applications (two literature examples) and (3) to demonstrate how a go/no-go workflow decision for ligand bias experiments varies based on the choice of operational model.

The main experimental components of this paper are divided into three parts to reflect progressive flow in our study. In the first part, the derivation of the intact operational model and general aspects of the data analyses are provided as they relate to both theoretical evaluation and practical applications. In the second part, the newly proposed method is theoretically evaluated via stochastic simulation estimation studies. In the third part, the findings from theoretical evaluation are further evaluated via practical applications in two examples. We have kept notation here the same as the convention. We note that R is conventionally used to represent two separate entities, the receptor R (shown in schematics and mass balance equations) and the transduction coefficient. To differentiate these components, when the receptor R is unbound to a ligand, A, we add the subscript ‘ub’ (R), and when referring to the total receptor numbers we use R. The transduction coefficient is indexed (if required) by either a numeric value (if more than one ligand is present), or by the pathway signal (e.g. ‘i’ for inhibitory pathway and ‘s’ for stimulatory pathways).

Part I

Derivation of intact operational model

Marginal operational model

The marginal operational model is illustrated in Fig. 1. In this representation, the different receptor conformational states are assumed to be independent. Within this model framework, the general operational model (Eq. 1, the same as Eq. 5 in van der Westhuizen et al., 2014) can be separately applied to each signalling pathway to fit the data.Here, E is the system maximal response, Basal is the baseline response in the absence of ligand and n is the Hill slope factor. These three parameters are system parameters. The ligand-specific parameters are K (the equilibrium dissociation constant), and R (the transduction coefficient, τ/K). Conventionally, the parameters K and R are transformed into logarithms (i.e., and 10).

Figure 1

The schematic plot of the marginal operational model and the intact operational model.

By definition, the ligand bias metric is calculated as the difference between transduction coefficients from two pathways for the test ligand normalised to those of the reference ligand (Eq. 2):Here, the superscript letter ‘REF’ indicates reference ligand. A subscript number “1” indicates pathway 1 and “2” indicates pathway 2.

The estimated standard error (SE) for test ligand is calculated using Eq. 3 (assuming independence – which for the standard post hoc analysis is reasonable).

Intact operational model

The intact operational model is illustrated in Fig. 1. Here we see that equilibria are linked among different receptor conformational states and there is mutual depletion of these receptor states, though there may exist a theoretically unlimited number of active states. We only consider two active receptor conformations in the current study but this can be generalised to a greater number of states. The intact operational model (Eq. 4 for pathway 1 and Eq. 5 for pathway 2) is derived in detail in Appendix 1. Due to identifiability issues, it is not possible to estimate individual K values given only functional assay data. Hence, the information for those values is integrated into (Eq. 6). Here, is proportional to the harmonic mean of individual K values. Note that is introduced to improve model identifiability but this should not prevent the researchers from estimating the individual K values when more information is at hand.Here, (Eq. 6) refers to the apparent equilibrium dissociation constant. R1 and R2 are the transduction coefficients (the ratio of transducer ratio and the equilibrium dissociation constant) for each of the two pathways.The intact operational model can be generalised to account for non-zero basal (Basal), non-unity slope factor (n), and estimable maximal system response (E), which we refer to as the general intact operational model (Eq. 7 for pathway 1 and Eq. 8 for pathway 2, detailed in Appendix 2). For the purpose of modelling, the parameters , R1 and R2 are transformed into logarithms (i.e., , and ).In order to jointly model all the functional assay data, Eqs 7 and 8 are combined into Eq. 9. Here, I and I are indicator functions. I is equal to 1 for producing effect of pathway 1 and I is equal to 1 for producing effect of pathway 2.

In order to directly estimate ligand bias metric, ΔΔlogR1−2, a re-parameterisation is performed to explicitly incorporate ΔΔlogR1−2 as an estimated parameter into the model. From Eq. 2, logR2 is expressed as Eq. 10:

Then, substituting Eq. 10 into Eq. 9 yields the intact operational model with ΔΔlogR1−2 as a directly estimated parameter (Eq. 11):

Thus, , , , logR1 and ΔΔlogR1−2 are the ligand specific parameters that are incorporated directly into the intact operational model. Hence, the ligand bias metric (ΔΔlogR1−2) can be estimated and the estimated standard error can be obtained from modelling output.

General aspects of data analysis

Model analysis

Parameter estimation was conducted using the software NONMEM 7.3.0 (ICON Development Solutions, Hanover, MD, USA)[9]. Model development was managed using Perl-Speaks-NONMEM 4.5.0, Pirana 2.9.0 and Xpose 4.5.3[10]. Due to the deterministic identifiability issue in the marginal operational model[11], it was not possible to estimate logK values of full agonists from direct fitting to concentration-response curves. We followed the convention to arbitrarily set logK of the most efficacious ligand to 0 to circumvent this problem when it was necessary[11]. Model selection was informed by the minimum objective function value (MOFV)[12]. Model evaluation was based on the visual inspection of goodness of fit plots: (1) observed responses vs. predicted responses, (2) weighted residuals vs. predicted responses and (3) weighted residuals vs. ligand concentrations.

Test for ligand bias

Statistical analysis was performed using a Wald test on the ligand bias metric ΔΔlogR. Since the statistical significance of ligand bias for each test ligand was screened individually, the data from the reference ligand was used multiple times for constructing the ligand bias metrics ΔΔlogR for each ligand. To avoid the issue with multiple comparisons (inflation of type I error), the Bonferroni correction[13] was used where the level of significance was set to α/n, where α was the desired overall significance level (i.e., 0.05 in this study) and n was the number of hypothesis tested.

Part II

Theoretical evaluation of the intact model

Since the intact operational model provides a more mechanism-based description of the biological system for functional selectivity, in current theoretical evaluation, we consider the intact operational model (Eq. 11) as the true model and then assess the influence of model misspecification (i.e., the marginal operational model) on the power of the ligand bias test via stochastic simulation estimation studies. In this circumstance we note that the intact operational model should not be inferior to the marginal operational model, rather, equivalent performance would favour the marginal operational model as a more parsimonious approach.

Power analysis process

The power analysis is divided into two steps. In the first step, we correct for false positive findings by calibrating the statistical test to achieve the appropriate type I error rate of 0.05. In the second step, the power was calculated for different levels of ligand bias with an appropriately calibrated type I error rate.

Calibration of the criterion for type I error (α = 0.05) via non-parametric approach

Pseudo-experimental data were simulated from the true model (i.e., the intact operational model) under the null hypothesis (H0:ΔΔlogR1−2 = 0), based on the given parameter values and study design in Table 1. The parameter values for Table 1 were chosen to represent a biologically possible system[5]. Here we simulated 1000 virtual experiments and estimated the parameters in each. In each experiment, twelve concentration-response curves corresponding to a pair of reference and test ligands in two signalling pathways (three replicates for every situation) were generated with random error. Then, the dataset generated from every simulation experiment was estimated with the candidate model. For each simulation experiment, we calculated the Wald statistic (the ratio of squared ΔΔlogR1−2 estimate and its variance, ). The 95th percentile from 1000 Wald statistics was used as the cut-off criterion for type I error (α) equal to 0.05.

Table 1

Parameter values and study design for the calibration of alpha criterion (α = 0.05).

Parameters	Reference	Test
Design
logA	From −11 to −4, increase by 1
Pathway 1
E m1	100	100
Basal 1	10	10
logR 1	7	6
Pathway 2
E m2	100	100
Basal 2	10	10
logR 2	7	6
logKA′	−5	−6
ΔΔlogR1−2	0
Random error
add.err1	5	5
prop.err1	10%	10%
add.err2	5	5
prop.err2	10%	10%

Calculation of power curve

The pseudo-experimental data were simulated from the true model (i.e., the intact operational model) under alternative hypothesis (H0: ΔΔlogR1−2 ≠ 0), based on the given parameter values and study design in Table 2 (the same as Table 1, except that ΔΔlogR1−2 was different from zero). For each ΔΔlogR1−2 value, we conducted 1000 simulation and estimation experiments. In each experiment, twelve concentration-response curves corresponding to a pair of reference and test ligands in two signalling pathways (three replicates for every situation) were generated with random error. Then, the dataset generated from every simulation experiment was estimated with both candidate models (the marginal and intact operational models). For each simulation experiment, we calculated Wald statistic . For each ΔΔlogR1−2 value, the power was calculated by taking the percentage of Wald statistics that were greater than the criterion for type I error equal to 0.05. Since, in this experiment, ligand bias was present then the proportion of statistically significant runs is an empirical approximation to the true power.

Table 2

Parameter values and study design for the calculation of power curve.

Parameters	Reference	Test
Design
logA	From −11 to −4, increase by 1
Pathway 1
E m1	100	100
Basal 1	10	10
logR 1	7	6.1, 6.2, 6.3, 6.4, 6.5, 6.6, 6.8, 7.0
Pathway 2
E m2	100	100
Basal 2	10	10
logR 2	7	6
logKA′	−5	−6
ΔΔlogR1−2	0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.8, 1.0
Random error
add.err1	5	5
prop.err1	10%	10%
add.err2	5	5
prop.err2	10%	10%

For power analysis of the intact operational model, the intact operational model (Eq. 11) was the candidate model and and was directly obtained from the output of NONMEM. For power analysis of the marginal operational model, the marginal operational model (Eq. 1) was the candidate model and and was calculated through post hoc analysis (Eqs 2 and 3).

Results

The cut-off criterion for type I error equal to 0.05 was 5.53 for the intact operational model and 3.67 for the marginal operational model. The theoretical criterion was 3.84. For intact operational model, there was a slight increase in the type I error rate. These cut-off criterion values were used in the power analysis for controlling of type I error. The power analysis results are presented in Fig. 2. It is seen that the intact operational model had greater power to detect ligand bias compared to the marginal operational model for all values of ΔΔlogR1−2. It was shown that ΔΔlogR1−2 had to be greater than 0.8 for the marginal operational model to achieve 80% power for the ligand bias, while, for the intact operational model, the requirement of ΔΔlogR1−2 was approximately 0.4.

Figure 2

The power curve for the marginal operational model and the intact operational model. The orange line indicates the power curve for the marginal operational model and the blue line indicates the power curve for the intact operational model. The red dashed line indicates 80% power.

Part III

Application of the intact operational model

The performance of the intact operational model was evaluated using two examples from the literature. In both examples signalling responses were measured, simultaneously, a desirable but not necessary characteristic for application of the intact model (to reduce assay-specific effects). In the first example, calculation of logR for the marginal operational model does not support continued experiments since logR could not be precisely determined for the single ligand under the marginal operation model, whereas estimation was sufficiently precise under the intact model to support a full analysis. In the second example, we work through to calculation of ΔΔLogR and arrive at a refined (and different) conclusion from the intact operational model.

Example I: the opposing effects of Gi and Gs signalling

Data description: In example I, data were extracted from the literature that described the opposing effects of Gi and Gs signalling pathways on α2-adrenergic receptor subtype C10 mediated adenylate cyclase activity in transfected Chinese hamster ovary (CHO) cells expressing different levels of receptors[14]. At low expression levels (1 pmol/mg), only the inhibitory effect of the Gi pathway was observed; at higher receptor levels (10 pmol/mg), an atypical biphasic response was observed due to stimulatory effect of Gs pathway.

Model analysis: Marginal operational model: In this example the two signalling pathways (Gi and Gs) resulted in a change in the same pharmacological response (cAMP). The marginal operational model was implemented as the sum of the general operational model (Eq. 1) for each pathway (Eq. 12 for low receptor expression level and Eq. 13 for high receptor expression level). With the assumption that each signalling pathway was independent, the marginal operational model was implemented without any constraint on the functional affinity in different pathways. Hence, K could be different in Gi and Gs pathways, denoted as K and K.

The subscript ‘i’ and ‘s’ indicates the Gi and Gs pathways, respectively. The superscript ‘L’ and ‘H’ indicates the low receptor expression level and high receptor expression level, respectively. We note here that the basal level of activity is common to both pathways.

By definition of the transduction coefficient (R) and the transducer ratio (τ), it was demonstrated that the transduction coefficient was proportional to the total receptor density (R) (Eq. 14).

Hence, the relationship between transduction coefficients in the Gi pathway at two receptor expression levels (1 pmol/mg and 10 pmol/mg) could be established (Eqs 15 and 16).

Similarly, the relationship between transduction coefficients in the Gs pathway at two receptor expression levels (1 pmol/mg and 10 pmol/mg) was also established (Eq. 17).

Then, substituting Eqs 16 and 17 into Eq. 13 yields the equation for quantifying pharmacological response at high receptor level (Eq. 18).

In this case, Eqs 12 and 18 were used to model the data from Example I, where; logK, logK, and were the drug specific parameters that would be directly estimated.

A preference of a ligand towards a particular pathway before normalisation to a reference ligand was termed the ligand preference profile. It could be regarded as an intermediate metric from the first normalisation process for calculating ligand bias metric (detailed in Appendix 3). Hence, the ligand preference profile between Gi and Gs signalling pathways was defined as Eq. 19:

For the marginal operational model, the ligand preference profile (logR) was calculated using Eq. 19 and the estimated standard error (SE) was calculated using Eq. 20:

Intact operational model: The intact operational model was implemented as Eq. 21 for low receptor level and Eq. 22 for high receptor level. Note here that the apparent equilibrium dissociation constant () was shared by Gi and Gs pathways.

Based on the definition of ligand preference between Gi and Gs signalling pathways (Eq. 19), are expressed as the difference of and logR (Eq. 23):

Substituting Eq. 23 into Eq. 17 yielded the expression of (Eq. 24):

Then, substituting Eqs 16, 23 and 24 into Eqs 21 and 22 yielded the intact operational model with logR as a directly estimated parameter (Eq. 25 for low receptor level and Eq. 26 for high receptor level).

In this case, , and logR were the ligand specific parameters. Hence, the profile metric (logR) is estimated directly and the estimated standard error can be obtained from the NONMEM modelling output.

As illustrated in Fig. 3, the prediction from the intact operational model could reproduce the effect of the α2-adrenoceptor agonist UK-14304 on adenylate cyclase in transfected CHO cells, both qualitatively and quantitatively. This was confirmed by the model evaluation plots (Figures S1 and S2).

Figure 3

The fitting result of the intact operational model for the effect of the α2-adrenoceptor agonist UK-14304 on adenylate cyclase in transfected CHO cells expressing α2-C10 receptors at 1.0 and 10.0 pmol/mg. Note, Gi and Gs pathways are acting under both conditions, however the Gi pathway dominates at the low receptor expression and Gs pathway can eventually overcome this at high receptor expression. Left: at expression level of 1.0 pmol/mg, Gi-protein mediated inhibition of adenylate cyclase played a dominant role. Right: at expression level of 10 pmol/mg, a ‘U-shape’ response was observed due to activation of Gs-protein. The red dots are the data grabbed from the literature[14]. The black lines are the simulation profiles from the intact operational model.

From Table 3, the estimated parameter values were generally comparable between the intact operational model and the marginal operational model. However, it was evident that the intact operational model rendered more precise parameter estimates than the marginal operational model. Some parameters (e.g., E and n2) from the marginal operational model yielded poorly estimated parameters with relative standard errors of more than 200%, indicating that the data (from a standard pharmacological experiment) was not sufficiently informative to allow precise estimation of these parameters. In particular, the ligand preference profile (logR = 2.22), calculated via a post hoc analysis from the marginal operational model, was imprecise with relative standard error of 69%. Given this level of imprecision there would be little value in continuing with further comparator ligands since discrimination between ligands is unlikely to be fruitful. Hence, by testing only one ligand and analysing it in this way it is possible to make a no-go decision until an improved experimental set up could be created. In contrast, with the same profile directly estimated from the intact operational model was much more precise with the relative standard error of 6.7%. Hence, it remains possible to detect ligand bias in upcoming experiments and further refinement of the experiment would be unnecessary if the intact operational model is used.

Table 3

Comparison of the estimation results from the intact operational model and the marginal operational model for the effect of α2-adrenoceptor agonist UK-14304 on CHO cells with different expression levels of α2C10 receptor.

Parameters	Intact operational model	Marginal operational model
	Estimated [RSE%]
Basal	103 [2.9%]	103 [2.9%]
E m, i	49.2 [8.8%]	48.7 [16.2%]
n i	1.47 [48.2%]	1.44 [52.0%]
E m, s	59.2 [23.3%]	69.9 [212%]
n s	2.27 [28.7%]	2.54 [302%]
logKA′	−5.32 [3.9%]	—
logK A, i	—	0 FIXa
logK A, s	—	−5.54 [57.6%]
logRiL	—	6.90 [1.7%]
logRsL	—	4.68 [32.5%]
logR i: s	2.28 [6.7%]	2.22 [68.7%]b

RSE, relative standard error.

A subscript letter ‘i’ indicates the Gi signalling pathway, ‘s’ indicates the Gs signalling pathway and ‘L’ indicates the lower receptor expression.

a: α2-adrenoceptor agonist UK-14304 appears to be full agonist of Gi pathways in both receptor expression levels (1 pmol/mg and 10 pmol/mg). Due to a deterministic identifiability issue of the marginal operational model, it is not possible to estimate logK values for full agonists from direct fitting to concentration-response curve. Conventionally, logK is arbitrarily fixed to 0 to solve this problem[3].

b: calculated from post hoc analysis.

Example II: simultaneous measurement of IP accumulation and AA release

Data description: Data from the simultaneous measurement of phospholipase C (PLC)-mediated inositol phosphate (IP) accumulation and phospholipase A2 (PLA2)-mediated arachidonic acid (AA) release after the activation of 5-HT2C receptors in the CHO-1C19 cell were extracted from[15].

Model analysis: For the marginal operational model, Eq. 1 was separately applied to each signalling pathway to fit the data. For the intact operational model, Eq. 11 was implemented to jointly model all the functional assay data. In order to circumvent the deterministic identifiability issue in marginal operational model, we followed the convention to arbitrarily set and [11].

The high efficacy ligand in both pathways, bufotenin (BUF), was chosen as the reference ligand[7]. For the marginal operational model, the ligand bias metric (ΔΔlogR1−2) and the estimated standard error were calculated via post hoc analysis (Eqs 2 and 3). For the intact operational model, the ligand bias metric was directly estimated and the estimated standard error could be obtained from NONMEM modelling output.

Test for ligand bias: The null hypothesis for this test was that ΔΔlogR1−2 was equal to 0. Since there were four comparisons in this example, the α value was adjusted to 0.0125.

As shown in Fig. 4, the intact operational model captured the signalling profiles of the 5-HT2C receptor agonists on PLA2-mediated AA release and PLC-mediated IP accumulation in CHO-1C19 cells. This was confirmed by the model evaluation plots (Figs S3 and S4).

Figure 4

The fitting result of the intact operational model for the effect of the 5-HT2C receptor agonists on PLA2-mediated AA release (left) and PLC-mediated IP accumulation (right) in CHO-1C19 cells. All the responses have been normalized to 5-HT maximal response. The dots are the data grabbed from the literature[15]. The lines are the simulation profiles from the intact operational model. Here, different colours indicate different types of ligands. Red: Bufotenin. Blue: DOI, (±)-1-(2,5-dimethoxy-4-iodophenyl)-2-aminopropane. Yellow: LSD, lysergic acid diethylamide. Green: Quipazine. Purple: TFMPP, 3-trifluoromethylphenyl-piperazine.

From Table 4, the parameter estimates from the intact operational model were comparable to the marginal operational model. Since the logR could be precisely estimated in the marginal operational model, considered as a go from a go/no-go decision, we moved forward to the full ligand bias analysis. It was noted that the ligand bias metrics calculated via post hoc analyses from the marginal operational model were generally less precise than those directly estimated from the intact operational model. For instance, the relative standard error of ligand bias for QUI from the marginal operational model was almost 400%, much higher than that from the intact operational model (41%). For the test of ligand bias (Table 5), the results of two ligands (QUI, and TFMPP) were comparable from both methods and no statistically significant bias would be concluded for both QUI and TFMPP. However, the intact operational model was able to identify two biased ligands (DOI and LSD) that were not classified as biased ligands by the marginal operational model.

Table 4

Comparison of the estimation results from the intact and marginal operational models for the effect of the 5-HT2C receptor agonists on PLA2-mediated AA release and PLC-mediated IP accumulation in CHO-1C19 cells with bufotenin as the reference ligand.

Parameter	Intact operational model					Marginal operational model
	Estimated [RSE%]
Basal 1	1.79 [21.2%]					1.83 [22.6%]
E m1	102 [3.7%]					98 [3.3%]
n 1	1.53 [5.0%]					1.55 [5.8%]
Basal 2	0 FIXa					0 FIXa
E m2	126 [9.6%]					121 [5.6%]
n 2	0.80 [6.2%]					0.78 [5.5%]
	BUF	DOI	LSD	QUI	TFMPP	BUF	DOI	LSD	QUI	TFMPP
logKA′	−6.05 [1.2%]	−6.46 [1.4%]	−7.77 [0.9%]	−5.24 [2.7%]	−5.97 [3.2%]	—	—	—	—	—
logK A1	—	—	—	—	—	−6.02 [1.4%]	−6.52 [1.0%]	−7.81 [0.9%]	−4.72 [12.1%]	0 FIXb
logK A2	—	—	—	—	—	0 FIXb	−5.98 [2.8%]	−7.66 [1.3%]	−5.45 [1.3%]	−5.86 [2.8%]
logR 1	6.50 [0.5%]	6.65 [0.7%]	7.31 [1.0%]	5.89 [1.0%]	6.37 [0.4%]	6.51 [0.5%]	6.70 [0.6%]	7.37 [1.1%]	5.85 [0.8%]	6.39 [0.4%]
logR 2	6.11 [2.2%]	—	—	—	—	6.14 [1.8%]	6.65 [1.2%]	7.28 [1.6%]	5.44 [1.8%]	5.85 [2.2%]
ΔΔlogR1−2	—	−0.476 [14.3%]	−0.393 [14.4%]	0.167 [40.7%]	0.102 [131%]	—	−0.32 [45.6%]c	−0.28 [65.3%]c	−0.04 [396%]c	−0.17 [103%]c

RSE, relative standard error; CV, coefficient of variation; Add. error, additive error; BUF, bufotenin; DOI, (±)-1-(2,5-dimethoxy-4-iodophenyl)-2-aminopropane; LSD, lysergic acid diethylamide; QUI, quipazine; TFMPP, 3-trifluoromethylphenyl-piperazine.

A subscript number “1” indicates the IP accumulation pathway and “2” indicates the AA release pathway.

a: The estimated basal effect of AA release pathway was not significantly different from 0. Therefore, it was fixed to 0 in the model refinement process to make the model more stable.

b: TFMPP appears to be full agonist of IP accumulation pathway and bufotenin appears to be full agonist of AA release pathway. Due to the identifiability issue of the marginal operational model, it is not possible to estimate logK values for full agonists from direct fitting to concentration-response curve from that pathway. Conventionally, logK is arbitrarily fixed to 0 to solve this problem[3].

c: calculated from post hoc analyses

Table 5

Test of ligand bias of 5-HT2C receptor agonists on PLA2-mediated AA release and PLC-mediated IP accumulation in CHO-1C19 cells with bufotenin as the reference.

ΔΔlogR1−2	Intact operational model		Marginal operational model
	Estimated [RSE%]	p-value	Estimated [RSE%]	p-value
DOI	−0.476 [14.3%]	<0.0125	NS	0.028
LSD	−0.393 [14.4%]	<0.0125	NS	0.126
QUI	NS	0.014	NS	0.800
TFMPP	NS	0.447	NS	0.330

DOI, (±)-1-(2,5-dimethoxy-4-iodophenyl)-2-aminopropane; LSD, lysergic acid diethylamide; QUI, quipazine; TFMPP, 3-trifluoromethylphenyl-piperazine.

A subscript number “1” indicates the IP accumulation pathway and “2” indicates the AA release pathway.

NS indicates no statistically significant bias.

Discussion

In this work, an intact version of the operational model was proposed as a more complete representation of the biological system. The intact operational model follows the concept of functional selectivity by explicitly implementing linked equilibria among different receptor conformational states. Hence, the apparent dissociation constant () is the same among different signalling pathways to reflect the natural correlation of C50 among different pathways (Figure S5). The intact operational model allows for joint modelling of all data from different pathways and provides a direct estimate of ΔΔlogR, which avoids propagation of errors in the post hoc analysis. At a practical level this means the bias parameters (i.e. ΔΔlogR) are more precisely estimated for any given experimental design and the intact operational model has greater power (more sensitivity) to detect weak ligand bias that might otherwise be missed by the marginal operational model. This latter feature provides the potential for more in-depth inferences from experiments.

According to the power analysis evaluation, both marginal and intact operational model work well for highly biased ligands (ΔΔlogR > 1), with the power for ligand bias test exceeding 80%. However, only the intact operational model is able to detect less strongly biased ligands (ΔΔlogR < 0.8) with confidence (power exceeding 80%). The result from this theoretical evaluation is consistent with the finding from one previous study[16]. It is apparent that the marginal operational model fails to distinguish true ligand bias from ‘error’ with confidence, especially when the ligand bias is weak. Our work further supports this finding and provides statistical and mechanistic support for the intact model in detecting ligand bias.

The concept of approaching a pharmacology experiment within a go/no-go workflow framework was demonstrated in current work. In Example I, it would be plausible to consider a no-go decision based on the fact that logR could not be precisely determined for the single ligand when using the marginal operational model and hence the future ability to demonstrate ligand bias would therefore be unlikely. In this case, it would be necessary for the investigator to reconsider the experimental conditions. However, application of the intact operational model to the same data yielded a precise estimate of the ligand preference profile logR and hence would be anticipated to support further work to estimate full ligand bias (ΔΔlogR) from the current experimental set up. In Example II, the values of logR could be precisely estimated with the marginal operational model and hence a logical go decision would be made. We see however, in this example that even working through the full experiment the marginal operational model did not identify bias in some pathways when the intact operational model did.

Though the intact operational model is a more mechanism-based model framework for functional selectivity, it is not as flexible as the marginal operational model in curve fitting. In the operational model, for low efficacy agonists (i.e., as τ approaches 0), the C50 approximates K (Eq. 27). As the values of K are the same for different signalling pathways (all equal to ) in the intact operational model, this precludes fitting concentration–response curves of a ligand behaving as partial agonists in two pathways with very different C50 values which would require equally distinct K[17]. In addition, an exact numerical evaluation (detailed in Appendix 4) was performed in a case with a non-unity slope factor, showing that the difference between C50 values could not exceed (approximately) 100 fold. Such a change in C50 values is unlikely to occur biologically, but rather, when this is observed it likely reflects that the responses of different signalling pathways are commonly measured under different experimental conditions, which may break the natural linkage among different signalling pathways.

Another model for quantifying functional selectivity is Rajagopal’s model[18-20]. In Rajagopal’s model, the operational model is separately applied to each signalling pathway with the constraint that K among different pathways are all fixed to a previously estimated value from a separate binding assay. This model shares the same mathematical form with the intact operational model (detailed in Appendix 5). Hence, the intact operational model could be regarded as one possible biological interpretation for Rajagopal’s model. However, it is noted that in Rajagopal’s model, this shared functional affinity is fixed rather than determined by the data[18]. The intact operational model does not require this constraint. Whether functional affinity should be fixed to a measured value from another binding assay or estimated as part of the current experimental set up is out of the scope of current work. But it is worth noting that most binding assays are performed on cell membranes, while signalling assays are performed on whole cells and therefore may not be representative of each other. One possible strategy would be to use a fully Bayesian approach with informative priors on parameters that would have otherwise been fixed. More theoretical and practical evaluations are warranted to address this.

The three-state model is another mechanism-based model for quantifying function selectivity[8,15]. Similar to the intact operational model, the linked equilibria among different receptor conformation states and the mutual depletion of these receptor states are explicitly implemented in the three-state model. Within this model framework, the EC50 values are the same for different signalling pathways, which makes it less flexible than the intact operational model. Moreover, the three-state model only accounts for the receptor binding and does not include signal transduction processes. Hence, it cannot explain the phenomenon often termed ‘receptor reserve’, and no ligand may behave as a full agonist in both pathways. Additionally, given only functional assay data, the three-state model is not structurally identifiable and only a simulation-based, heuristic search for parameter values is possible (see for example[15]), which limits its utility in quantifying functional selectivity. Contrary to the empirical basal effect in the intact operational model, the three-state model mechanistically incorporates constitutive activity. In the future, it is desirable to combine these two model structures to render more insight into the underlying mechanism of functional selectivity.

As demonstrated in this work, the intact operational model works well when different signalling pathways are measured simultaneously. However, in practice, it is not uncommon that the responses from different signalling pathways are measured under quite different experimental conditions. The experimental conditions for the functional assay to study one pathway may eliminate the responses from other pathways. Under this circumstance, the association between different active receptor states can be disentangled and the intact operational model can be reduced into the marginal operational model (Appendix 6). In this sense, the marginal operational model can be regarded as a special case of the intact operational model when the responses from different pathways are measured under different conditions. However, in order to gain the maximum mechanistic insight into functional selectivity by application of the intact operational model, it is desirable that the links between different signalling pathways are preserved by keeping experimental conditions as consistent as possible between assay types.

For the accessibility of the intact operational model, our primary goal was to implement it in GraphPad Prism (a standard software in the area of experimental pharmacology). However, due to the nature of the model this was not possible within Prism. In this work, we used a more general and flexible analysis platform, in this case NONMEM. Though NONMEM has been widely used in the area of pharmacometrics[21], there are only few applications in analytical pharmacology[22]. There are other platforms like NONMEM offering flexibility in modelling, such as Monolix (developed by Lixoft) and Phoenix NLME (developed by Cetara), which could be used for analysis of this sort of data. (We provide the NONMEM estimation control streams as a Supplement).

In conclusion, the intact operational model can be applied to circumstances where the marginal operational model has been used. The intact operational model is more sensitive to identifying biased ligands (i.e. has greater power), and provides a more precise estimate of the operational model ligand bias metric (ΔΔlogR). The intact operational model may provide a valuable step to describe and improve understanding of the underlying mechanisms of functional selectivity.

Appendix

Supplemental Information