Kinetic tools for the identification of ligand–receptor interaction mechanisms

Tools for the identification of ligand–receptor interaction mechanisms were developed by mathematical modeling of the influence of a ligand on the kinetics of a reporter ligand binding with a receptor. This approach allows kinetic differentiation between ligands of both rapid and slow binding modes, but also distinguishes compounds that share the same binding site with the reporter ligand or bind non-competitively to a distinct binding site. In order to simulate the kinetic behavior of this system, a mathematical model comprising ordinary differential equations was derived and solved numerically.


INTRODUCTION
Implementing various physical chemistry concepts to explain physiological phenomena has led to an understanding that the dissociation constant of a drug-receptor complex (K d ) is one of the most fundamental characteristics that describe the pharmacological properties of drugs, as well as the bioactivity of various toxins and biosimilars in general. Nowadays, this concept, initially based on simple drug-receptor interaction equilibria, has progressed in different directions, including the idea of a multi-step binding mechanism introduced by Strickland et al. [1]. Following this, the apparent effectiveness of the drugreceptor complex formation can be amplified through the additional drug-receptor complex isomerization step [1]. This way of modulating a drug's overall effectiveness is especially important if the isomerization step is a relatively slow process [2]. Moreover, this knowledge provides new possibilities for the advancement of drug design that focuses on the isomerization process. Recently, the importance of this kinetic aspect of the drug-receptor interaction mechanism has once again been highlighted by the introduction of the term "drug residence time" in the receptor complex, which represents metaphorically the slowness of the ligand-receptor complex dissociation rate [3].
All of these experimental studies were based on the determination of the interaction of a labeled reporter ligand with the receptor by filtrating samples through GF/B filters. This technique yields reproducible results if the reporter ligand A forms a slowly dissociating complex (RA) that can be reliably separated from the excess of radioligand (Scheme 1).

Scheme 1. Reporter ligand A binding to receptor R.
In this scheme, the first binding step is rapid, having association and dissociation rates with half-lives of the order of seconds, and is followed by a slow isomerization step, characterized by half-lives of the order of minutes. The rate equation for this reaction can be solved, and for excess of ligand A over the receptor concentration it yields a hyperbolic plot of the observed rate constant (k obs ) of complex (RA) formation versus ligand A concentration: where the rate constants k 12 and k 21 correspond to Scheme 1 and K A is the equilibrium dissociation constant of the complex RA, K A = k 10 /k 01 . In experiments, the constant k 21 represents the observable offrate of the reporter ligand [4]. Based on this reaction mechanism, a kinetic assay procedure was developed for characterizing the interaction of another non-labeled ligand B with receptor R using ligand A as a reporter ligand [11]. This method significantly extends our perspectives to help understand the mechanisms of the action of various bioactive compounds, which interact with the receptor protein, and certainly allows more detailed investigations into the binding modes of agonistic, partially agonistic, and antagonistic ligands.
Until recently, this kinetic analysis has rarely been used because of the complexity of kinetic models that describe the possibilities of simultaneous binding and isomerization steps for two ligands. In order to analyze the kinetic behavior of this complex assay system, a kinetic model for a rather general ligandreceptor interaction scheme (Scheme 2) was compiled, taking into consideration different modes of ligand binding and isomerization steps and the simulated influence of ligands A and B on the kinetics of the reporter ligand binding. This general reaction scheme analyzed in this study is presented in Scheme 2. To describe the kinetic behavior of this complex process, a system of ordinary differential equations was derived (see Section 2), and the numerical solution of these equations yielded a series of characteristic plots, demonstrating the dependence of the rate of the reporter ligand binding on the concentrations of ligands A and B. These plots can be used as tools for distinguishing different ligand binding modes.

