Considerations for improved performance of competition association assays analysed with the Motulsky–Mahan's “kinetics of competitive binding” model

Background and Purpose Target engagement dynamics can influence drugs' pharmacological effects. Kinetic parameters for drug:target interactions are often quantified by evaluating competition association experiments—measuring simultaneous protein binding of labelled tracers and unlabelled test compounds over time—with Motulsky–Mahan's “kinetics of competitive binding” model. Despite recent technical improvements, the current assay formats impose practical limitations to this approach. This study aims at the characterisation, understanding and prevention of these experimental constraints, and associated analytical challenges. Experimental Approach Monte Carlo simulations were used to run virtual kinetic and equilibrium tracer binding and competition experiments in both normal and perturbed assay conditions. Data were fitted to standard equations derived from the mass action law (including Motulsky–Mahan's) and to extended versions aiming to cope with frequently observed deviations of the canonical traces. Results were compared to assess the precision and accuracy of these models and identify experimental factors influencing their performance. Key Results Key factors influencing the precision and accuracy of the Motulsky–Mahan model are the interplay between compound dissociation rates, measurement time and interval frequency, tracer concentration and binding kinetics and the relative abundance of equilibrium complexes in vehicle controls. Experimental results produced recommendations for better design of tracer characterisation experiments and new strategies to deal with systematic signal decay. Conclusions and Implications Our data advances our comprehension of the Motulsky–Mahan kinetics of competitive binding models and provides experimental design recommendations, data analysis tools, and general guidelines for its practical application to in vitro pharmacology and drug screening.


| INTRODUCTION
The kinetics of drug interactions with their targets are thought to influence pharmacological effects in patients, especially in situations where binding equilibrium is not reached (Copeland, 2016;Copeland, Pompliano, & Meek, 2006;Paton, 1960;Paton, 1961;Swinney, 2006). Long target residence times can lead to occupancies outlasting the pharmacokinetics of a compound and, consequently, to sustained efficacy after the free drug has been cleared from the system. Moreover, differentiated kinetic profiles for primary and off-targets can result in kinetic selectivity despite similar steady-state affinities for both targets (Copeland et al., 2006;Dahl & Akerud, 2013). Thus, characterisation of screening hits, leads, and candidate compounds regarding their target binding kinetics (BK) has become a standard in drug discovery laboratories.
The growing interest in BK has promoted the development of assay technologies dedicated to measure association (k on ) and dissociation (k off ) rate constants, as well as binding affinities, (K D = k off /k on ; Antoine et al., 2016;de Witte et al., 2018;Guo, Mulder-Krieger, IJzerman, & Heitman, 2012;Schiele, Ayaz, & Fernández-Montalván, 2015;Stoddard et al., 2018;Stoddart et al., 2015;Swinney et al., 2014;Sykes et al., 2017;Xia, de Vries, IJzerman, & Heitman, 2016). Among them, competition association assays, monitoring the simultaneous binding of a labelled tracer (we have referred to the labelled compound as tracer, although any type of fluorescent or luminescent probes and radioactive tracers can be used for competition association assays) and an unlabelled compound to the target molecule, have gained popularity due to their relatively high throughput and ability to deal with both soluble and membrane-bound proteins (Schiele et al., 2015;Stoddart et al., 2015;Xia et al., 2016). In these assays, the rate constants for test compounds' interactions with their molecular targets are typically obtained from fitting the experimental data to the "kinetics of competitive binding" equation initially reported by Motulsky & Mahan (1984). This analytical procedure uses a differential equation system where the compound dose (I), tracer concentration (L), tracer k on (k 1 ), and tracer k off (k 2 ) values are known (fix) parameters, to solve for the unknown test compounds' k on (k 3 ) and k off (k 4 ) values. To this end, k 1 and k 2 are determined in direct binding experiments ideally performed under similar conditions as the competition assays.
In recent years, our laboratory routinely performed high throughput homogenous TR-FRET-based competition association assays (also known as kPCA: kinetic probe competition assay, Schiele et al., 2015) to evaluate hundreds of compounds interacting with several dozens of targets (see Bosma et al., 2019;de Witte et al., 2018;Georgi et al., 2018;Heroven et al., 2018;Nederpelt et al., 2016;Schiele et al., 2015). Typically, k on , k off , and derived affinities (K D,kin ) of test compounds showed excellent agreement among replicates, and with reference values. However, under our experimental conditions, we encountered two recurring situations in which analysing data with the Motulsky-Mahan model became challenging: The first issue (Case 1) was characterised by k on values being determined precisely in replicate measurements, but with a high variability of the corresponding k off and K D,kin (CV > 50%). In these examples, K D,kin always differ from the steady-state affinity values (K D,eq ) calculated from an independent equilibrium probe competition assay (ePCA). The second problem (Case 2) manifested itself by k on and k off parameters determined with poor precision (CV > 50%) despite acceptable data quality, and the fact that the corresponding K D,kin was consistent with K D,eq . While in Case 1, the k off was obtained by multiplying the on-rate by the K Deq from ePCAs (Schiele et al., 2015;Georgi et al., 2018) no meaningful BK values could be extracted from measurements described by Case 2 .
These are examples of decreased precision and accuracy and here it is necessary to define what we mean by these terms: Precision is the closeness of a set of measured replicate values to each other, and accuracy is the agreement of the measured value compared to the true value (JCGM200:2012 International vocabulary of metrology -Basic and general concepts and associated terms (VIM).3 edn.). The true value is the value that would be obtained by a perfect measurement but is indeterminable in reality. Nonetheless, the mean values acquired under repeatability conditions (=same assay components, measuring, and evaluation procedure) will be similar to the true value as long as there is no systematic error. • This study enables more precise and accurate kinetic parameter quantification for existing and novel medicines.
GEORGI ET AL. of competitive binding" model, to ask whether they were random or reproducible events, linked to specific experimental conditions. To address this question, we generated and evaluated thousands of simulated data points in a variety of assay set-ups via Monte Carlo (MC) simulations. The results described here provide explanations for the limitations described above, as well as recommendations for better design and analysis of these experiments, especially if they are run using the standard instrumentation of current drug screening laboratories.

| MC analyses
To test the performance of the non-linear multiparametric Motulsky-Mahan model within the parameter space, we performed MC simulations (Metropolis & Ulam, 1949) of thousands of pseudo-interactions with known input parameters, including the association and dissociation rates of the compounds, and with a random scatter in the measured signals using GraphPad Prism versions 6.07 and 7.00 for Windows (GraphPad Software, San Diego California USA, www.graphpad, RRID:SCR_002798). Our kPCA set-up (described in Schiele et al., 2015) served as reference for the experimental parameters, and-whenever applicable-deviations from the standard procedure were introduced as indicated in Table S1 to simulate systematic or random errors. These in silico experiments were subsequently evaluated by global curve fitting of the obtained signal traces to Motulsky-Mahan's "kinetics of competitive binding" model (Motulsky & Mahan, 1984) or others, depending on the scope of the simulation. For assessment of accuracy, we calculated the relative error of the output variables compared to the true value (input variable). For evaluation of precision, we calculated the coefficient of variation (CV) for all data obtained under repeatability conditions. Details of each pseudo-experiment are specified in Table S1, as well as in the corresponding figures. In addition, example GraphPad Prism files with the simulations are provided as Data S4 and Data S5 in the Supporting Information file.

| Models and equations
Most mathematical models used for our analyses are readily available in GraphPad Prism, and some were generated using the software's equation editing tool. Induced fit and irreversible interactions were modelled with COPASI 4.1.9 (Hoops et al., 2006), RRID:SCR_ 014260 and are provided as Data S2 and Data S3 in the Supporting Information file. Previously described models and their literature sources are described in the supplementary methods section. Models generated for this study are described below.
Normalised Motulsky-Mahan "kinetics of competitive binding" The percentage of equilibrium reached was approximated as follows:

| Data normalisation and statistics
For the normalised "Motulsky-Mahan" model, kinetic traces were normalised as percentage of tracer binding (0%: no tracer binding = background signal = maximum inhibition control; 100%: tracer binding control = no compound binding) and inverted to obtain percentage of compound binding.
Precision was quantified by determination of the CV, which is the quotient of the SD and the mean: CV [%] = SD/mean × 100%.
Accuracy was quantified as relative error of the determined mean value compared to the input value for the simulations: Z'-factor was calculated as described below (Zhang, Chung, & Oldenburg, 1999). The mean of all Z'-factors at each time point at equilibrium of the BK assay was used to approximate the Z'-factor of the assay.
where SD p , SD n , mean p , mean n : SD and mean of positive p (here: 100% tracer binding) and negative n (background signal, here: 0% tracer binding) controls.
The binding parameters obtained by the two 'signal decay' methods described in the supplementary method section (multiplication with signal drift term or normalisation of the "kinetics of competitive binding" equation) were compared to literature radioligand binding data by (a) Spearman correlation calculations (coefficient r and two-tailed P value) and (b) Bland-Altman analysis (Bland & Altman, 1986, 1990; mean log differences) using GraphPad Prism.

| Nomenclature of targets and ligands
Key protein targets and ligands in this article are hyperlinked to corresponding entries in http://www.guidetopharmacology.org, the common portal for data from the IUPHAR/BPS Guide to PHARMACOLOGY (Harding et al., 2018),

| Precision and accuracy of rate constant determination is significantly influenced by target residence times of test compounds'
The Case 1 described in the introduction prompted us to ask if onrates are precisely fitted when the corresponding off-rates and affinities cannot be determined; the underlying question being whether it is justified to calculate off-rates as proposed in previous work Schiele et al., 2015). Likewise, to better understand if and to what extent Case 2 results can be trusted, we decided to explore the range of rate constants that can be quantified with acceptable precision and accuracy under our experimental conditions for competition association assays.
The results of the corresponding MC analyses (described in Table S1, Experiment 1) are shown in Figure  combinations of 10 4 -0.1 and 10 9 -0.01 respectively. Further simulations showed that this particular problem is considerably alleviated by using at least 10-fold lower or higher doses of test compounds (not shown).
Figure 1a also shows that there is an excellent overlap of measured and fitted competitive association traces for compounds sharing the same on-rate and with dissociation rates of 0.0001 s −1 or slower. Consequently, the on-rates obtained are consistent with the input data, with relative errors of <4% and CV < 7%, but the affinities and off-rates are calculated with a significant lack of precision and accuracy (errors and CVs of >>500% and >75% respectively, Figure  These results suggest that on-rates in this situation are determined with high precision and accuracy and that the reason why off-rates and affinities cannot be determined is most likely related to the long residence times of the test compounds. In the following sections, we will use the term "slow off-rate problem" to describe this in silico reconstitution of real life's Case 1.
Compounds placed along an isoaffinity diagonal also show agreement in simulated and fitted competitive association traces when dissociation rates are faster than 0.1 s −1 . Accordingly, the output affinities are consistent with the true value (relative errors <0.7% and CV < 9%); however, the on-and off-rates drift to extremely high values (Figure 1b), reflected in relative errors of >>500% and CV > 298% respectively. This simulation matches the experimental issue depicted in Case 2 and will be referred to as the "fast off-rate problem." 3.2 | Accurate off-rate calculations from on-rates and steady-state affinities depend on sufficiently fast tracer BK and appropriate equilibration times The results described above suggest that in the Case 1 situation, onrates together with affinities from independent equilibrium endpoint experiments can indeed be used to calculate the off-rate. To determine the minimum incubation time required for precise and accurate K D,eq determination, we simulated our standard equilibrium probe competition assay procedure (Schiele et al., 2015) at equilibration times ranging between 1 and 24 hr and assessed the affinities of a 12-compound subset from the 35 described above, assuming similar tracer BK ( Figure S1 and characterising the influence of assay parameters such as the total measurement (observation) time, the sampling frequency (measuring interval), the tracer concentration, and its BK parameters in the precision and accuracy of parameter determination by the kinetics of competitive binding model. To this end, we considered hypothetical compounds for which either the "slow off-rate" or the "fast off- The rate plots represent the input binding kinetic parameters (left panel, and large symbols in all plots) for the 35 compounds (one colour per compound) as used for the simulation along with the corresponding output rates (small symbols in middle and right graphs) calculated by using the Motulsky-Mahan model. The diagonals in the plots correspond to the isoaffinity lines. The right panel zooms into the range of the input parameters, which span a physiologically meaningful range of the binding parameters. In contrast, the middle panel allows visualisation of all generated output parameters rate problem" occurs (i.e., k on [M −1 s −1 ], k off [s −1 ]: 10 6 , 10 −4 or 10 6 , 10 −1 respectively) and performed further MC simulations as described in Table S1, Experiments 2 to 5.
The results of these experiments show that increasing the observation time from 100 to 3,600 s decreased stepwise the relative error and CV for k off and K D determinations for compounds with the "slow off-rate problem" (k off ≤ 0.0001 s −1 ) from >500% and around 150% to <4% and <22% (Figures 2, S3A, and S4A). In contrast, compounds exhibiting first signs of the "fast off-rate problem" were less affected by the observation time, suggesting that the total observation time determines the longest residence time quantifiable in "kinetics of competitive binding" experiments and slower off-rates can be reliably determined by using longer incubation times.
Decreasing the measuring interval from 100 to 1 s progressively reduced the relative error and CV for k on and k off determinations for "fast off-rate problem" compounds (k off ≥ 0.1 s −1 ) from >500% and >400% to around 1% and <7% (Figures 3, S3B, and S4B). Interestingly, the relative error and CV for k off and K D of the slow-dissociating compounds (k off ≤ 0.0001 s −1 ) were also affected, but the slow off-rate problem could not be completely prevented, with CV changes from    Figures S6-S7). In agreement with our observation for slowly dissociating compounds, the on-rates for an irreversible binding compound could be also quantified with high accuracy ( Figure S7B). Interestingly, deviations from the true tracer concentration and BK, as well as compound concentration errors in a single well, accounted for the biggest accuracy issues ( Figure S7A), with different dependencies of effect, direction, and magnitude for fast-and slow-dissociating compounds.
3.5 | Simultaneous determination of association and dissociation kinetics in a single titration experiment favours precision and accuracy of tracer characterisation Tracer binding parameters are commonly determined in kinetic association experiments followed by data fitting to a pseudo-first-order exponential rate equation. Sometimes, a displacement (chase) experiment is performed separately, and the off-rate can be determined by fitting the curve to a single phase exponential decay model (Schiele et al., 2015). Recently, we described an alternative approach where binding of increasing tracer doses is followed by displacement by an excess of "cold" tracer within the same experiment, and the kinetic traces are globally fitted to an 'association-then-dissociation' model accounting for various tracer concentrations (de Witte et al., 2018;Nederpelt et al., 2016). In order to compare the precision and accuracy of these methods, we conducted MC simulations for a series of hypothetical tracers spanning a broad range of on-and off-rates (examples are shown in Figure 5a) and included a random fluctuation of the fluorescence signals. Data were then evaluated with the respective multivariable non-linear models (Table S1, Experiments 8 and 9) and fit performances were compared. Despite similar overall performance, the analysis suggests that the approach combining association and Monte Carlo simulations and analyses similarly or identical to those shown in (a) were performed 100 times respectively. The graphs represent the input binding kinetic parameters (horizontal and diagonal dashed lines) along with the output parameters (grey dots) calculated by using either the "association kinetics" or the "association-thendissociation" model. Not all output parameters are inside the y-axis limits dissociation is beneficial when dealing with fast-dissociating tracers ( Figure 5b). "Kinetic association assays" showed lower precision and accuracy for determination of slowly dissociating compounds (10 −5 s −1 ), an issue that can be prevented by longer incubation times (not shown). For tracers with low affinities (<10 −7 M) and fastdissociation rates, the corresponding kinetic rate constants could not be determined accurately. With the "kinetic association-thendissociation assay," all off-rates <1 s −1 (and the corresponding on-

FIGURE 6
Comparison of two kinetics of competitive binding models dealing with systematic signal decay. (a) The left graph shows an example of a competitive tracer-compound binding trace with systematic signal decay (from kPCA experiment, Bosma et al., 2019) and the corresponding fit (solid line) derived from the Motulsky-Mahan model equation multiplied with a signal drift term where K Drift = 0.0028 ± 0.0004. The signal is also decaying in the control well without compound. The right graph shows the associated normalised traces: The percentage of binding was calculated by normalisation between the tracer binding control (0% compound binding) and the background signal (100% compound binding). The corresponding fit (solid line) was calculated by using the normalised Motulsky-Mahan equation. (b-c) For nine compounds, the resulting binding parameters from both fits were compared to those obtained from radioligand binding experiments (Bosma et al., 2019) represents the correlation plots. Spearman correlation coefficients are indicated below the graphs. Panel (c) shows Bland-Altman plots. The solid horizontal lines indicate the mean log difference for all data points-which is also given as mean ± SD below the graphs. The dashed lines represent the upper and lower 95% limit of agreement rates) could be determined with high precision and accuracy, regardless of the affinity of the tracer, and performance issues were only observed for faster dissociating probes. Moreover, the associationthen-dissociation signal traces can be interpreted more intuitively, allowing for easy recognition of slow-and fast-dissociating tracers ( Figure 5a).
3.6 | Systematic fluorescence signal decay can be addressed both analytically and with raw data normalisation approaches Fluorescence-based competition association assays are sometimes affected by signal decay over time, which is not attributable to the pharmacological events being investigated (Bosma et al., 2019;Nederpelt et al., 2016;Schiele et al., 2015). We previously hypothesised this loss of fluorescence to be linked to photobleaching and introduced an equation in which the Motulsky-Mahan model was extended to account for this phenomenon (Schiele et al., 2015). With further expansion of our assay portfolio and expertise, we have noticed that fluorescence losses also depend on the nature of the target, the background, buffer and tracer warhead. Therefore, here, we developed a more general approach consisting on a modified Motulsky-Mahan model that can be used to fit normalised kinetic traces. Figure 6a shows that both correction procedures generate better fits to literature data for which comparative radioligand experiments with no signal drift were available (Bosma et al., 2019).
Parameters (k on , k off , and K D ), obtained with both models, are in good agreement with data from orthogonal assays (r > 0.92, mean log difference < 0.15, Figure 6b-c). Nevertheless, the correlation is slightly better for normalised fluorescence data (mean log difference: 0.10 ± 0.25 vs. 0.12 ± 0.35 for photobleaching model). Both methods are now available in the Genedata Screener® software, (Dubrovskiy et al., 2016) and can be easily implemented in comparable analysis software.

| DISCUSSION AND CONCLUSIONS
Differences in the clinical outcomes of drugs addressing a particular medical need are often defined by their molecular mode of action and-more specifically-by their target interaction kinetics (Swinney, 2004;Swinney, 2006). Not surprisingly, this parameter is nowadays monitored in drug discovery programmes all the way from the hit identification and lead optimisation phases up to the candidate evaluation and selection (Georgi et al., 2017). To this end, the Motulsky-Mahan "kinetics of competitive binding" model is widely applied, with recent efforts aiming at screening assay formats with increased robustness and throughput (Antoine et al., 2016;Guo et al., 2012;Schiele et al., 2015;Stoddard et al., 2018).
In this study, we have addressed challenges to the Motulsky- I. When the residence time of the compound is fivefold longer than the observation time, neither k 4 nor K D can be precisely fitted. In contrast, k 3 is determined with high precision and accuracy. While observation time can be increased to circumvent this problem; this solution decreases the throughput and will be often limited by the stability of assay components. Alternatively, off-rates can be calculated by multiplying the precise Motulsky-Mahan on-rate and the affinity determined in an independent equilibrium assay (k off = K D × k on ). In this context, attention should be paid to the fact that the calculated steady-state affinity can be influenced by the incubation time, as well as by ligand and target depletion (assay wall) issues (Aranyi, 1980;Easson & Stedman, 1936;Heise, Sullivan, & Crowe, 2007;Klebl, Müller, & Hamacher, 2011).
Looking for discrepancies when comparing K D,eq and K D,kin as recommended in the GraphPad Curve Fitting Guide (Motulsky, 1995 GraphPad Software, Inc.) is a helpful exercise to assess whether this approach is justified. Future developments of the model could include global fitting of competitive association (kPCA) traces and steady-state competition (ePCA) curves..
II. If k 4 is faster than the measurement frequency, neither k 3 nor k 4 can be correctly fitted, but the affinity is still determined with utmost precision and accuracy. This feature can be used as empirical tool to identify fast-dissociating compounds. Importantly, this implies that a good correlation of K D,eq and K D,kin is not a sufficient quality criterion for k 3 and k 4 determinations. According to the GEORGI ET AL.  (Bosma et al., 2019). In that paper, the choice of slow-and fast-dissociating tracers was suggested to be the main factor influencing experimental outcomes. Here, we provide an additional explanation to the differences observed, relying on the distinct measuring intervals and observation times chosen. Along these lines, here, we remind readers that k 1 and k 2 accuracies are critical for the quality of k 3 and k 4 determinations and demonstrate that experimental tracer characterisation approaches involving "association-then-dissociation of multiple tracer concentrations" can increase the accuracy of parameter determination. Also in relation to (Bosma et al., 2019), where the HTRF-based assay was affected by a systematic signal drift, here, we propose that normalisation of the Motulsky-Mahan model and signal traces is a valuable alternative to cope with this issue, especially in cases where photobleaching might not be the only source of fluorescence decay or the latter is not mono-exponential. Of note, photobleaching is often more complex than a mono-exponential process (Rigaut & Vassy, 1991;Song, Hennink, Young, & Tanke, 1995), and the signal decrease can be related to other phenomenon such as cell sedimentation or protein target degradation. Having said this, the downside of a normalised model with no parameter for B max is the potentiation of issues associated with noisy signal traces and fast-dissociating compounds.
In conclusion, this study provides new insights into the dependence of the Motusky-Mahan model on experimental parameters.
Experimentalists are encouraged to use the MC simulations presented here for detailed evaluation of their own procedures (adjustments may be necessary if they differ significantly from the kPCA assay set-up, Schiele et al., 2015) and to address the remaining open questions, such as the influence of ligand and target depletion. Along these lines, further studies could include the expansion of these analyses to recent variants of the "kinetics of competitive binding" equation that incorporate alternative scheduling of the labelled and unlabelled ligand addition (Hoare, 2018;Shimizu, Ogawa, & Nakayama, 2016) or two-state binding (Guo et al., 2018). All in all, this work should enable the optimised design and analysis of competition association assays.

ACKNOWLEDGEMENTS
This study was partly undertaken within the framework of the "Kinet-