A rebinding-assay for measuring extreme kinetics using label-free biosensors

In vitro kinetic measurements allow mechanistic characterization of binding interactions and are particularly valuable throughout drug discovery, from confirmation of on-target binding in early discovery to fine-tuning of drug-binding properties in pre-clinical development. Early chemical matter often exhibits transient kinetics, which remain challenging to measure in a routine drug discovery setting. For example, characterization of irreversible inhibitors has classically relied on the alkylation rate constant, yet this metric fails to resolve its fundamental constituent rate constants, which drive reversible binding kinetics and affinity complex inactivation. In other cases, extremely rapid association processes, which can approach the diffusion limit, also remain challenging to measure. To address these limitations, a practical kinetic rebinding assay is introduced that may be applied for kinetic screening and characterization of compounds. The new capabilities afforded by this probe-based assay emerge from mixed-phase partitioning in a flow-injection configuration and have been implemented using label-free biosensing. A finite element analysis-based biosensor model, simulating inhibition of rebinding within a crowded hydrogel milieu, provided surrogate test data that enabled development and validation of an algebraic model for estimation of kinetic interaction constants. An experimental proof-of-principle demonstrating estimation of the association rate constant, decoupled from the dissociation process, provided further validation.

A rebinding-assay for measuring extreme kinetics using label-free biosensors John G. Quinn In vitro kinetic measurements allow mechanistic characterization of binding interactions and are particularly valuable throughout drug discovery, from confirmation of on-target binding in early discovery to fine-tuning of drug-binding properties in pre-clinical development. Early chemical matter often exhibits transient kinetics, which remain challenging to measure in a routine drug discovery setting. For example, characterization of irreversible inhibitors has classically relied on the alkylation rate constant, yet this metric fails to resolve its fundamental constituent rate constants, which drive reversible binding kinetics and affinity complex inactivation. In other cases, extremely rapid association processes, which can approach the diffusion limit, also remain challenging to measure. To address these limitations, a practical kinetic rebinding assay is introduced that may be applied for kinetic screening and characterization of compounds. The new capabilities afforded by this probebased assay emerge from mixed-phase partitioning in a flow-injection configuration and have been implemented using label-free biosensing. A finite element analysis-based biosensor model, simulating inhibition of rebinding within a crowded hydrogel milieu, provided surrogate test data that enabled development and validation of an algebraic model for estimation of kinetic interaction constants. An experimental proof-of-principle demonstrating estimation of the association rate constant, decoupled from the dissociation process, provided further validation.
In vitro kinetic assays are valuable throughout drug discovery, particularly during hits-to-leads progression, by providing mechanistic discrimination of artifactual binding 1-3 from tractable binding modes. Routine kinetic measurements also allow compounds to be optimized towards a desirable target-specific kinetic profile 4-6 allowing fine tuning of compound properties, including target engagement and residence time for enhanced clinical efficacy 7 . Indeed, routine measurement of direct binding kinetics using real-time label-free biosensors 8,9 provides a practical means of leveraging kinetics for compound prioritization yet the transient kinetics of early chemical matter largely remains beyond the limit of detection. To date, transient kinetics, defined here as affinity complexes that fully dissociate in < 1 s, are measured using low throughput stop-flow based methods 10 , which are impractical for analysis of compound collections in a drug discovery setting. A biosensor-based approach has recently been reported 11 , addressing this limitation but specific system customizations are required to enable routine application. Probe-based kinetic competition assays 12 may be implemented using surface plasmon resonance (SPR)-based biosensors 13 , or other equivalent flow-injection-based biosensors such as grating coupled interfereometry 14 . However, the estimation of transient kinetics and extremely rapid association processes (i.e. association rate constant (k a ) exceeds 5 × 10 7 M −1 s −1 ) for small molecule inhibitor interactions remains challenging using such biosensors 15 and here we introduce a probe-based kinetic rebinding assay to address this unmet need. The rebinding assay is relatively insensitive to both bulk-refractive index mismatches and baseline drift and may be particularly valuable for estimating extreme kinetics of mechanistically complex inhibitors such as irreversible inhibitors. Biophysically realistic virtual biosensors 11,15,16 built using finite element analysis-based computational modeling have provided realistic surrogate data for validation of mechanistic binding models in the past and was adopted for development and validation of the rebinding assay. Advection in bulk flow and diffusion/reaction within the extended hydrogel matrix were modeled as coupled domains of defined volume, where species advection and reaction were computed. Importantly, a finite element-based biosensor model is far too mathematically complex to allow practical estimation of kinetics from actual experimental data. Instead, the virtual instrument simulated rebinding progress curves over extensive parameter ranges to produce surrogate experimental data that enabled the development and validation of an algebraic model. Non-linear least-squares fitting of the algebraic model to both surrogate experimental data, and actual experimental data, returns reliable where R max − R ∝ P, R max is the biosensor response at surface saturation, R is the response at any given time and k t (s −1 ) is the mass transport coefficient 17 . k t defines the rate at which B may enter, or escape, the mass transport boundary layer that forms over the sensing region. The escape time τ = 1/k t is typically in the μs-to-ms regime and confers extraordinary sensitivity to transient kinetics in an inhibition of rebinding format. An expression for k t that accounts for mass transport resistance through the flow cell and through the hydrogel may be defined by the expression where ν c = maximum flow velocity at center of flow channel (m/s), D = diffusion coefficient of analyte in bulk liquid (m 2 /s), h = flow cell height (m) and l = length (m) of functionalized sensing region upstream and including the optically interrogated region. The incorporation of a hydrogel transport resistance term T γ is required to account for hydrogel transport resistance 18 . This term is defined by the height of the hydrogel H gel relative to the mean free path taken by B before being bound, where T γ = Tanh(γ)/γ with γ = H gel / (D gel .K part /(k a ′.P)) 0.5 , K part = hydrogel partition coefficient (unit less) and D gel = diffusion coefficient of B within the hydrogel (m 2 /s). The mass transport coefficient may be expressed in terms of biosensor response as k t ′ = 10 9 .M r_B .k t , where 10 9 is a unit scaling factor (g m/mol) and M r_B is the molecular weight of B (g/mol). Equations (1) and (2) can be used to guide experimental design and imply that rebinding increases with increasing k a ′, increasing P, increasing hydrogel thickness and decreasing flow rate. Equation (2) allows approximation of k t , where hydrogel parameters (e.g. H gel , D gel , K part ) are available, although it is generally estimated as a globally constrained parameter when fitting a mechanistic kinetic model to binding progress curves 15 . The simulated binding curves in Fig. 2 (plot to the left) show that mass transport resistance results in slowing of both association-and dissociation-phase curves as a function of increasing mass transport resistance i.e. increasing Da. However, injection of an excess of A during the dissociation phase ( Fig. 2 (plot on right)) inhibits rebinding of B, restoring the true dissociation rate of AB, while partial inhibition occurs at lower concentrations of A. Inhibition of rebinding model. An optimal balance between kinetic measuring range and sensitivity to inhibition requires relatively low levels of transport limitation (Da < 10) which also maintains monophasicexponential behavior, as shown in the simulated inhibition of rebinding cures in Fig. 3a. Such monotonic behavior also depends on rapid development of a quasi-steady-state with respect to competing pathways acting upon