MATHEMATICAL MODEL
The mathematical model for Scheme 2 is based on the following system of ordinary differential equations.
This system of ordinary differential equations is derived on the basis of the concepts of physical chemistry and the principles of mass conservation. Coefficients k ij describe the corresponding reaction rate constants in Scheme 2. Unknowns denoted by square brackets represent the molar concentrations of several complexes composed of receptor (R), reporter ligand (A), and ligand (B) depicted in Scheme 2. Initial values and parameters are presented in Appendix 1; constants describing the reporter ligand A were taken from the paper by Kukk and Järv [10]; the values of k 01 , k 03 , k 15 , and k 35 are diffusion-limited rate constants.
In most cases, a one-term exponential fit is sufficient to describe the time curves of binding; when a more complex situation is observed, a sum of two exponents should be used. Adding a third exponential term has little or no effect. Binding data were fitted with the MATLAB fit function with option exp2, which represents a sum of two exponentials: where [(RA)] is the concentration of the isomerized ligand-receptor complex bound at time t, exponent coefficient d is the apparent rate constant of the reaction, and c is the specific binding at equilibrium. The contribution of the first factor of Eq. (2) is usually very small, and is neglected in this paper. In cases where ternary (RAB) complex formation is possible, since both complexes are separated by slow methods, e.g. filtration. Equation (2) is an analog of a function describing ligand binding time curves used before: where B t is the binding at time t, B eq is the binding at equilibrium, and k obs represents the observed rate constant describing the rate of binding [4].
The most important outcomes of simulations are time curves of (RA) binding and the coefficients of exponential fit, which represent the observed rate constants and the specific binding at equilibrium. These curves are generated at different concentrations of the ligand, and the effects can be seen when plotting the inhibitor concentration against the observed rate constants. These plots reveal if the ligand has the ability to form a slowly dissociating step and what the on-and off-rate constants of the slowly dissociating step are, the dissociation constant of the fast step (K A ), and the inhibition constant of the other ligand. Relative rate coefficients were calculated by dividing the apparent rate coefficient by the maximum rate coefficient at extremely high concentrations of the reporter ligand.
MATLAB version 9.1.0 with Curve-Fitting Toolbox (2016b, Mathworks, USA) was used for the simulations. The ordinary differential equation system describing Scheme 2 was solved numerically with the MATLAB ode23s function designed to solve stiff differential equations. The interpolation between curves to yield a surface was achieved using the Delaunay triangulation.

Reporter ligand A
Interaction of a reporter ligand A with a receptor R must obey a specific kinetic mechanism, where the overall two-step binding process should include a fast complex formation followed by its slow isomerization (states 0, 1, 2; denoted as superscripts in Scheme 2). This two-step reaction scheme is a precondition for all receptor assays, where the ligand-receptor complex is separated from the excess of the reporter ligand by slow methods, e.g. filtration, and the concentration of the isomerized complex is determined.
The results of the simulations shown in Fig. 1 correspond to kinetic curves of the (RA) complex formation and were further analyzed using Eq. (2), which is the general function for the model used. Coefficient c is the specific binding of the reporter ligand, and coefficient d corresponds to the rate constant. Coefficients c and d depend on the concentration of ligand A and level off at high ligand concentrations, as shown in Fig. 2. Figure 2a corresponds to the common ligand binding curve used for the calculation of K d from the specific binding versus the concentration of A ([A]) plot. Data shown in Fig. 2b can be fitted by Eq. (1), which allows calculation of the apparent on-rate and off-rate rate constants for the slowly dissociating complex, as well as the K A value, where K A = k 10 /k 01 .  Taken together, the set of kinetic data shown in Figs. 1 and 2 can be used to identify ligands that induce the slow isomerization step of the receptor-ligand complex. In the rest of this study, a similar approach for more complex reaction mechanisms is adopted.

Competitive binding, rapidly dissociating ligand-receptor complex
Simultaneous interaction of reporter ligand A and another ligand B with the receptor R may follow competitive or non-competitive binding mechanisms. In this section, ligand B competes with ligand A for the same binding site but does not induce a slowly dissociating step. This situation is described by considering the formation of states 1, 2, and 3 (Scheme 2). For the rate constant k 30 , two different values (2 and 200 s −1 ) are used to mimic the binding effectiveness of B (Appendix 1, Fig. 1A). The results of this simulation are shown in Fig. 3, where coefficients c and d are plotted against the concentration of B, while the concentration of ligand A was fixed in both cases. It can be seen that the plot of coefficient c versus [B] in Fig. 3a yielded common displacement curves, shifted due to the different binding effectiveness of B. The marked IC 50 values should be corrected using the Cheng-Prusoff equation to obtain the K i values for B, which are equal to the ratio of rate constants k 30 /k 03 [12]. The plot of d versus [B] (Fig. 3b) reveals two descending hyperbolas, corresponding to the two models, where the binding effectiveness of B is different. This plot is characteristic for the competitive mechanism, where A and B share the same binding and B binds reversibly, while A induces the isomerization step that is responsible for the slow off-rate of the reporter ligand.
It can be concluded that the rate of the reporter ligand binding (coefficient d) is always decreasing in the presence of a rapidly dissociating B, which binds reversibly and competes with A for the binding site.

