Multi-temperature experiments to ease analysis of heterogeneous binder solutions by surface plasmon resonance biosensing

Surface Plasmon Resonance (SPR) biosensing is a well-established tool for the investigation of binding kinetics between a soluble species and an immobilized (bio)molecule. While robust and accurate data analysis techniques are readily available for single species, methods to exploit data collected with a solution containing multiple interactants are scarce. In a previous study, our group proposed two data analysis algorithms for (1) the precise and reliable identification of the kinetic parameters of N interactants present at different ratios in N mixtures and (2) the estimation of the composition of a given mixture, assuming that the kinetic parameters and the total concentration of all interactants are known. Here, we extend the first algorithm by reducing the number of necessary mixtures. This is achieved by conducting experiments at different temperatures. Through the Van’t Hoff and Eyring equations, identifying the kinetic and thermodynamic parameters of N binders becomes possible with M mixtures with M comprised between 2 and N and at least N/M temperatures. The second algorithm is improved by adding the total analyte concentration as a supplementary variable to be identified in an optimization routine. We validated our analysis framework experimentally with a system consisting of mixtures of low molecular weight drugs, each competing to bind to an immobilized protein. We believe that the analysis of mixtures and composition estimation could pave the way for SPR biosensing to become a bioprocess monitoring tool, on top of expanding its already substantial role in drug discovery and development.


Boundaries on the Identified Parameters
We impose the following boundaries on the parameters to be identified: 1. Equilibrium contributions must be between 0 and 1: 0 ≤ , , ≤ 1 ∀ , 2. The dissociation rates are positive: , , 1 > 0 ∀ 3. The dissociation rates increase with temperature. From the Eyring equation ), it follows that: 4, < 0 ∀

Constraints on the Identified Parameters
We consider three types of constraints on the identified parameters: 1. The sum of the equilibrium contributions is always 1 for all mixtures and temperatures ( • linear equality constraints): In terms of the adjusted parameters , 1 , , we have: , =1 2. The association rates , increase with temperature. The expression of this constraint with respect to the identified parameters is described lower. This leads to • ( − 1) nonlinear inequality constraints.

3.
The affinity , varies monotonously with respect to temperature (either increases or decreases, depending on if the dissociation or association rate temperature dependence is greater). The expression of this constraint according to the identified parameters is described lower. This leads to • ( − 2) nonlinear inequality constraints.

Constraints on the Association Rates
The association rates , increase with temperature. One can express , , via the Eyring equation-derived definition of the dissociation rates with a reference temperature 1 : , , = , , , , 1 1 exp ( 4, ( 1 − 1 1 )) We take two temperatures and . If > , it follows that , , > , , ∀ : We know that the affinity can be obtained from the equilibrium contribution: , , = , 1 , ( , By inputting this definition in the inequality, we obtain: As , 1 and , are not temperature dependent, one only has to verify that the numerator , = , 1 , ( , , ) 1 , only increases or decreases with temperature. We take three temperatures > > to obtain the following constraint: This gives a total of • ( − 2) nonlinear inequality constraints.

Starting Point of the Optimization Routine
The order of magnitude of the , , is known as they are required to be comprised between 0 and 1 and sum to one. Hence, an obvious starting would be 1/ for all analytes. An appropriate starting point for , , 1 and 4, is not as obvious, a priori. SPR biosensors can measure dissociation rates ranging from approximately 10 -4 to 0.5 s -1 , and the 4, are a function of the enthalpy of the dissociation pseudo-reaction. Given that the dissociation phase is long enough to reach a null response for all sensorgrams, we suggest a starting point based on the data and a coarse assumption. This will give in turn a coarse starting point, but of a reasonable order of magnitude.
We take the integral of the dissociation phase (ranging from the beginning of the dissociation phase to the end of the sensorgram ) of the normalized SPR sensorgrams ( , ): If is such that the SPR response is approximately zero towards the end of the sensorgram, we obtain: We express the previous equation with respect to a reference temperature 1 to obtain: We consider the case where all analytes share the same , , 1 = and 4, = . We obtain: ) with respect to by considering all sensorgrams (all mixtures and temperatures) leads to estimates of and which are themselves coarse estimates of the , , 1 and the 4, , respectively. This procedure leads to starting values of the right order of magnitude for the suggested optimization routine.