Figure 2.
Simulated binding response curves for rebinding and inhibition of rebinding. Three mass transport limited binding curves for 10 nM B binding to P tethered within a hydrogel are shown on the left. The response curves were normalized to the maximum saturation response R max . The three simulated curves correspond to k t values of 1 × 10 7 s −1 , 1 × 10 8 s −1 and 1 × 10 9 s −1 producing Da values of 100, 10, and 1, respectively. The simulation was performed using a two-compartment model (see Methods for more details) and the interaction parameters were R max = 20 RU, k a ′ = 5 × 10 7 M −1 s −1 and k d ′ = 0.05 s −1 . Inhibition of rebinding during the dissociation phase at fixed operating conditions is depicted on the right, where A was injected (injection time = t 0 ) at 100 µM and replicated over a serial fivefold range in k a from 1 × 10 5 M −1 s −1 to 7.81 × 10 9 M −1 s −1 to produce dose-dependent inhibition of rebinding. The response observed at the onset (t 0 ) of the inhibitor injection is R 0 . Inhibition curves were simulated using the virtual instrument and finite element analysis engine described in the Method section. .P within the hydrogel that is partitioned by three rate coefficients re-association k a ′.P (s −1 ), inhibition k a .A (s −1 ) and hydrogel escape k t (s −1 ). Therefore, under quasi-steady-state conditions the inhibition curves follow an approximate exponential decay where the change in response is given as  Injection of A produces an inhibition rate that increases k off by lowering rebinding such that k off ≈ k d ′ when fully inhibited. Rather than attempting to estimate transport and reaction fluxes from first principles, the transition state-based model described in Fig. 1 employs a phenomenological rebinding factor α, where k off = k d ′.α and is suitable for estimation of kinetic constants by fitting The rebinding factor α is the degree to which dissociation is slowed due to rebinding and is given by the ratio of rebinding k a ′.P relative to hydrogel escape β. Therefore, α = 1 when rebinding does not exist, or when rebinding is fully inhibited, otherwise α < 1. The partition function f accounts for loss in inhibition of rebinding due to unbinding of AB before escaping the hydrogel. The rate constants associated with f are high relative to k d ′ allowing a quasi-steady-state to be assumed. The response-to-concentration factor G expresses R max in terms of a concentration of P and for many SPR-based biosensors G = 100 RU/g/L. The protein is assumed to be distributed homogenously within the hydrogel and there are experimental methods 19 to estimate this parameter for higher accuracy. In practice, P is selected to possess moderate k d ′, where 1 s ≤ 1/ k d ′ ≤ 300 s, in order to support higher throughput. The assay tolerates wide variation in k a ′ since the reaction flux (k a ′.P) may be modulated by the concentration of P yet higher values (e.g. k a ′ > 1 × 10 5 M −1 s −1 ) are advantageous as this reduces the concentration of B required to achieve a given occupancy.