Competitive binding, slowly dissociating ligand-receptor complex
In this case, ligand B competes with the reporter ligand A and also induces the isomerized complex (RB). Therefore, the reaction scheme includes states 0, 1, 2, 3, and 4 (Scheme 2). Ligands of high and low binding effectiveness were characterized by different k 30 values (2 and 200 s −1 ). The kinetic curves of (RA) formation in the presence of ligand B are shown in Fig. 2A in Appendix 1. Analysis of these kinetic curves reveals that the rate constant of (RA) formation increases in the presence of B (Fig. 4b). This can be explained by the fact that R is consumed in two isomerization processes. At the same time the specific binding of ligand A decreases, demonstrating the competitive nature of A and B binding with the receptor. As a result of this competition, common displacement curves can be observed (Fig. 4a). It can be observed that the isomerization step, characterized by the ratio k 21 /k 12 , enhances the effectiveness of ligand binding.
The kinetic behavior of the overall process, simulated at various concentrations of A and B, is shown in Fig. 5b. This 3D plot allows characterization of the binding effectiveness of B (k 30 /k 03 ) and estimation of the rate constant k 34 of RB isomerization. It is important that these parameters, which characterize the interaction of B with R, can be obtained from the kinetics of the binding of A with R in the presence of B.
Contrary to what is shown in Fig. 5a, in the case of isomerization of the RB complex, the reaction rate of reporter ligand binding increases with increasing concentration of B (Fig. 5b), and this effect can be used to explicitly distinguish the two reaction mechanisms. Higher concentrations of A decrease the influence of ligand B, and the maximal value of the relative rate can be reached at high concentrations of ligands B or A.

Competitive binding, transition between rapidly and slowly dissociating ligand-receptor complexes
The two mechanisms discussed above present extreme cases, with clear differences between fast and slow binding modes of B. However, it is also interesting to consider intermediate situations between these two extremes. Therefore, kinetic curves were simulated with increasing k 43 values. The results of these calculations, where the concentration of A was fixed and the concentration of B was varied, are shown in Fig. 5c. It can be seen that the rate of reaction depends significantly on the ratio of the rate constants k 21 and k 43 , which describe de-isomerization of (RA) and (RB). At sufficiently high k 43 values, where deisomerization of (RB) is fast, the rate of (RA) formation slows down at high concentrations of B. The shape of this plot resembles the results shown in Fig. 5a, illustrating the rapid binding mode of B, while at low k 43 values the analogy with the results shown in Fig. 5b is obvious.

Non-competitive ligands A and B
Ligands A and B bind in a non-competitive manner if their binding sites do not overlap and formation of the ternary complex RAB can occur. In order to investigate the formation of this ternary complex, it is not sufficient to study the effect of ligand B on the kinetics of (RA) formation, but variation of both ligand concentrations is necessary. Three possibilities of non-competitive binding are considered:  RB and RAB do not isomerize (Scheme 2, complexes 1-3, 5);  RAB isomerizes (Scheme 2, complexes 1-3, 5, 6);  RB and RAB isomerize (Scheme 2, complexes 1-6).
In the first case, illustrated in Fig. 5d, increasing concentrations of ligand B always decrease the reaction rate, and this effect cannot be overcome by increasing the concentration of ligand A. The effect of B is more drastic than in the case of the competitive inhibitor (Fig. 5a), and reaction rate always decreases to a common plateau, independently of the concentration of A.
The second (non-competitive) case (Fig. 5e) is quite different from that shown in Fig. 5d: the relative rate is independent of the concentration of B, and the height of the plateau depends on the concentration of reporter ligand A. No ligand B concentration can increase the rate of (RA) formation. One should bear in mind that the (RAB) complex is also detected in this case. This means that at some point the specific binding might not decrease with increasing concentration of B, since at high concentrations of B (RAB) is formed instead of (RA).
Third, if both RB and RAB isomerize, the reaction accelerates with increasing concentration of either ligand A or ligand B (Fig. 5f). At any concentration of B, sufficiently high concentration of A will always accelerate the reaction, until the maximum relative rate is reached. This 3D plot is similar to competitive binding with RB isomerization (Fig. 5b), only differing in the steeper approach to the maximum relative rate. Meanwhile, the effect on the specific binding is less drastic than in the case of competitive binding because of the contribution of (RAB) to the specific binding.
In summary, these simulations demonstrate that all three options, which assume formation of the ternary complex RAB, can be analyzed and differentiated if the kinetics of the reporter ligand binding can be studied at a sufficiently wide range of concentrations of both ligand A and ligand B.

