Unlocking latent kinetic information from label-free binding

Transient affinity binding interactions are central to life, composing the fundamental elements of biological networks including cell signaling, cell metabolism and gene regulation. Assigning a defined reaction mechanism to affinity binding interactions is critical to our understanding of the associated structure-function relationship, a cornerstone of biophysical characterization. Transient kinetics are currently measured using low throughput methods such as nuclear magnetic resonance, or stop-flow spectrometry-based techniques, which are not practical in many settings. In contrast, label-free biosensors measure reaction kinetics through direct binding, and with higher throughout, impacting life sciences with thousands of publications each year. Here we have developed a methodology enabling label-free biosensors to measure transient kinetic interactions towards providing a higher throughput approach suitable for mechanistic understanding of these processes. The methodology relies on hydrodynamic dispersion modeling of a smooth analyte gradient under conditions that maintain the quasi-steady-state boundary layer assumption. A transient peptide-protein interaction of relevance to drug discovery was analyzed thermodynamically using transition state theory and numerical simulations validated the approach over a wide range of operating conditions. The data establishes the technical feasibility of this approach to transient kinetic analyses supporting further development towards higher throughput applications in life science.

Modeling dispersion gradients. Analyte dispersion gradients form sigmoidal rise/fall profiles and may be modeled using an ODE, or an analytical expression 1 such as; . Equation (S1) is a biophysical expression relating the relevant physicochemical parameters and accommodates model fitting over a range of experimental conditions. Here we assume constant experimental conditions and avoid the computationally intensive Erf(z) of equation (S1) by using the following logistic equation to model dispersion, where a, b, c, d are gradient fit coefficients. The fit coefficients may be pre-calculated by fitting equation (S2) to bulk refractive index (RI) dispersions but may also be estimated during fitting of the threecompartment model to binding curves. The time for the analyte dispersion to rise by 95% is tgrad = a19 1/b and for the fall is tgrad = c19 d . From a practical perspective tgrad is readily controlled by flow rate selection and by employing a range of injection loops geometries. In the limit of tgrad <<  the uniform analyte concentration assumption holds such that the three-compartment model effectively reduces to a twocompartment model. When tgrad >>  then occupancy develops slowly relative to  allowing systematic errors, such as baseline drift, to impact the data to a greater degree but this remains negligible for transient dispersion injections because tgrad < 1s.
Influence of mass transport limitation. When the surface reaction flux coefficient kr = ka (Rmax-R) is in the same order as kt, or higher, then mass transport limitation MTL = kr/(kr+kt) will be significant and a compartment model that includes kt is then required for reliable kinetic analysis.