Estimation of k a for non-transient inhibitors.
When analyzing non-transient binders we assume k d < < k t and therefore f ≈ 1 and can be neglected. The kinetics of BP formation, namely k a ′ and k d ′, are predetermined by conventional direct binding kinetics at low surface density of probe prior to characterizing inhibitors. These kinetic constants are then held constant when analyzing dose-dependent inhibition curves, allowing k a , k d and k t to be readily determined by global fitting Eq. (4). Although not essential, a zero-inhibition curve, where A = 0, may be included to allow k t to be estimated in the absence of inhibition, where α = k t /(k t + k a ′.P). Pre-estimation of k a ′ and k d ′ by conventional direct binding kinetics may be avoided by employing the soluble probe as an inhibitor species in the rebinding assay while also maintaining it as the surface-bound probe. In this case, AB becomes a fully soluble form of BP and hence both k a and k a ′ govern the same interaction occurring in homogenous phase and heterogeneous phase, respectively. The associated kinetic rate constants are related through molecular weight-dependent diffusion scaling, where k a ′ ≈ k a /( M r_P / M r_B ) 1/3 and substitution into Eq. (4) allows k a ′, k d ′ and k t to be estimated from global fitting. These parameters are then held constant when fitting unmodified Eq. (4) to inhibition curves recorded for the inhibitor panel allowing estimation of k a . We generated surrogate experimental data over a wide range in k a , at a fixed concentration of A, in order to determine the relative error and confidence intervals associated with k a -estimation for non-transient inhibitors, where k d < < k t , producing the data show in in Fig. 3. For each curve set the upper and lower limit curves correspond to k off at zero inhibition and k d ' at full inhibition, respectively. These limits define a twofold wider response widow for a tenfold higher k a ′ leading to an increase in measuring range with increasing mass transport limitation (Fig. 3b,c).

Experimental proof-of-principle for estimating k a of non-transient inhibitors.
The probe, in soluble form, was employed as a surrogate inhibitor in order to cross-validate parameter return. In this particular case, the homogenous phase k a may be related to the heterogeneous phase k a ′ by normalizing for differences in collision frequency through diffusion rescaling, where k a ′ = k a /(M r_P /M r_B ) 1/3 , and the experimental data is show in Fig. 4. Direct binding kinetics returned k a = 2.76 ± 0.008 (× 10 6 ) M −1 s −1 for soluble-probe binding to immobilized-target ( Fig. 4a) and after diffusion re-scaling is k a = 1.0 × 10 5 M −1 s −1 , which is within 15% of k a ′ = 8.5 ± 0.004 (× 10 4 ) M −1 s −1 obtained for the reverse format (Fig. 4b), where soluble-target was bound to immobilized-probe. More importantly, k a = 2.76 ± 0.008 (× 10 6 ) M −1 s −1 for direct binding of probe and is within 10% of k a = 3.07 ± 0.01 (× 10 6 ) M −1 s −1 returned from inhibition of rebinding (Fig. 4d), thereby cross-validating the results from these reversed assay formats. The assay is compatible with a throughput of > 500 inhibitors/day using a biacore 8 K + (Cytiva Inc), assuming four concentrations per inhibitor and may be performed in singleton for higher throughput screening (e.g. fragment libraries).
Estimation of transient kinetics. In principle, it is possible to maintain inhibition of rebinding when k d > k t through rapid alkylation of the inhibition complex, where k inact > k d and k inact > k t , leading to a partitioning function z = (1 + k inact /k t ). However, such high k inact values are highly unfavorable in drug discovery due to their non-specific alkylation potential and therefore we consider only k inact < 1 s −1 , allowing z to be neglected. Therefore, estimation of k a and k d for transient irreversible binding remains identical to reversible inhibitors. To illustrate this, we added an irreversible inhibition complex (AB*) to the computational model such that formation of an irreversible inhibition complex (AB*) proceeds at rate constant k inact , where dAB*/dt = k inact .AB. Inhibition curves at six k a values were replicated (n = 8), at four transient binding levels (0.1 ≤ k d ≤ 10 s −1 ) with, and without, inclusion of irreversible alkylation (k inact = 1 s −1 ) for a total of forty eight separate conditions. As shown in Fig. 5a, the resulting forty eight inhibition curves superimpose almost perfectly at each k a value, indicating that the rebinding assay may be expected to return k a estimates using Eq. (4) without inclusion of a partition function www.nature.com/scientificreports/ As shown in Fig. 5b, inhibition of rebinding will decrease for transient inhibitors when k t is in the same order, or less than k d and follows a partition function f = 1 1+k d ./k t , where total hydrogel escape is β = k t + f.k a .A . It also implies that the measurable k d -range may be increased by employing conditions that promote a wider range in k t , (see Eq. (2)) and modulating secondary non-specific transient interactions between the hydrogel and the target. The surrogate experimental data in Fig. 5c,d assume a dissociation phase beginning at 10% of saturation R 0 = 0.1.R max and exhibit low systematic error (< 5%) over 2-orders for simultaneous estimation of k a and k d . Higher R 0 /R max results in higher measurement error, as shown in Fig. 5e,f. Multiple dissociation phase curves were generated over a range in R 0 /R max , were response-normalized (Fig. 5e) and show a maximum divergence of ~ 6% occupancy at R 0 /R max = 1 with negligible divergence for R 0 /R max ≤ 0.25. Propagation of this systematic deviation into error in kinetic parameter return was evaluated by fitting Eq. (4) to surrogate experimental data over a wide range in R 0 /R max , as shown in Fig. 5f. The lowest systematic error was observed for k a -estimation when k d was held constant, resulting in a maximum of 13% underestimation at R 0 /R max = 1, while < 4% error  Plots (c, d, f) include error bars ± parameter fitting error but are low in magnitude and are not visibly resolved from each plotted data point at these font settings. Note that parameter-fitting error is a measure of the amount of information in the curves available to define a unique parameter value and is distinct from parameter recovery error, which compares the estimated value to the true value. It is possible to return low parameter fitting error while parameter recovery error may be relatively high. www.nature.com/scientificreports/ was observed at R 0 /R max < 0.2. This error is related to divergence of the dissociation curves, when R 0 /R max > 0.2, causing a ~ twofold increase in systematic deviation when both k a and k d were fit simultaneously, though this remains within an acceptable error range (< 1.3-fold error) for drug discovery applications. These results are expected as the magnitude of heterogeneous rebinding regimes 20 exhibiting multiphasic dissociation increases when dissociation traverses a wider occupancy range.