DISCUSSION
Previous attempts at analyzing ligand-receptor kinetics, which were made by solving systems of differential equations analytically [13] or by solving differential equations numerically [14,15], did not consider isomerization phenomena. As seen from the results of the present simulations, this extra step has a significant influence on the kinetic behavior of ligand-receptor interactions. In this study, a method for investigating complex mechanisms, summarized in Scheme 2, is proposed.
Ligands capable of forming a slowly dissociating isomerized ligand-receptor complex can be determined by using a constant reporter ligand concentration and monitoring the effect of ligand B on the rate coefficient. If the observed rate of reaction increases with increasing concentration of the unlabeled ligand, isomerization is expected; when the opposite is observed, only fast complexes are possible (Fig. 5a versus 5b). Formation of an isomerized ternary complex can be assessed by comparing the observed reaction rate constant dependence (Fig. 5d versus 5e) on the concentration of the ligand in question. Unless the isomerized ternary complex is very rapidly dissociating, the apparent reaction rate coefficient does not increase with increasing unlabeled ligand concentration, and a slowly dissociating isomerized ternary complex is highly probable. The reaction rate remains unchanged (Fig. 5e) or increases (Fig. 5f) if the ternary complex isomerizes with increasing concentration of B. The opposite is observed when the ternary complex does not isomerize (Fig. 5d).
When non-competitive binding of A and B is suspected, the concentration of the reporter ligand should be varied while keeping the concentration of the unlabeled ligand constant. This should be carried out at different concentrations of the unlabeled ligand, and the effect of the reporter ligand concentration on the apparent reaction rate should be monitored. If either ligand B or the ternary complex does not isomerize, non-competitively binding ligand B can be differentiated by the fact that different relative rate plateaus are reached depending on the concentration of B, and this cannot be overcome by high concentrations of ligand A (Fig. 5d). This is analogous to the distinction between different inhibition types in enzyme kinetics.
The simulations demonstrate that the rate of the ligand off-rate from its binding site can very effectively increase the apparent affinity of any ligand. The presence of this slow step can be identified by kinetic analysis, and the simulation results provide clear guidelines for it. In general, all mechanisms that involve isomerization of the receptor complex with ligand B exhibit an increase in the rate of the reporter ligand binding in the presence of B.

CONCLUSIONS
Solving a numerically ordinary differential equation system, derived for a complex ligand-receptor binding mechanism, provides a set of kinetic tools that can be used for differentiating distinct ligand-receptor complexes. These complexes are involved in rapid equilibria and in slow ligand-receptor complex isomerization steps that determine a fast or slow binding mode, sometimes referred to as short or long drug residence time. If a receptor interacts simultaneously with two ligands, this kinetic approach allows for differentiating between competitive and non-competitive binding mechanisms.

ACKNOWLEDGMENTS
This work was funded by Grants IUT20-15 and IUT20-57 of the Estonian Ministry of Education and Research. The publication costs of this article were covered by the Estonian Academy of Sciences. Competitive binding, rapidly dissociating ligand-receptor complex Fig. 1A. Effect of competitive ligand B concentration with higher (a) and lower (b) effectiveness of (RA) formation, rapidly dissociating ligand-receptor complex.  Competitive binding, slowly dissociating ligand-receptor complex Fig. 2A. Effect of competitive ligand B concentration with higher (a) and lower (b) effectiveness of (RA) formation, slowly dissociating ligand-receptor complex.