Confidence Intervals on the Kinetic Parameters
The confidence intervals can be computed based on the standard error 7,9,15. For the ℎ parameter, the standard error is given by: Where [ −1 ] , is the ℎ element on the diagonal of the inverse of the Hessian matrix. 2 is a function of the number of identified parameters = • (2 + 1) and of the number of data points : The computation of the Hessian matrix is detailed in the Supplementary Materials of 15 for the case of a single injection temperature. The extension to multi-temperature experiments is given lower. Once the standard error of a given parameter has been obtained, the 100(1 − )% confidence interval for that parameter is computed via Student's T distribution:

Hessian Matrix with analytes and Injection Temperatures
The Hessian matrix can be approximated by the following sum on all time steps of all available sensorgrams (sensorgrams corresponding to different mixtures and temperatures): The derivatives above can be evaluated by solving the following ODEs along with the system of ODE in : Here we present a way to compute the necessary gradient

Where
, and , are the sections of vectors and that correspond to temperature . We consider only those sections because kinetic rates at different temperatures have no effect on the predicted response at temperature and hence corresponding derivatives will be zero. Maximal responses, however, are not temperature dependent. We consider a third matrix: For each time point of each sensorgram, we obtain the gradient , . For injection temperature , it contains the subsections computed using , , , and and zeros for every element corresponding to kinetic rates at other temperatures. Computing the gradient at every time step of every sensorgram in the data set (multiple mixtures and temperatures) is necessary to compute the Hessian matrix.
If bulk effect parameters ( , one per sensorgram) were added to the fitted model, additional elements must be added to , as the parameter vector then contains as many additional elements as there are sensorgrams in the data set. For a given sensorgram , the associated bulk effect parameter is noted , , and , , is equal to 1 for time steps during the association phase and 0 for time steps during the dissociation phase.

