TitrationAnalysis: a tool for high throughput binding kinetics data analysis for multiple label-free platforms

Mathematical model for tool development

The TitrationAnalysis tool was built upon Mathematica v12.0 and can be easily adapted for Mathematica v13.0. The package as well as example input and output files can be accessed at https://github.com/DukeCHSI/TitrationAnalysis and at https://zenodo.org/record/7998652³³.

The TitrationAnalysis tool uses Equation 1 and Equation 2 shown below to fit the sensorgrams to a 1:1 Langmuir binding model. The tool provides the option to handle non-regenerative titrations (alternatively known as single cycle kinetics) that do not include a regeneration step between cycles.

The linear equation for fitting association data and dissociation data are:

Here R_t is the response at time t. C_i is the molar concentration of analyte in cycle i, R_max is the maximal response feasible. k_a is the association rate constant, k_d is the dissociation rate constant and t_asso is the absolute time when association ends. In non-regenerative titration fitting, $t_0^i$ fits for the extrapolated time where the response is 0 for analyte cycle i; in regenerative titration fitting, $t_0^i$ becomes a fixed value t₀, corresponding to the absolute time when the association starts. In local R_max fitting, $R_{max}^i$ fits for R_max value for analyte cycle i, and becomes a non-local R_max in the case of global R_max fitting. $R_{shift}^i$ is optional and fits bulk shift at the start of the association. This bulk shift is typically due to a mismatch between the analyte buffer and the running buffer used for collecting baseline and dissociation data, and will therefore typically disappear when association ends. This causes a signal disconnect both at the beginning and at the end of the association phase. $R_{drift}^i$ is optional and accounts for quick change in signal at the beginning of dissociation, due to factors such as the loss of non-specifically bound analyte. To avoid over-parameterization, $R_{drift}^i$ term will be dropped if $R_{shift}^i$ term is included. In practice, Equation 2 is modified to Equation 3, which produces identical kinetics and R_max estimations and has more stable fitting performance than Equation 2. Table 1 summarizes the parameter details.

Standard error estimation

Standard errors for R_max, k_a and k_d estimations are calculated from “NonlinearModelFit”, the Mathematica module used for data fitting using Equation 1 and Equation 3. The parameter optimization was done through the minimization of sum of square error. The standard error for K_D estimation was calculated through error propagation using the standard errors of k_a and k_d through Equation 4:

where the symbol Δ before k_a, k_d and K_D represents the standard error of the corresponding value.

Implementation of the TitrationAnalysis tool

Figure 1 demonstrates the installation and execution of the TitrationAnalysis tool. The user will need to install the package under the name “KineticsToolkit” and use the command Get["KineticsToolkit`"] to activate the package inside Mathematica. Then the TitrationAnalysis tool can be implemented for the appropriate platform (Figure 1D). A series of pop-up windows will guide the user through data import and global settings before sequentially fitting sensorgrams and generating output files. More details on implementation can be found at https://github.com/DukeCHSI/TitrationAnalysis.

Figure 1. The TitrationAnalysis tool can batch process sensorgram fitting and automatically generate reports.

This figure shows an overall schematic of the installation and execution of the TitrationAnalysis tool. Panel A shows the installation of KineticsToolkit, the overall package containing TitrationAnalysis tool. Panel B shows the input commands to execute the TitrationAnalysis tool and the output at the end of TitrationAnalysis tool execution. Panel C shows an example of PDF report pages automatically generated after fitting analysis. Panel D shows the available modules that can be called within the TitrationAnalysis tool to import and analyze data collected on different instruments.

Label-free platforms adaptability of TitrationAnalysis

TitrationAnalysis tool is designed to directly import a large amount of reference subtracted data exported from commercial software provided for different label-free platforms, with no or minimal reformatting. The tool has the ability to handle data from three different instruments for measuring binding kinetics data (Figure 1): Carterra LSA for high-throughput SPR_i, Octet Red384 for high-throughput BLI, and Biacore T200 for standard throughput SPR.

Operation of the TitrationAnalysis tool