Scientific Reports
Limit of detection and parameter return for rebinding assay relative to competitive kinetics. Monte Carlo simulations seeded with pairs of pseudo-random kinetic values were generated to compare the solution phase competitive kinetic binding model of Motulsky-Mahan 21 with the rebinding assay given by Eq. (4). For each assay format, the iso-response contours on the top left-hand corner define the sensitivity limit and indicate a broad measuring range. For competitive kinetics shown in Fig. 6a, the diagonal iso-response contours are affinity isotherms that dominate affinity space while vertical iso-response contours are confined to an affinity region composed of tightly bound inhibitors (bottom right-hand corner), indicating k a -driven inhibition. Conversely, as shown in Fig. 6b, vertical iso-response contours are observed for the rebinding assay indicating fully k d -independent k a -determination over the majority of affinity space. Affinity isotherms are confined to extremely transient affinity space because such transient complexes approach steady-state faster than the inhibitor can escape the hydrogel and become subject to partition function f. The kinetic measuring range of both formats were evaluated using Monte Carlo simulations 22 where each respective kinetic model was back-fit to a large set of simulated curve sets produced from each respective parent model and is shown in Fig. 6c-g. For competitive kinetics, the k a -correlation plot in Fig. 6c shows that kinetic parameters are poorly defined over a broad range in k a when k d is transient, while the remainder of the simulations returned reliable k a estimates, as shown in Fig. 6d. Furthermore, the k d -correlation plot in Fig. 6e also shows poor k d estimates for transient binders. In contrast, the k a -correlation plot for the rebinding assay shown in Fig. 5f indicates that k a remains well defined over the full k a -range and included highly transient binders. In addition, the associated k d -correlation plot shown in Fig. 6g indicates that reliable k d estimates are returned for transient binders, consistent with the results obtained for surrogate data generated from the full computational model and fitted to Eq. (4), as shown in Fig. 5d,f.