Confidence Intervals on the Fractions and Concentrations
A method to evaluate the confidence interval on the estimated fractions is given in 15 for the case of known total concentration. It leads to asymmetrical intervals so that the constraints on the can be taken into account (contained between 0 and 1). We may use the same general method here, while considering , as an additional identified parameter. The method is based on the Fisher F statistic 15,51. Taking into consideration the objective function in (26), we have: Where (̂) is the value of the objective function corresponding to the optimum point (estimated fractions and concentration) and ( )| = 0 is the value of the objective function obtained by optimizing with an added constraint = 0 . The boundaries of the 100(1 − )% confidence interval are such that: is a quantile of the Fisher law with degrees of freedom in the parentheses. To find the upper boundary of the confidence interval of parameter , we restart the optimization by disturbing parameter such that ′ =̂+ ∆ . This is done multiple times with ∆ being progressively larger, until ( )| = ′ satisfies . When this is the case, the current ′ corresponds to the upper boundary of the confidence interval of . The lower bound can be found by using negative values for ∆ . A bisection method may be used to find the proper ∆ , as is detailed in the Supplementary Materials of 15. 43.95 ± 0.08 Table S1: Kinetic parameters identified from single-and multiple-analyte injections of mixtures of known composition. Multi-analyte fits were obtained by fitting the multi-analyte model independently using any combination of three mixtures (A-B-C-D). 95% confidence intervals are given underneath the identified parameters. The estimated affinity is also reported. Association rates ( ) are reported in 10 4 s -1 M -1 , dissociation rates ( ) are reported in 10 -2 s -1 , affinities ( ) are reported in 10 6 M -1 and maximal responses ( ) are reported in RU. A star (*) next to a parameter value indicates that a constraint on the corresponding parameter was active in the fit.  Table S2 : Identified thermodynamic parameters obtained from single-analyte experiments and by fitting the multi-analyte model to any combination of three mixtures (A-B-C-D) at four temperatures. The activation enthalpy and activation entropy were obtained by linearizing the Eyring equation using the identified association and dissociation rates at 12, 16, 20 and 24 o C for each analyte. They are given for both the association (∆ * and ∆ * ) and dissociation (∆ * and ∆ * ) pseudo-reactions. The standard reaction enthalpy and entropy were obtained by linearizing the Van't Hoff equation using the identified affinities of each analyte.  Table S3: Kinetic parameters identified from multiple-analyte injections of mixtures of known composition. Multi-analyte fits were obtained by fitting the multi-analyte model using mixtures B, C and D. Starting points and final estimates for part 2 of the suggested algorithm are reported whether part 1 was performed prior (our algorithm) or not (suboptimal starting point, expected order of magnitude of the parameters). The estimated affinity is also reported. Association rates ( ) are reported in 10 4 s -1 M -1 , dissociation rates ( ) are reported in 10 -2 s -1 , affinities ( ) are reported in 10 6 M -1 and maximal responses ( ) are reported in RU. 2 values are reported below the data set name, in parentheses.  Table S4: Kinetic parameters identified from multiple-analyte injections of mixtures of known composition. Multi-analyte fits were obtained by fitting the multi-analyte model independently using any combination of two mixtures (A-B-C-D). 95% confidence intervals are given underneath the identified parameters. The estimated affinity is also reported. Association rates ( ) are reported in 10 4 s -1 M -1 , dissociation rates ( ) are reported in 10 -2 s -1 , affinities ( ) are reported in 10 6 M -1 and maximal responses ( ) are reported in RU. A star (*) next to a parameter value indicates that a constraint on the corresponding parameter was active in the fit.  Table S5: Identified thermodynamic parameters obtained from single-analyte experiments and by fitting the multi-analyte model to any combination of two mixtures (A-B-C-D) at four temperatures. The activation enthalpy and activation entropy were obtained by linearizing the Eyring equation using the identified association and dissociation rates at 12, 16, 20 and 24 o C for each analyte. They are given for both the association (∆ * and ∆ * ) and dissociation (∆ * and ∆ * ) pseudo-reactions. The standard reaction enthalpy and entropy were obtained by linearizing the Van't Hoff equation using the identified affinities of each analyte.

Thermodynamic Parameters Identified with Two-Mixture Data Sets
-Deviations from Single-Analyte Experiments Figure S1: Deviation between the thermodynamic parameters (∆ * , ∆ * , ∆ * , ∆ * , ∆ 0 and∆ 0 ) derived from single-analyte experiments and parameters derived by fitting the multi-analyte model at four temperatures with any combination of two mixtures (A-B-C-D). A) Average deviation computed over all parameters for all analytes in all fits. B) Average deviation computed over all fits (except C-D) for all parameters of all analytes. Figure S2: Performance of the composition estimation algorithm with known concentration. A) Calculated fractions with respect to actual fractions of the four analytes in each of the 8 mixtures detailed in . These fractions were estimated from sensorgrams at 12 o C only with kinetic parameters identified from data set B-C-D. B) Mean absolute error of estimated fractions for each analyte at each temperature. The composition estimation algorithm was used independently for each set of kinetic parameters corresponding to each data set (all combinations of 2 and 3 mixtures) and each temperature, and mean absolute errors were averaged across all data sets. C) Mean absolute error (averaged across all data sets) with respect to the performance indicator. Annotations on the graph indicate the corresponding temperature. D) Mean absolute error of the fractions with respect to the data set that was used to identify the kinetics depending on if the concentration was considered known or added as a parameter to be estimated. Those were taken at 12 o C only. In this figure, 'singles' refers to single-analyte fits. Figure S3: Profile likelihood of part 1 of the parameter identification algorithm with data set B-C. Each parameter was disturbed 21 times to fixed values ranging from 0.75 times to 1.25 times its optimal value (X-axis of the plots). For each disturbance, all the other parameters were reoptimized. The effect of disturbing each parameter on the value of the objective function is reported on the Y-axis of the plots. The red lines indicate the boundaries of the 95% confidence intervals of each parameter, as obtained with the method presented in section 0. A star (*) marks the location of the optimal fitted parameter. The confidence intervals of all parameters are finite in both directions.