Here we provide a general overview of how the users may typically operate the TitrationAnalysis tool, as shown in Figure 2. The minimal system requirements for using the TitrationAnalysis tool is the same as those for using the Mathematica environment overall: Windows 10 or higher, 19 GB of disk space and 4 GB of RAM (https://support.wolfram.com/6479).

Figure 2. Flowchart of TitrationAnalysis tool data importing and results exporting.

Typically, the commercial software is capable of data reference subtraction, zero analyte concentration cycle (blank cycle) subtraction and data smoothing. User is expected to do the aforementioned steps as data pre-processing and export of the pre-processed data. For data obtained on Biacore T200 and Octet Red384, the blank cycle subtraction can be done during automatic fitting using TitrationAnalysis tool if data for a zero analyte concentration cycle is provided and therefore is optional during pre-processing. After pre-processing, the time points and their corresponding response values can be exported as data tables in various file formats.

When calling a specific instrument module, the user will be prompted to provide the pre-processed data and a corresponding sample information spreadsheet. The content of the spreadsheet will be described in detail in section titled “Sample information and analysis preference settings”. If the content of the spreadsheet matches what is included in the pre-processed data, the tool will automatically fit each sensorgram sequentially and generate report files after all fittings are done.

Report files include a PDF report and a standalone report of parameter estimates. Each page of the PDF report will correspond to one sensorgram series associated with a given ligand surface and include the fitted sensorgram overlaid with underlying data, fitting residuals, a dose-response plot (Response at the end of association phase versus log of analyte concentration) and a summary of parameter estimates. Alongside of the PDF report, an additional report in .csv format will also lay out the details of kinetics parameter estimates, and associated standard errors of R_max, k_a, k_d, and K_D. The standalone .csv report can be readily used to calculate the relative standard error of each kinetics parameter estimate, the averages of the estimates among replicates and the percent coefficient of variation among replicates of the same kinetics parameter.

Sample information and analysis preference settings

A user-prepared spreadsheet with sample information and analysis preference is to be provided so that the TitrationAnalysis tool can correctly import and analyze as well as export fitting results. The spreadsheet can be in .csv or single tab Excel format. The information and preferences that are expected to be included in the spreadsheet are summarized in Figure 3.

Figure 3. Summary of information and preferences to be provided for each sensorgram prior to automatic fitting.

Carterra LSA software is capable of simultaneously collecting titration data for up to 384 spots on a single sensor chip, and can have pre-processed data on all spots exported collectively as a single file. Biacore T200 software is capable of exporting all titration cycles from a specific channel as a single file. Octet Red384 software is capable of exporting data from each sensor as a single file. For Carterra LSA, the TitrationAnalysis tool requires the user to list sample information for all spots, with each spot appearing once, and choose what subset of ligands to be included in fitting. For Biacore T200 and Octet Red384, the user is only required to include relevant sensorgrams, and the same sensorgram can appear multiple times with varying analysis preferences.

TitrationAnalysis internal workflow

After matching the instrumentation with the user provided information spreadsheet, for each titration series sensorgram, the TitrationAnalysis tool extracts data points based on the sample locations user has provided. Then the following steps will be executed to prepare data for fitting:

Step 1 and 3 are depicted in Figure 4.

Figure 4. The TitratonAnalysis tool can do automatic baseline alignment and analyte concentration range down-selection before fitting.

Fitting process of a set of simulated regenerative titration data (A–D) and a set of simulated non-regenerative titration data (E–H) is shown. Panels A and E show the titration data prior to baseline alignment. The color of each titration cycle and the corresponding analyte concentration in nanomolar is shown in the legend. Panels B and F show the titration after automatic baseline alignment. Panels C and G show the automatic selection of a subset of the analyte concentrations that best approximate the dose response linear range, highlighted as points in red in dose response curves, matching kinetics traces in red in panels B and F. Panels D and H show the resulting fitted sensorgram overlaid on top of underlying data points. Titration data were simulated using Equation 1 and Equation 3 with k_a = 1×10⁵ (M^-1 s^-1), k_d = 1×10^-4 (s^-1) and R_max = 100 (RU).

The Mathematica module “NonlinearModelFit” is used to call the kinetics model and conduct fitting. Depending whether the user chooses to include bulk shift, to fit for global R_max or fit for regenerative cycles, the correct variation of kinetics model will be called for fitting.

Materials

CH31^31,34 and AE.A244 gp120^32,35,36 were produced by the Duke Human Vaccine Institute, Duke University. The transfection was done using 293 cells or CHO mammalian cells with plasmids for recombinant expression. The proteins were quality controlled for purity, including using SDS-PAGE, Western Blot and size exclusion chromatography.

Carterra LSA data collection

Kinetics titrations were performed using HC30M sensor chips (Carterra, Part# 4279) at 25°C. Aqueous solutions were delivered onto the sensor chip using the Carterra LSA microfluidic modules, including a 96-channel print-head (96PH) and a single flow cell (SFC).

Goat anti-Human IgG Fc antibody (Millipore, Cat# AP113) was first immobilized onto the chip through amine-coupling. Briefly, the chip surface was activated using 100 mM N-Hydroxysuccinimide (NHS) and 400 mM 1-Ethyl-3-(3-dimethylaminopropyl) carbodiimide hydrochloride (EDC) (Cytiva, Cat# BR100050, mixed 1:1:1 with 0.1 M MES buffer at pH 5.5) for 600 seconds. Then anti-Human IgG Fc (in 10 mM Sodium Acetate at pH 4.5) was immobilized onto the activated surface for 900 seconds at 50 µg/ml, followed by an injection of 1 M Ethanolamine-HCl at pH 8.5 for 600 seconds to quench unreactive esters. The chip was then exposed to two 30 seconds injections of 10 mM Glycine at pH 2.0. The anti-Human IgG Fc immobilization steps were done using SFC and 10 mM MES buffer at pH 5.5 with 0.01% Tween-20 as running buffer. CH31 was then captured by the anti-Human IgG Fc at 10 µg/ml for 600 seconds using the 96PH, with 1X HBSTE buffer (10 mM HEPES pH 7.4, 150 mM NaCl, 3 mM EDTA and 0.01% Tween-20) as running buffer and antibody diluent.

A two-fold dilution series of the antigen was prepared, with the top concentration for AE.A244 gp120 being 1µM. The antigen was then injected onto the chip surface from the lowest to the highest concentration sequentially without regeneration using SFC, preceded by 8 cycles of buffer injection for signal stabilization. For each concentration, the time-length for the data collection of baseline, association and dissociation was respectively 120 seconds, 300 seconds and 750 seconds. 1X HBSTE was used as titration running buffer and sample diluent.

The titration data collected were first pre-processed in the Kinetics (Carterra) software, including reference subtraction using empty spots on the sensor chip, blank subtraction and data smoothing. The data were analyzed within Kinetics software as well as exported and analyzed using the TitrationAnalysis tool.

Biacore T200 data collection

Kinetics titrations were performed using a CM5 sensor chip (Cytiva, Cat# BR100530) at 25°C. The activation of the carboxymethylated-dextran gold surface was achieved by injecting 200/50 mM EDC/NHS (Cytiva, Cat# BR100050) solution in ultrapure water pH 7.0 injected at 5 µL/min for 400 seconds. Following the activation step, a 50 µg/mL solution of anti-human IgG Fc (Millipore, Cat# AP113) in sodium acetate (NaOAc) pH 5.0 (Cytiva) was injected over the activated surface at 5 µL/min. Anti-human IgG Fc was injected in the sample channel for one injection of 120 seconds to reach ~7700 RU, and was injected in the reference channel for three injections of 200 seconds total to reach ~ 6900 RU. After covalent modification of the sensor surface, a quenching solution of ethanolamine pH 8.5 (Cytiva) was injected over the surface for 600 seconds to cap any residual active NHS esters.

PBS 1X pH 7.4 was used for the running buffer during titration. During the kinetics assay, one flow cell channel with only anti-human IgG Fc served as a reference channel to monitor and subtract binding responses due to non-specific interactions. 190–380 RU of CH31 at 5 µg/mL was captured onto the chip for each cycle at 5 µL/min for 60 seconds. The optimized capture of CH31 was followed by baseline monitoring for 60 seconds and the injection of AE.A244 gp120 for 180 seconds. Then a dissociation step was performed using an injection of running buffer for 600 seconds. Following the dissociation step, regeneration of the anti-human IgG Fc surface was performed using 1 injection of glycine•HCl pH 2.0 (Cytiva) at 30 µL/min for 40 seconds. The flow rate for association and dissociation was 30 µL/min.

The kinetics traces were reference subtracted using the responses of the reference channel in each cycle and blank subtracted using a zero-concentration cycle. Then the kinetics constants k_a, k_d and K_D values were determined using Biacore T200 evaluation software and TitrationAnalysis tool.

Octet Red384 data collection

BLI measurements were made using ForteBio biosensors (Fortebio - Sartorius). Both the Data Acquisition 12.0 and Data Analysis 12.0 software packages used were United States Food and Drug Administration’s (FDA) Title 21 Code of Federal Regulations (CFR) Part 11 (FDA Title 21 CFR Part 11) compliant versions. All data collection were performed at 25°C using settings of Standard Kinetics Acquisition rate (5.0 Hz, averaging by 20) at a sample plate shake speed of 1000 rpm. CH31 was loaded onto Anti-human IgG Fc Capture (AHC, Part# 18-5060) sensors with a Δλ = 0.5 nm loading threshold. The AHC sensors loaded with CH31 were then dipped into 1x kinetics buffer (ForteBio, Part# 18-1105) for 60 seconds to obtain baseline and then dipped into wells containing AE.A244 gp120 at different concentrations in 1X kinetics buffer to monitor antibody association. The dissociation step was monitored for 900 seconds by dipping Ab-bound sensors back into the wells used for baseline measurements to facilitate inter-step correction.

Antigen specific binding responses were obtained by double referencing; subtracting responses of blank AHC sensors tested in parallel and 1X kinetics buffer. The specific binding responses were fitted using ForteBio Data Analysis 12.0 software and TitrationAnalysis tool.

Titration data fitting

All sensorgrams were fitted using 1:1 binding model. For fitting of Biacore T200 and Octet Red384 data using TitrationAnalysis tool, data was thinned to one data point per second before fitting. During the fitting for data from Carterra LSA and Octet Red 384, signal shift at the beginning of dissociation was not fitted for when using either the commercial software (Carterra Kinetics software and Data Analysis 12.0) or the TitrationAnalysis tool. During the fitting for data from Biacore T200, signal shift at the beginning of dissociation was fitted for when using both the commercial software (Biacore T200 Evaluation Software) and the TitrationAnalysis tool.

Parameter estimates from TitrationAnalysis matched well with outputs of commercial software

In order to assess the quality of results generated by the TitrationAnalysis tool, we collected binding titration data of HIV-1 envelope glycoproteins AE.A244 gp120 binding to HIV-1 neutralizing monoclonal antibody (mAb) CH31. The binding titration data were collected on Carterra LSA, Biacore T200 and Octet Red384. Non-regenerative cycles were used when collecting data on Carterra LSA.

Figure 5 and Table 2 show the comparison of fitted sensorgram and parameter estimates. To shorten the fitting time, the data collected on Biacore T200 and Octet RED384 were thinned to 1 Hz (one data point per second) when fitting using the TitrationAnalysis tool. The data collection frequency of Carterra LSA was ~ 0.5 Hz (about 1 data point per 2 seconds).

Figure 5. Comparison of fitted sensorgrams obtained using commercial software with the TitrationAnalysis tool fitted sensorgrams.

Each panel from A to F shows the binding of AE.A244 gp120 to mAb CH31. Data collected on Biacore T200 are compared in A (Biacore T200 evaluation software) and B (TitrationAnalysis); data collected on high-throughput SPR_i platform Carterra LSA are compared in C (Carterra Kinetics software) and D (TitrationAnalysis); data collected on Octet RED384 are compared in E (ForteBio Data Analysis software) and F (TitrationAnalysis).

For Octet RED384 platform data, the kinetics estimates and the associated standard errors between the commercial software and the TitrationAnalysis tool are essentially indistinguishable. For the other two platforms, the estimates also closely resemble between commercial software and the TitrationAnalysis tool.

For Biacore T200 platform data, the kinetics estimates of the TitrationAnalysis tool showed less than 4% differences when compared to estimates from the commercial software, with the largest being k_a (3.2%). With the commercial software fitting data obtained at 10 Hz and the TitrationAnalysis tool fitting data obtained at 1 Hz, the standard errors from the TitrationAnalysis fitting only showed modest increase (2.13 – 3.22 fold) and were negligible compared to the estimates (%CV < 1%).

For Carterra LSA platform data, the kinetics estimates for AE.A244 binding showed less than 6% differences between the two fitting methods. Standard errors of k_d and K_D values associated with TitrationAnalysis tool fit were smaller than those associated with commercial software fit. The standard error of k_a from both fits were comparable.

Replicate measurements of the same interaction yielded similar kinetics estimates

Next, we assessed whether the TitrationAnalysis tool can be used to compare replicate measurements of the same interactions. Multiple replicate sensorgrams of AE.A244 gp120 binding to mAb CH31 were collected on the Carterra LSA platform and compared.

Figure 6 and Table 3 show the comparison of fitted sensorgram and parameter estimates among the replicates. There are small variations among the values, with the largest being 2.11 fold difference among the k_a values of AE.A244 binding to CH31. The level of standard errors is reproducible for the different replicates. Here, the TitrationAnalysis high throughput capacity is utilized, and all replicates of a given gp120 species was analyzed in a single run. During this fitting run using the optimized analyte concentration range, the average fitting time for each sensorgram was about 6 seconds.

Figure 6. Comparison of fitted sensorgrams for replicate measurements.

Panels A–D show fitted sensorgrams using TitrationAnalysis tool for replicates of AE.A244 gp120 binding to mAb CH31. All data were collected on Carterra LSA.

The fitting output is not affected by the change in users and machines

To assess whether the fitting output can be reproduced by multiple users implemented on different computers running different Mathematica versions, two users were asked to independently analyze the exact same four sets of Carterra titration data shown in sensorgrams in Figure 5 and Figure 6. The testing was done on two separate computers, with one user using Mathematica 12.2 and another using Mathematica 13.0. The fitting results were compared to the estimates done using Mathematica 12.0 shown in Table 2 and Table 3, and were shown to be highly reproducible (Table 4). The parameter estimates and the associated errors were indeed identical and independent of the specific computer and software version used.