Discussion
The measurement of transient inhibitor kinetics in a practical drug discovery setting has not been enabled despite three decades of real-time, label-free technologies. For example, early chemical matter is often transiently bound, with rapid development of a steady-state plateau that is devoid of kinetic information. In general, steady-state dose response curves are prone to artifacts and it is self-evident that mechanistic discrimination of high quality hits based on exceeding a threshold k a would be valuable in prioritizing compounds. The time required to fully displace the volume of the flow cell (> 0.1 s) defines the upper kinetic limit of available biosensors and system modifications that overcome this limit are not yet commercially available 11 . This longstanding unmet need has been addressed here without ultra-fast injection/detection systems by exploiting analyte rebinding within a crowded receptor environment, mixed-phase partitioning in a flow injection configuration and development of an algebraic model (Eq. (4)) of these interdependent processes.
A practical rebinding model for flow injection-based biosensors does not exist and is complicated by the many physical parameters that define the overall system. For example, thick hydrogels containing high concentrations of binding sites may be required for the rebinding assay and a hydrogel resistance term Tγ 17 is needed (see Eq. (2)) to account for the associated increase in mass transport resistance. Equations (1-3) show that rebinding is critically dependent on k t yet it remains impractical to pre-estimate because it requires precise estimates of hydrogel parameters (e.g. H gel , D gel , K part ) that are usually unavailable. In conventional biosensing, rebinding is an interference that compromises kinetic analyses and injection of a soluble form of bound ligand during the dissociation phase has occasionally been employed to inhibit rebinding in an attempt to recover the true dissociation rate constant from simple 1:1 binding curves. Indeed a semi-analytical numerical model was reported 23 for this rate recovery application but was unsuitable for more general applications. While Eq. (4) has been developed for kinetic analysis of inhibitor-target interactions, it may also be employed for this rate recovery application by pre-estimating k a , k d and k t while solving for k d ′.
Equation (4) assumes a phenomenological encounter complex 24 that is partitioned between alternative reaction paths and follows from a recently reported self-rebinding model 25 for bulk solution phase interactions. Equation (4) extends beyond self-rebinding by accounting for rebinding within a receptor-crowded hydrogel and with mixed-phase partitioning in a flow injection configuration. In common with the well know two-compartment model 15 , k t accounts for mass transport resistance but in our rebinding configuration target is partitioned between surface rebinding and escape as a function of k d / k t , thereby providing extraordinary sensitivity to transient kinetics (e.g. ≤ ms-time domain). This obviates the need for ultra-fast opto-electronics and enables k d -independent estimates of k a for non-transient complexes. Partitioning is strongly dependent on hydrogel dimensions and relative spacing of interactants. Briefly, the escape rate of B and AB through the diffusion boundary layer is k t and the probability that B..P reforms BP before exiting a hydrogel of height H gel is a function of the mean free distance d traveled by liberated B before being rebound 17 and is given by P r = 1-exp[-(H gel /d)], where d = (D gel .K part /( k a ′.P)) 0.5 . In the case that d < < H gel , as for a thick hydrogel and/or a high reaction flux coefficient (k a ′.P), then P r ≈ 1 and reformation of BP is favored. However, when d > > H gel then P r ≈ 0 and escape of B from the hydrogel is favored.
In agreement with Eqs. (1) and (2), surrogate-rebinding data showed that k t increased exponentially with decreasing hydrogel thickness H gel and increased at higher flow velocities with expected ν c 1/3 scaling. The rebinding model returns k a estimates that are fully independent of steady-state, or k d , while k d < < k t , and this condition was adopted for analysis of the surrogate-rebinding assay data in Fig. 3. The two curve sets in Fig. 3a   www.nature.com/scientificreports/ identical simulation parameters other than k a ′, giving a tenfold difference in mass transport limitation. The analysis shows that the rebinding assay measured inhibitor kinetics over a wide range (2.5-3.5 orders) at a single test concentration (Fig. 3a,b). The analysis supports operation in a moderate mass transport limited regime since this allowed inhibition to occur at lower inhibitor concentrations (Fig. 3b) thereby avoiding compound solubility artifacts. This elevated parameter return error (Fig. 3b,c) but k a estimates were nevertheless returned with < 1.2-fold relative error and with narrow (< 5%) confidence limits (assuming 95% confidence interval). While the primary objective of the current work was to develop, and validate, Eq. (4) using surrogate data from the computational model, we also demonstrate an experimental proof-of-principle as shown in Fig. 4. Indeed, k a estimated from Eq. (4) was in good agreement with estimates from direct binding and was cross-validated within the rebinding assay itself by also introducing the probe as the inhibitor species. Experimentally, this proof-of-principle shows that the rebinding assay format is comparable to direct SPR binding in terms of experimental complexity and data analysis. As already mentioned, when k d < < k t , then AB is relatively stable such that it remains independent of both k d and any slower coupled kinetic reaction rates. The data in Fig. 5a implies that k a may be measured for complex interaction mechanisms, such as irreversible inhibition, providing an opportunity to complete full mechanistic analysis with estimation of the fundamental kinetic rate constants (i. e. k d , k a , k inact, ). However, when k d = k t then affinity partitioning results in 50% loss in inhibition (Fig. 5b) and transient partitioning ultimately defines the limit of detection of the assay. However, such partitioning is sensitive to transient kinetics allowing k a and transient k d to be measured over 2-orders (Fig. 5c,d). Furthermore, while the assay performs best when R 0 < 0.1.R max , it was found that holding R 0 = R max (Fig. 5e,f) incurred < 1.3-fold relative error in kinetic estimates, an acceptable tolerance for many drug discovery applications. Solution-phase competitive binding kinetics may be described by the analytic model of Motulsky-Mahan 20 and assumes formation of a "hot" inhibition complex containing a tracer compound to indirectly report the evolution of the inhibition complex. This format can be replicated for surface sensitive biosensors 13 by injecting a mixture containing soluble probe and inhibitor over a target-coated surface to generate a resolvable binding response, assuming a significant refractive index difference exists for binding of the competing probe relative to inhibitor binding. In solution-phase competitive kinetics, the reactions evolve towards steady-state, whereas the rebinding format maintains a k a -driven regime over a wide range. To explore these contrasting properties, sensitivity analysis was performed for both formats using Monte Carlo-like simulations seeded with random kinetic parameter values, and the results are shown in Fig. 6a,b. The resulting affinity plots show that the rebinding assay reports k a independently of affinity since the affinity isotherms that are typical of competitive kinetics (Fig. 6a) are absent for rebinding (Fig. 6b), except in the case of transient binders, which are detected at 33-fold higher sensitivity. Furthermore, transition from an affinity-dependent regime to an affinity-independent regime can only occur for solution-phase competitive kinetics when the inhibition complex is extremely stable (k d < 1 × 10 −4 s −1 ). For the rebinding assay, the affinity-independent regime dominates and can be further extended by increasing k t and can be accomplished by limiting self-self-interactions, non-specific hydrogel interactions, and enhancing convective/diffusive mass transport within the hydrogel-flow cell system. Monte Carlo simulations were also performed to test the limits of kinetic parameter return. The data showed that solution phase affinity (Fig. 6c-e) did not return reliable kinetics for transient binding (k d > 10 s −1 ), whereas the rebinding format accurately returned both k a and k d for such transient binders (Fig. 6f, g), while reporting affinity independent k a for non-transient binding. In summary, a rebinding assay exploiting a flux balance of target partitioned between rebinding and escape from within a hydrogel film was developed. The method is well suited to resolving both k a and k d of transient affinity complexes with adequate throughout for practical applications, which remains challenging in a contemporary drug discovery setting. An experimental proof-of-principle demonstrated estimation of k a that was independent of both k d and steady-state thereby establishing the feasibility of measuring extremely rapid association kinetics for non-transient binding complexes. A comprehensive application of this method in drugdiscovery is currently in progress and will be the subject of a publication in the near future.