(S3)
where uc is the maximum velocity through the center of the microchannel, D is the analyte diffusion coefficient and can be estimated from the Einstein-Stokes-Sutherland equation 2 , h is the height of the flow cell and l is the length of the sensing region. Additional target-coated areas upstream are undesirable as this increases %MTL for a given Rmax but when present this additional analyte binding sink requires an additional area correction factor.
Data analysis and statistics. Binding curves can be fit with a mechanistic interaction model by numerical integration of a coupled set of ordinary differential equations (ODE) combined with nonlinear squares curve fitting. Microsoft Excel and Biaevaluation (GE Healthcare Bio-Sciences AB) were employed for data processing. The computer generated image of the full injection-binding process shown in Fig. 1a was generated in Comsol. Model fitting to binding curve sets was performed using Biaevaluation (GE Healthcare Bio-Sciences AB) and all images of fitted binding curve sets were taken directly from Biaevaluation. Graphpad Prism version 6 (GraphPad Software, Inc. 7825 Fay Avenue, Suite 230, La Jolla, CA, 92037, USA) was employed to plot Fig. 1b and Fig. 1d, parameter estimation error in Fig. 2 and Fig. 4 and for fitting thermodynamic models in Fig. 3b and presenting the transition state reaction plot in Fig. 3c. Microsoft excel was used to plot the evolution of the analyte concentration profile across the height of the flow cell at two flow rates (Fig. 4a). Curve fitting programs enable 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 3 . 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 and approaches the baseline noise for the best fits. 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. A large systematic residual deviation can indicate inappropriate model fitting as observed in Fig. 2b.
Curves generated by numerical simulation follow deterministic algorithms and reproduce without error and therefore do not require replicates. However, simulation over a sufficiently wide parameter range is important in evaluating simpler models specified as ODEs and we choose ranges that were representative of practical experimental ranges. For a well-optimized assay, label-free binding technology can reproduce binding curves to a very high accuracy as observed in Fig. 3a. Each curve set contains duplicate curves at 6 concentrations. Duplicates are not averaged into an average curve but instead the model is fit simultaneously to all curves in the set producing 12 simulated curves that when well fit superimpose onto the experimental curves. All interaction parameters are constrained to a single global value over the entire set providing a robust test of the model. The resulting parameter estimates therefore represent the sum of multiple replicated measurements varied in concentration at a given temperature. Biaevaluation 3.0 was used to fit interaction models by writing a custom model using the general model fit editor in the form of the coupled ODE for a two-compartment model as given in Fig. 1c. and using equation (S2) to model dispersion and equation (S3) to pre-calculate kt.
Correlation analysis. The correlation of model parameters was determined using Graphpad Prism by simulating a set of model curves using the three-compartment model where ka = 1 x 10 6 M -1 s -1 , kd = 10 s -1 , Rmax = 60 RU and kt = 1.15 x 10 8 m/s. The covariance matrix was automatically calculated by the program and was plotted as correlation coefficients in figure 2i.
Numerical simulations. This complex mathematical discipline requires specialists working with dimensionless numbers in order to collapse the near infinite possibilities of parameter space into the simplest governing theory applicable to a given regime of operation 4 . However, we confine our simulation studies to model specific experimental conditions rendering the added level of abstraction imposed by dimensionless analysis unnecessary. We employ terminology that is familiar to the drug discovery community in order to aid understanding. 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 biosensor system given in Fig. 1a 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 thereby reducing the model to a cross-section through the microchannel. The microchannel has height h = 50 M and length l = 21 mm with a short flow cell at the exit region. A mesh composed of > 200 k elements was applied with exponentially closer spacing at the wall boundaries and the number of elements was increased until no detectable change was observed in the simulation output. The serpentine microchannel was included in Fig. 1 and Fig. 2 in order to generate a dispersion profile thereafter dispersion was defined using an analytical dispersion for all other data sets. Equation (S2) was used because the transient kinetics ( = 10 ms) being simulated demanded a very high mesh density to ensure that the simulation was stable over all parameter ranges evaluated and resulted in impractically long computational times when applied over the entire serpentine microfluidic channel. The incompressible form of the Navier-Stokes equation was used to solve the two-dimensional velocity profile through the channel, assuming steady-state, a constant flow rate and atmospheric pressure. The velocity at the walls uwall = 0, the inlet velocity uinlet = 0.1 m s -1 , then 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 a dilute species to solve for the analyte concentration field in the bulk flow.

∇. (−D∇c) + . ∆ = (S6)
Here D is the diffusion coefficient, c is analyte concentration and R is a reaction term. 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 suitable scaling function (e.g. rectangular, sigmoidal, step, wave function). Target was assumed to be bound at the wall of the flow cell adjacent to the microchannel exit and required a surface reaction flux at that wall. The associated surface reactions were modeled by a surface molar flux, Nt,i according to = − ∇ (S7) where Ds,i is the surface diffusion coefficient for species i at concentration cs,i. The governing equation for surface reaction is where Rs is the surface flux balance (mol m -2 ) with Ds = 0. The reversible binding reaction at the surface was given by a mechanistic 1:1 pseudo-first order interaction model identical in form to equation (1) given in terms of biosensor response but is presented here in terms of species concentrations where cmax is the total concentration of target bound to the sensing surface, cs is the concentration of affinity complex at any given time.
The binding flux of analyte to the target-coated sensing surface was balanced by a coupled flux loss from the bulk liquid. The time-dependent change in analyte accumulation was found from a surface flux balance at the sensing surface where the simulation was performed in time-stepping mode. The accumulation of affinity complex was expressed in terms of an equivalent biosensor response R(t) = cs(t).Mr.10 6 . Parameter values employed in numerical simulations are given in table S1.  Table S1. Parameter values employed in numerical simulations.