Methods
Experimental proof-of-principle for estimating k a of non-transient inhibitors. Assays were conducted using a Biacore S200 (GE Healthcare Bio-Sciences AB, SE-751 84, Uppsala, Sweden) with analysis temperature set to 20 °C. All binding curves were acquired and plotted at 40 Hz and the baseline noise was approximately 0.03 RU (root mean square) thereby producing a high signal-to-noise ratio for all binding curves. All reagent coupling kits and sensors were from GE Healthcare. A biotinylated-avi-tagged 22 kDa target protein was expressed recombinantly and purified in-house using standard protocols. The probe molecule (M r_P = 737 Da) was a PEGylated compound with moderate affinity for the target where a terminal primary amine on the PEG linker allowed coupling to a CM5 sensor chip through standard EDC/NHS covalent linkage chemistry. All experiments were performed using an assay buffer containing 50 mM 4-(2-hydroxyethyl)-1-piperazine-ethanesulfonic acid (HEPES), pH 7.5, containing 0.15 M sodium chloride, 0.2 mM tris(2-carboxyethyl)phosphine (TCEP), 0.1% polyethyleneglycol (M n ~ 4 kDa) and 1 mg/ml carboxymethylated dextran (M n ~ 10 kDa).
Kinetics of probe binding to tethered target. Target was captured onto a sensing surface of a series S SA sensor chip. Probe was injected for 30 s in duplicate over a serial doubling dilution range from 1 µM to 31 nM prepared in assay buffer. The binding curves were double referenced and fit to a two-compartment 1:1 interaction model for estimation of k a and k d . All parameters were fit globally.
Probe coupling. Channels 2 and 4 of a series S CM5 sensor chip were activated in-situ for 8 min using standard EDC/NHS. The chip was undocked, rinsed with buffer and 40 μl of probe solution containing 1 mM probe, diluted in a 1:1 (v/v) solution of DMSO:1 M HEPES( pH 7.5), was pipetted onto the sensing surface and incubated for 2 h at room temperature. The surface was rinsed in buffer and blocked with 1 M ethanolamine for 10 min. The chip surface was rinsed with 100% DMSO, 20 mM NaOH and ultrapure water, dried, and re-docked. www.nature.com/scientificreports/ Kinetics of target binding to tethered probe. Target protein was diluted in assay buffer and was injected for 60 s in duplicate over a serial-doubling dilution range from 1 µM to 31 nM. The binding curves were double referenced and fit to a two-compartment 1:1 interaction model for estimation of kinetics. All parameters were fit globally.
Rebinding assay curves. A fresh probe-coated CM5 sensor chip was prepared, using the above probe coupling method but with a 20 min probe-solution contact time, which lowered the R max relative to the previous 2 h exposure. 500 nM target was injected over channel 2 (probe-coated) at 30 μL/min for 60 s followed by dissociation. This was repeated (8 replicates) and the target was allowed dissociate from the surface between replicate cycles. The inhibitor (i.e. soluble probe) was injected over channels 1 and 2, in duplicate, at 0 μM, 0.016 μM, 0.6 μM and 4 μM with assay buffer as diluent. This inhibitor injection commenced for each of these target-loading injections when the dissociation response had decreased to 118 RU. Fig. 2. Mass transport limited binding curves for both association and dissociation phases were stimulated using Biaevaluation 4.1.1 (Cytiva life sciences) by selecting a twocompartment mechanistic interaction model and numerical integration of the associated coupled set of ordinary differential equations (ODE) given as Simulation details associated with Fig. 6. The simulation parameters for each assay format were matched in terms of R 0 , R max and kinetic range. Kinetic parameter values were chosen pseudo randomly over many orders in both k a and k d to producing large sets of simulated response curves. The range in simulated k a values were defined by the limit 4 ≤ Log(k a ) ≤ 9 as this range represents diffusion limited association reactions with an upper limit bounded by the maximal diffusion rate in aqueous phase and a lower limit bounded by the onset of conformationally gated binding. The k a was varied over 10-orders spanning the upper and lower limits of practical importance in drug discovery, with limit -6 ≤ Log(k d ) ≤ 4. Competitive kinetic inhibition curves were simulated over 6-serial tenfold dilutions of inhibitor from 1 mM in order to compensate for intrinsically lower measuring range. Rebinding curves were simulated over three serial tenfold dilutions of inhibitor from 1 mM, where each was repeated at two injection flow rates. Mote Carlo simulations were generated using Graphpad Prism version 9.0.0 (Graphpad software LLC, 7825 Fay Avenue, Suite 230, La Jolla, CA, 92037, USA), with baseline noise of 0.06 RU added to response curves to mimic experimental baseline noise. Each simulated curve set was then back-fit to its parent model in order to test the accuracy of parameter return, measuring range and detection sensitivity. Both models were applied using matched constraints when fitting. k a and k d were fit globally and k t was fixed. The number of simulations included in each plot varies because the simulation range exceeded the desired ranges a limitation of the adapted Monte Carlo routine. All parameters were fixed in both simulations other than k a and k d . In a drug discovery setting, the minimal resolvable response change defines the limit of detection for weak binding compounds and extending this range is of considerable value. Only the highest concentration of 1 mM was required for sensitivity analysis, where the response was measured over an average of 4 points at the end of each 5 min injection and a cut-off response of 0.5 RU (signal-to-noise ratio = 8) was selected. The resulting responses were plotted as response isotherms on an affinity space plot where x-axis = Log(true-k a ) and y-axis = Log(true-k d ).

Simulation details associated with
Virtual instrument. The specifications given apply to all simulation data unless otherwise stated. The three-dimensional hydrogel is well approximated using a two-dimensional geometry because contributions near the flow cell walls can be neglected. This is the case because the flow cell is thin relative to its width. The data collection rate was 1 Hz and baseline noise equivalent to 0.03 RU (root mean square) was added to response curves to mimic actual instrument performance.
Hydrogel modeling. Hydrogels are highly complex environments making it difficult to accurately determine species diffusion rates and related mass transport effects. For example, diffusion might be > tenfold slower at high protein concentrations and some size exclusion may also occur further slowing diffusion. For example, a mass equivalent to > 10,000 RU of protein can readily be bound within a hydrogel. This would correspond to a volume fraction of > 10% (v/v), assuming a specific volume ϕ V (ml/g) = 0.754 and would be expected to increase viscosity > tenfold and decrease diffusion by at least this factor 26 . A 200 nm thick hydrogel was modeled as an area containing a homogenous density of hydrogel grafted to the sensing surface that decreases rapidly at the hydrogel-liquid interface according to a hydrogel density function density = 1-Exp(-100*(1-z)), where z is the www.nature.com/scientificreports/ dimensionless relative height of the hydrogel. The concentration of P is assumed to be scaled by the hydrogel density and the diffusion coefficient of all species tethered to the hydrogel (i.e. P, BP) is assumed to be zero. Diffusion of all unbound species inside hydrogel is assumed to be twofold lower due to a twofold increase in viscosity within the hydrogel relative to the bulk liquid and soluble species are subject to molecular weightdependent partitioning. Therefore, parameters related to mass transport of soluble species inside the hydrogel are defined as follows. Diffusion coefficient of A = D A = 5 × 10 −10 (m 2 /s), diffusion coefficient of B = D B = D A (M r _ B / M r _ A ) 1/3 , diffusion coefficient of A inside hydrogel = D gel = 2.D A .K part , where the hydrogel partition coefficient for A = K part = Exp(-10 −3 * M r_A 2/3 ), diffusion coefficient of all soluble species containing B (i.e. B, AB, AB*) inside hydrogel = D gel = 2.D B .K part , where the hydrogel partition coefficient for B = K part = Exp(-10 −3 * M r_B 2/3 ). The initial conditions for the simulation include addition of tethered P some fraction of which is in the form of affinity complex BP before the onset of the inhibitor injection. The inhibitor injection was simulated as a sample pulse entering from one end of the rectangular flow cell and exiting at the opposite end, where the tethered hydrogel film is located at one of the flow cell walls and is parallel to the direction of flow.
Finite element analysis. Coupled ODEs were solved numerically coupled to the master equations, which are partial differential equations (PDEs) governing flow, advection, diffusion, and reaction for a 20 µM thick flow cell housing a sensing region containing a hydrogel film functionalized with P. The entire geometry was discretized in space and solved over incremental time periods to generate surrogate experimental data. Comsol multiphysics 5.1 (COMSOL AB, Tegnérgatan 23, SE-111 40, Stockholm, Sweden) was used to perform all numerical simulations. A computational model replicating the flow injection-based biosensor system depicted in Fig. 1 was created. Typical microfluidic channels employed in biosensors have high aspect ratios where side walls are far apart relative to the top/bottom walls allowing microchannel width to be neglected reducing the model to a cross-section through the microchannel. The two dimensional flow cell geometry housing a hydrogel film grafted to the flow cell sensing region was meshed with > 14 k elements. This mesh was optimized until no detectable change was observed in the simulation output and included a higher density of elements at the hydrogel interfacial boundaries arriving at 6762 elements over the hydrogel domain. The incompressible form of the Navier-Stokes equation was used to solve the two-dimensional velocity profile through the channel, assuming steady-state, at constant flow rate and at atmospheric pressure. The initial velocity at the walls u wall = 0, the inlet velocity was variable and was defined by solving for the velocity vector field over the full domain where ρ is the density, p is the pressure and μ is the dynamic viscosity.
The flow velocity vector field was coupled to the steady-state advection/diffusion equation for each dilute species to solve for the analyte concentration field in the bulk flow domain and hydrogel domain.
Here D is the diffusion coefficient, c is the concentration of a given species and R is a reaction term associated with that species. Initially the analyte concentration in the microchannel c = 0. At the inlet the initial analyte concentration profile along the microchannel height was defined by multiplying the concentration by a rectangular function to simulate a continuous injection of sample for a given contact time. P was assumed to be distributed within the hydrogel domain tethered to the sensing surface. The associated reactions within the hydrogel domain between soluble species with P distributed within the hydrogel and the formation of non-tethered affinity complexes AB were defined as ODEs coupled to the advection/diffusion Eq. (6) and are given as where k a and k d , are the forward and reverse kinetic interaction constants for the AB complex and k a ′ and k d ′, are the forward and reverse kinetic interaction constants for the BP affinity complex. This set of ODEs describe reversible affinity complex formation and formation of an irreversible complex (AB*) required the following ODEs The time-dependent change in species accumulation was found from a species flux balance over the hydrogel domain and the simulation was performed in time-stepping mode in order to produce reaction progress curves. The accumulation of affinity complex was expressed in terms of an equivalent biosensor response where 100 RU = 1 mg/ml of a given species.
Curve fitting statistics. Microsoft Excel and Biaevaluation (GE Healthcare Bio-Sciences AB) were employed for data processing. Graphpad Prism version 9.0.0 (Graphpad software LLC, 7825 Fay Avenue, Suite 230, La Jolla, CA, 92,037, USA) was employed for all plots other than Fig. 6 a,b, which were ploted using DPlot www.nature.com/scientificreports/ Version 2.3.5.7 (HydeSoft Computing, LLC). Graphpad Prism enabled fitting of binding interaction data to interaction models by nonlinear regression, and the associated statistical methods to confirm goodness of fit and confidence in parameter estimates are well established 27 . Statistical parameters such as the standard error of the fit (SE) associated with a given parameter returned in the fit were used to report confidence in parameter estimates. The SE is a measure of the information content of the data and specifies the degree to which the curves define the parameter value from the fit. Values < 5% indicate high confidence and values > 10% indicate that the parameter is poorly defined. The goodness of fit between a model curve and an experimental curve is described by χ 2 when the number of data points is high and by a regression coefficient R 2 when the number of values is low. %χ 2 is the square of the averaged residual response difference expressed as a percentage of maximum response recorded for the curve set. Typically high quality fits will produce χ 2 values < 5%. Occasionally χ 2 may be within acceptable limits but the fit may remain questionable if residuals are not distributed randomly. Curves generated by numerical simulation follow deterministic algorithms and therefore do not require replicates. Global parameter fitting refers to constraining a parameter value to a single global value over the entire curve set.