Isothermal titration calorimetry and surface plasmon resonance analysis using the dynamic approach.

Biophysical techniques such as isothermal titration calorimetry (ITC) and surface plasmon resonance (SPR) are routinely used to ascertain the global binding mechanisms of protein-protein or protein-ligand interaction. Recently, Dumas etal, have explicitly modelled the instrument response of the ligand dilution and analysed the ITC thermogram to obtain kinetic rate constants. Adopting a similar approach, we have integrated the dynamic instrument response with the binding mechanism to simulate the ITC profiles of equivalent and independent binding sites, equivalent and sequential binding sites and aggregating systems. The results were benchmarked against the standard commercial software Origin-ITC. Further, the experimental ITC chromatograms of 2'-CMP + RNASE and BH3I-1 + hBCLXL interactions were analysed and shown to be comparable with that of the conventional analysis. Dynamic approach was applied to simulate the SPR profiles of a two-state model, and could reproduce the experimental profile accurately.

Introduction

Biophysical techniques such as isothermal titration calorimetry (ITC), and surface plasmon resonance (SPR) are used to thermodynamically and kinetically characterize the binding mechanism of the protein-ligand or protein-protein interactions, respectively [[1], [2], [3], [4]]. ITC measures the heat released or absorbed during the protein-ligand interactions [[1], [2], [3],5,6], whereas, SPR measures the change in reflective angle of the incident light caused by surface waves called ‘plasmon polaritons’. Since plasmon polaritons, are sensitive to binding events occurring on a surface, it proportionally affects the angle of the reflected light [7,8].

In ITC experiments, the thermogram with asymmetric gaussian-like peaks are numerically integrated and normalized with respect to the titrated ligand to obtain ‘normalized delta heat’ data (NDH) [9,10]. Assuming a kinetic model, NDH data is then analysed to obtain the stoichiometry, binding equilibrium, and enthalpy constants. On the other hand, the thermogram (time domain data without integrating the peaks) can be directly analysed as shown by Dumas etal, to obtain the kinetic rate constants [11]. Interestingly, the delay in ligand dilution after each injection and the heat released or absorbed due to binding events were modelled as a first order ‘instrument response’. While doing so, the instrument response is considered independent of the kinetic binding mechanism. Here, in this article, we have incorporated both the instrument response and the binding mechanism within an integrated kinetic framework.

In a conventional SPR experiment, the response unit (RU) of the SPR sensogram exhibits, three distinct phases such as association, dissociation and regeneration [7]. These phases are either analysed separately/piecewise or in an integrated manner to obtain the kinetic rate constants. Dynamic approach, as applied to ITC can also be easily extended to model the SPR profiles of different binding mechanisms [[12], [13], [14], [15]]. Yet again the advantage of dynamic approach is that the instrument response can be seamlessly integrated within the kinetic framework of the binding mechanism thereby simplifying the data analysis. A representative of different kinds of binding mechanisms such as, single set of equivalent sites, two sets of equivalents sites with sequential binding mode, two sets of equivalents sites with parallel binding mode, and aggregation have been considered here [[16], [17], [18], [19], [20], [21], [22], [23], [24], [25], [26], [27]]. The simulation of ITC profiles for all these mechanisms were realised and found to be consistent with that of the previous reports [9]. Further, we have also analysed the experimental ITC data of 2′-CMP + RNASE and BH3I-1 + hBCLXL, and determined the kinetic and thermodynamic parameters of the binding process. The results were consistent with the earlier reports [28,29]. Dynamic approach based simulation of the SPR profiles for a single binding site also yielded accurate profiles similar to the reported experimental profiles.

Theory

Protein-ligand/Protein-protein binding mechanism can be as simple as a single site binding or as complex as multi step sequential binding. Previous modelling approaches have taken into account the instrument repsonse, independent of the kinetic mechanism; wherein, the correction for ligand dilution and heat detection were made outside of the dynamics of the binding mechanism. Firstly, a detailed comparison of different modelling approaches, such as lumped modelling, sequential kinetic modelling, and the current parallel kinetic modelling, using a simple single site binding mechanism, is provided, Sec (2.1) & SI (1.1–1.4); Fig 1A-H. Secondly, how dynamic approach can be applied to model complex cases such as a single set of equivalent binding sites Sec (2.2.1), two sets of equivalent sites with sequential Sec (2.2.2) or parallel binding modes SI (2.1), and aggregation SI (2.3) are presented, Fig (2). Extension of this approach to SPR for a single binding case, with and without ligand leakage is detailed in SI (3), Fig (4).

Fig. 1

Comparison of the thermogram and the NDH data obtained for a single binding site mechanism using four different approaches. (A,B) without instrument response; (C,D) with instrument response based on lumped modelling; (E,F) with instrument response based on kinetic modelling in a sequential manner; (G,H) with instrument response based on kinetic modelling in a parallel manner.

Fig. 2

Simulation of the ITC thermogram and its corresponding NDH data for different binding mechanisms (A,B) M equivalent single site binding; (C,D) M, N, two equivalent independent/parallel binding sites (FEOTF54); (E,F) M, N, two equivalent sequential binding sites (PROTDB). (G,H) M, N, O, R, four equivalent sequential binding sites (PERSSON). In the NDH plots of B, D, F, H, the open circle represents the NDH data points obtained independently through simulation based on algebraic model and the smooth line represents the NDH data obtained from integrating the simulated thermogram shown in A, C, E, G, respectively. The parameters used to simulate both algebraic and dynamic profiles (Table 1) were obtained through origin-ITC software by fitting the experimental data to appropriate models provided therein.

Fig. 4

Simulation of the SPR sensogram using dynamic approach for a single binding site mechanism. (A) Without any leakage of ligand during the dissociation phase (B) with leakage of ligand during the dissociation phase. The concentrations of the ligand used for each instance of the simulation is labelled above its respective traces.

Table 1

Parameters used to simulate the thermograms of different models using dynamic approach. These parameters are in turn obtained from fitting experimental data to equilibrium models available in Origin-ITC and are summarised in SI Table S5.

	RNAHH	FEOTF54	PROTDB	PERSSON
Model	One sites	Two sites (independent)	Sequential binding sites (2 sites)	Sequential binding sites (4 sites)
Parameters
Kinetic constants (k)	5.59 × 104 Hz(k_1)1.0 (Hz,k_−1)	1.18 ± 0.40 × 1010 Hz(k_1)1.0 Hz (k_−1)3.46 ± 0.91 × 106 Hz(k_2)1.0 Hz (k_−2)	4.13 × 107 Hz(k_1)1.0 Hz (k_−1)1.40 × 105 Hz(k_2)1.0 Hz (k_−2)	2.39 × 103 Hz(k_1)1.0 Hz (k_−1)114 Hz (k_2)1.0 Hz (k_−2)2.26 × 103 Hz(k_3)1.0 Hz (k_−3)22.3 Hz (k_4)1.0 Hz (k_−4)
Stoichiometric constant	1.02 (mindex)	1.06 (mindex)0.941 (nindex)	1.0 (mindex;fixed)1.0 (nindex;fixed)	1.0 (mindex;fixed)1.0 (nindex;fixed)1.0 (oindex;fixed)1.0 (rindex;fixed)
Thermodynamic constants(ΔHcals/mol)	-1.354 × 104 cal/mol (ΔH_bind)	767.3 cal/mol (ΔH_bind_1)1.203 × 104 cal/mol (ΔH_bind_2)	-8194 cal/mol (ΔH_bind_1)-3121 cal/mol (ΔH_bind_2)	314.8 cal/mol (ΔH_bind_1)4769 cal/mol (ΔH_bind_2)834.9 cal/mol (ΔH_bind_3)6962 cal/mol (ΔH_bind_4)
Injection, Cell parameters	20 (Inj no)4 μL(Inj vol)1345 μL (Cell vol)	17 (Inj no)5 μL (Inj vol)1411 μL (Cell vol)	20 (Inj no)4 μL(Inj vol)1320 μL (Cell vol)	52 (Inj no)2 μL (1st),5 μL, rest Inj vol)1320 μL (Cell vol)
Instrument response	0.33 Hz (k_lig)0.33 Hz (k_+h,k_−h)	0.33 Hz (k_lig)0.33 Hz (k_+h,k_−h)	0.33 Hz (k_lig)0.33 Hz (k_+h,k_−h)	0.33 Hz (k_lig)0.33 Hz (k_+h,k_−h)
Other parameters	651 μM (Prot)21160 μM (Lig)	31.4 μM (Prot)1560 μM (Lig)	494 μM (Prot)20700 μM (Lig)	6360 μM (Prot)315000 μM (Lig)
Integration parameters per peak in the thermogram	0–100 s (each inj)50 (data points)	0–100 s (each inj)50 (data points)	0–100 s (each inj)50 (data points)	0–100 s (each inj)50 (data points)

Dynamic model with instrument response

In this section, we outline the kinetic based modelling approach, which incorporates the instrument response within the binding mechanism. The kinetic mechanism can be proposed as follows:Where, , is the concentration of the injected ligand and , is the concentration of the ligand available for binding after dilution. Eqn (1) accounts for the instrument response of the ligand dilution from the center of the ‘cell’ into the bulk of the solution:

Eqn (2) accounts for the binding mechanism of the available ligand, with the protein, , to form the complex, and , , are the forward and reverse kinetic rate constants of the complex formation, respectively:

Eqn (3) accounts for the instrument response related to the heat released or absorbed caused by binding process. Here, the heat response rate () is related to the time delay , observed between the instance of heat generation and the instance of heat signal detected by the sensor/detector. As a simplification, we will always assume that, (), in the following discussion. The detected heat is represented in terms of ‘power’ in the thermogram.We have equated the ‘power’ pertaining to the binding process to be proportional to rate of change of , rather than itself. Here, is considered as a fictitious species of the protein complex , as detected or ‘seen’ by the instrument with an intrinsic response delay. The rate constant, , associated with the interconversion of to has the same unit as that of first order rate constants (). Following the same line of thought, Eqn (4) accounts for the instrument response associated with the heat of ligand dilution. Here, the ‘power’ due to dilution is made proportional to the rate of change of rather than itself; , is again a fictitious species associated with the instrument response of :

We have made an assumption here that the heat of dilution is primarily arising from the ligand rather than any buffer component or co-solvents. This assumption is valid because, in most of the experiments, the concentrations of the buffer and the co-solvents are maintained identical in both the cell and the injected ligand solutions, so as to avoid the heat of dilution due to buffers and co-solvents. The rate equations for the above kinetic mechanism can be written as Eqns (5), (6), (7), (8), (9), (10):

The above set of differential equations, Eqns (5), (6), (7), (8), (9), (10), can be solved simultaneously to obtain, , , at different time points. Given the time profile of , the time profile of the ‘power’ due to binding () can be calculated as follows, Eqn (11):Where, and are the concentrations of complex, at time periods, and , respectively. is the change in enthalpy and is the total volume of the cell. Given the time profile of , the time profile of the ‘power’ due to ligand dilution () can be calculated as follows Eqn (12):Where, and are concentrations of, at and , respectively.

The simulated ITC chromatogram, is nothing but the sum of ‘power’ due to binding and dilution of ligand, and . Comparison of four different approaches used to calculate the ITC time profile including the method described here SI (1.1–1.4) are shown in Fig (1 A-H). The parameters used to simulate these profiles are provided in SI Table (S1).

Complex binding mechanisms

In this section, we provide the explicit kinetic mechanism for complex binding cases such as,

equivalent single site binding (E.g. RNAHH; Fig. 2A and B).

, two sets of equivalent sites with sequential binding (E.g. PROTDB; Fig. 2E and F).

Other complex cases such as

, , two sets of equivalent sites and independent binding (E.g. FEOTF54; Fig. 2C and D).

, , , , four sets of equivalent and sequential binding (E.g. PERSSON; Fig. 2G and H).

Aggregation of proteins.

are explained in the supplementary SI (2.1–2.3).

equivalent Single binding site

In this mechanism, the ligand binds to a single site available in the protein, but the binding site can exist in equivalent conformational forms. At the outset the ligand ‘available’ for the binding process is modelled as shown in Eqn (1). The kinetic mechanism of equivalent sites binding can be outlined as in Eqn (13):

In the above representation of the complexes, to , the superscript index denotes, to , different equivalent sites; the subscript denotes the stoichiometry of ligand bound to the protein. The ‘power’ required to simulate the ITC profile, can be modelled as follows (Eqn (14)):Where, , , are the forward and reverse heat response rate constants of the heat detector, respectively, and are assumed to be equal to each other. In the above Eqns(14), the to , are the species, whose change in concentration generate the change in heat profile. While writing the mass balance for the total ligand or protein, we include only the species that are directly involved in the binding process (,, to ),and exclude the species concerned with the instrument response such as , to . For example, the mass balance of the ligand would be, and for protein it would be, . Additionally, we also consider a separate mass balance for the available and non-available ligand, that is independent of the binding process; . Based on the above proposed kinetic mechanism, the rate equations can be framed for each species and their time profile determined using numerical integration methods.

Before framing the rate equation, we make an assumption that the binding affinity for all the ‘’ equivalent binding sites are equal. This also implies that the concentration of all the bound species such as to are all equal. If we represent, , then the sum of all the bound complexes can be written as, . While writing the rate equations for the bound complex, only will be considered as a single entity, which represents the sum of all the equivalent bound forms.

The rate equation for the above kinetic mechanism, Eqns (1), (13), (14), can be written as Eqns (15), (16), (17), (18), (19), (20):Where, represents the total number of equivalent conformations, , represents the enthalpy due to binding, , represents the enthalpy due to dilution, , represents the total cell volume. The final ITC chromatogram () is the sum of ‘power’ due to binding ( as in Eqn (21)) and ‘power’ due to ligand dilution ( as in Eqn (22)) and is represented as, , where:

In the all the numerical simulations mentioned in this work we used the ‘difference’ form of the above equations Eqn(21, 22) rather than its differential form, as shown below Eqns (23), (24):

Simulations of ITC thermograms were carried out for a two-state system at three data acquisition frequencies, such as 10, 50 and 100 data points per injection. Though the ‘differential model’ and the difference model’ yielded an identical results at higher frequency (50 or 100 pts/injection), large discrepancy (from the intended result) was observed for the differential model compared to that of the difference equation model at lower frequency (10 pts/injection) SI (1.6). The simulation parameters and the corresponding ITC thermogram of a M equivalent binding site model is presented in (Table 1) and Fig 2A,B, respectively; matlab code to simulate thermogram using dynamic approach; and NDH data using algebraic method is provided in SI (1.7). ODE15() solver which is suitable for stiff problems (where significant changes occur within short period of time) was used for all the numerical integration mentioned in this work.

, , two sets of equivalent binding sites with sequential binding

In this mechanism, there are two binding sites available for the ligand to bind to the protein in a sequential manner. The first binding site can exist in equivalent conformational forms. The second binding site can exist in equivalent conformational forms. At the outset the ‘ligand available’, , for the binding process is modelled using Eqn (1). The kinetic mechanism of the first binding step is represented as below, Eqn (25):Where, , , are the forward and reverse kinetic rate constants of the first binding process, respectively. In the above representation of the complexes, to , the superscript index has two place holders separated by comma. The first place holder is alloted for the first binding site and the corresponding index represent one of the possible conformations or equivalent sites. The second place holder is alloted for the second binding site and the corresponding index represent one of the possible conformations or equivalent sites. An index value of , in any of the place holder indicates that the partcular binding process has not occurred yet. The subscript ‘1’, denotes that the stoichiometry of the ligand bound to the protein in all these complexes is ‘1’.

The second binding step which follows the first binding step, can be represented as follows, Eqn (26):

In the above representation of the complexes, to , the superscript index denotes 1 to n equivalent sites, for the second binding process; the subscript ‘2’ denotes that the stoichiometry of the ligand bound to the protein in all these complexes is ‘2’.

As explained in the case of equivalent single site binding, the concentrations of all the bound complexes after first binding steps can be represented as , and the sum of all the bound complexes for the first binding process can be written as, . Similarly, the concentrations of all the bound complexes after second binding step can be represented as , and the sum of all the bound complexes can be written as, . While writing the rate equations, for the bound complexes, only and will be considered, which represent the sum of all the equivalent bound complex forms corresponding to the first and second binding processes, respectively. The heat response equilibrium can be written as Eqn (27):

The rate equation for the above proposed kinetic mechanism can be written as, Eqns 28–35:

The final ITC chromatogram () is the sum of the ‘power’ due to binding (: Eqn (36) and : Eqn (37)) and ‘power’ due to ligand dilution (: Eqn (38)) and is represented as , where:

The simulation parameters and the corresponding ITC thermogram based on Eqns 28–38 is presented in (Table 1) and Fig 2E,F. As mentioned in M equivalent single site binding case, here too, we have used the ‘difference form’ of the differential Eqns 36–38, for all of the numerical simulations.

Results

Standard experimental data belonging to four different mechanisms such as (1) single site (RNAHH) (2) two sequential sites (PROTDB) (3) four sequential sites (PERSSON) and (4) two independent sites (FEOTF54) as provided in Origin-ITC software package were subjected to the conventional NDH based analysis [9,10]. The equilibrium constants and thermodynamic parameters obtained through optimization in Origin-ITC are tabulated in SI Table (SI 5). The same parameters were used to simulate the thermograms of RNAHH (Fig (2 A,B)), FEOTF54 (Fig (2C,D)), PROTDB (Fig (2 E,F)) or PERSSON (Fig (2 G,H)) using dynamic models such as (1) M equivalent single site, (2) M, N equivalent parallel/independent binding sites and (3) M, N (or M, N, O, R) equivalent sequential sites (2 or 4 sites) respectively, Fig (2), Table (1). The NDH data resulted from the integration of the simulated thermograms were identical to that of the NDH values obtained independently from conventional algebraic approach using the same set of fit parameters SI (1.7). Additionally, simulations for protein dimerization and octamer formation were also performed with the aggregation parameters as tabulated in SI Table (S2).

The experimental ITC data of 2′-CMP + RNASE and BH3I-1 + hBclXL were analysed using M equivalent single site model and M, N equivalent two site sequential model, respectively, Table 2 & Fig 3A,B. The model fitted well to the experimental data and the results were comparable to that of the previously reported values [28,29]. In our current analysis, because of large residual values in initial fits for BH3I-1 + hBCLXL, we let the injection volumes to be varied during the optimization. While analysing, 2′-CMP + RNASE and BH3I-1 + hBCLXL, we fixed the protein, ligand and cell volumes as constants. Probability density based (PDF) based global sensitivity analysis for all fit parameters suggested that the model is equally sensitive to all the parameters used in the model, SI(4) & SI Fig (S3) [30,31]. Fitting of a complex ITC model to an experimental data can yield several equivalent sets of optimized solutions for the parameters. In contrast to optimizing all the parameters as carried out here, some of the parameters can be experimentally determined and fixed during optimization to minimize fitting bias.

Table 2

Fit parameters for 2′-CMP + RNASE and BH3I-1 + h BCLXL based on dynamic approach.

	2′-CMP + RNASE	BH3I-1 + hBclXL
Model	M Equivalent single site model	M,N Equivalent two sequential model
Parameters
Kinetic rate constants (k)	1.15 ± 0.06 × 103 Hz (k_on)2.0 ± 1.0 × 102 Hz (k_off)	3.09 ± 0.72 × 103 Hz (k_on_1)45.0 ± 4.0 × 102 Hz (k_off_1)2.56 ± 1.12 × 103 Hz (k_on_2)1.39 ± 0.25 Hz (k_off_2)
Equilibrium constants (K_eq)	47.55 ± 3.66 × 103(K_1)	6.82 ± 1.67 × 103(K_1)1.84 ± 0.81 × 103(K_2)
Stoichiometric constant	1.36 ± 0.01 (mindex)	1.75 ± 0.26 (mindex)2.55 ± 0.54 (nindex)
Thermodynamic constants	-1.02 ± 0.00 × 104 cal/mol (ΔH_bind)-1.34 ± 0.37 × 102 cal/mol (ΔH_dil)	-2.32 ± 0.00 × 104 cal/mol (ΔH_bind_1)-1.81 ± 0.00 × 104 cal/mol (ΔH_bind_2)4.99 ± 1.14 × 103 cal/mol (ΔH_dil)
Instrument response	7.11 ± 1.29 Hz (k_lig)6.47 ± 1.29 Hz (k_h)	15.26 ± 0.06 × 103 Hz (k_lig)3.86 ± 0.12 × 103 Hz (k_h)
Corrections to the Injection vol		0.02 ± 0.03 μL3.81 ± 1.20 μL4.29 ± 2.06 μL2.09 ± 1.80 μL1.17 ± 1.70 μL0.44 ± 1.63 μL-0.08 ± 1.58 μL-0.51 ± 1.55 μL-0.92 ± 1.51 μL-1.34 ± 1.46 μL-1.76 ± 1.41 μL-2.14 ± 1.37 μL-2.45 ± 1.33 μL-2.77 ± 1.28 μL-3.06 ± 1.24 μL-3.28 ± 1.21 μL-3.51 ± 1.17 μL-3.67 ± 1.15 μL-3.72 ± 1.14 μL-3.76 ± 1.13 μL-3.82 ± 1.12 μL-3.81 ± 1.12 μL-3.80 ± 1.12 μL-3.81 ± 1.11 μL-3.58 ± 1.14 μL-3.57 ± 1.13 μL-3.35 ± 1.15 μL-3.15 ± 1.17 μL-2.79 ± 1.22 μL

Fig. 3

The experimental data and its model fit for (A) 2′-CMP + RNASE system using M equivalent single site binding, (B) BH3I-1 + hBCLXL using M, N, two sequential binding sites.

Simulation of SPR profile was carried out for a simple single site binding model (or two state model). The dynamic model was modified to account for the baseline anomalies often encountered in SPR profiles. By introducing a leakage factor in the ligand channel during dissociation, the concentration dependent residual changes in the baseline could be modelled accurately (SI Table (S3)) and Fig (4). The model used for ligand leakage () during dissociation consisted of a constant part (basal value: ) and an exponentially decaying function which is dependent on time (t); ; where, , is the leakage rate. The simulation profiles were consistent with that of the earlier reported experimental profiles [7,8,32].

Discussion

ITC thermogram is sensitive to various parameters such as protein concentration, ligand concentration, cell volume and injection volume [33]. The dilution effect of protein/ligand due to injected ligand, can be calculated as explained in SI (1.2.3). In real time situation, the dispersion of the highly concentrated ligand injected into the cell from the centre of the cell in to the bulk of the solution follows fick's law based diffusion formulated in the form of partial differential equation (PDE) Fig (5) [34,35]. In such a model, the concentration of ligand at each spatial location (X, Y, Z coordinates) varies with time, till the solution becomes homogeneous at equilibrium time point. This spatial and temporal dependence of the concentrations of the ligand and the protein can be simplified by taking into account the efficient stirring of the injector. If the mixing is effective, the spatial dependence can be ignored and the PDE reduces to an ordinary differential equation (ODE) with dependence only on time. Under such circumstances, the dilution can be considered to be instantaneous, implying that equilibrium concentration is achieved instantaneously for the bulk solution, after the injection of ligand solution.

Fig. 5

A comparison of the ligand dilution effect as addressed by ODE and PDE based dynamic modelling. The left and the right most figures represent the initial and final condition of the ligand concentrations immediately after injection and final equilibrium states, respectively. The darker shades represent higher concentration. In the upper scheme (ODE model), we assume an instantaneous mixing of sample being injected over discretized period of injection. Whereas, in the lower scheme (PDE model) we assume that the homogenization of the injected ligand is both time and spatial dependent.

Conventional analysis (or NDH analysis) of ITC data requires integration of each peaks in the thermogram and representing it as a function of protein to ligand ratio. In time domain analysis (dynamic analysis), the thermogram is directly fitted to the model as explained in Sec 2.2.1, 2.2.2, to obtain, thermodynamic and kinetic parameters. A closer inspection of the thermogram in Fig. 1 A,C,E,G, clearly shows that though the time profile of with and without instrument response (left panel) differ significantly, the integrated data (right panel) is identical. This suggests that the heat energy represented by the area under the curve is conserved across time profiles of varied instrument response. In other words, NDH analysis is insensitive to instrument response seen with ligand dilution or heat detection.

In dynamic analysis, when an experimental thermogram is fitted against a kinetic model, the initial few data points may not fit well to the model. This is primarily due to the fact that, only a fraction of the data points contribute to the asymmetric gaussian like thermogram profile, whereas, the rest of the data are baseline. Hence, during model fit, a weightage factor is calculated for each data point which is proportional to the difference of its value from the baseline. The model fit improved significantly, when the optimization was carried out with weightage factor. After optimization, probability density function (PDF) based sensitivity analysis was carried out on the fit parameters using SAFE toolbox [30,31]. Basically, two PDF's were calculated, namely, (1) conditional where one of the parameters is fixed and the rest of the parameters are varied and (2) unconditional, where all the parameters are varied. The difference between the conditional and unconditional PDF is a measure of sensitivity of that parameter to the given model. Based on this approach, all the optimized parameters were found to be invariably sensitive and significant to the models used in the analysis of RNASE and hBCLxl.

Dynamic analysis of 2′-CMP + RNASE experimental data yielded an equilibrium constant of 4.75 × 104 and an enthalpy of -1 × 104 cal/mol. The fit parameters were comparable to that of the previous NDH based analysis ( M−1 and cals/mol) [9]. BH3I-1 + hBCLXL on the other-hand was fitted to a , equivalent two site sequential binding site model. Previously, NDH based analysis carried out on the same data using a two-site sequential model with and , yielded, , M−1, which is comparable to the dynamic fit values of, , M−1 [28,29]. Similarly, the enthalpies from both the models were also comparable (NDH/(dynamic analysis):- , cals/mol and , ) cals/mol) [28,29]. Fitting of a complex ITC model to an experimental data can yield several equivalent sets of optimized solutions for the parameters. To increase the confidence level in the fit parameters, it is important to estimate or determine as many fit parameters as possible through other experiments, and treat these parameters as constants in the final model analysis.

Aggregation mechanisms are usually complex and difficult to model, due to the various binding modes available for the intra and inter molecular complex formation [25,27,36]. However, any complex aggregation mechanism can be modelled as a series of sequential binding mechanism, with several steps in between the monomer and the final aggregate polymer. In this study, we compared two aggregation cases, where a protein forms dimer (); and an octamer (). When the degree of polymerization for the protein increased from 2 to 8, distinct lag phase at the outset of the thermogram could be clearly seen [27]. In general, the ITC profile of the aggregation of a macromolecule in the presence of a ligand will exhibit an initial lag phase (flat line) correlated with the critical aggregation concentration (CAC), and is related to the non-cooperative binding. Subsequently, the co-operative phase steps in causing a rapid increase or decrease in Q value (heat) due to the intra-micellar aggregation. The higher order association between polymeric protein and polymeric ligand is characterized by the inter-micellar aggregation occurring at critical micellar concentration (CMC) [27]. Though our simulations are carried out only for the case of macromolecule aggregation, comparison could still be made with the first two phases of the experimental ITC profiles of SDS + PEG systems, exhibiting a characteristic lag and cooperative binding phases.

Conventional data analysis of SPR comprises of three distinct phases namely, association, dissociation and regeneration [7]. Since these three regimes are independent, these can be individually analysed using three separate analytical expressions. But dynamic approach integrate all these regimes seamlessly through a single set of ODE's, which, not only accounts for kinetics but also the instrument response due to ligand dilution and detector. A common feature that is often observed in the SPR profile is the concentration dependent residual baseline during dissociation/regeneration phase; which had been modelled here successfully by considering ligand leakage.

Conclusion

Here, we have explicitly, modelled the general binding mechanisms often encountered in ITC experiments such as independent, sequential and aggregation. Using dynamic approach, the instrument response due to ligand dilution and heat detection were incorporated within the kinetic framework so as to simulate the experimental profiles accurately. Furthermore, the experimental thermogram of 2′-CMP + RNASE and BH3I-1 + hBCLXL were analysed and the thermodynamic and kinetic parameters that are consistent with earlier studies were obtained. The same approach has been extended to simulate the SPR profiles of a single site binding mechanism.

Supporting Information contains explicit derivation of the models such as, thermogram without instrument response (Sec 1.1), with instrument response, (lumped model (Sec 1.2), kinetic sequential model (Sec 1.3) and kinetic parallel model (Sec 1.4). Equivalence of all these models are also detailed (Sec 1.5). Mathematical expressions for the dilution effect of protein and ligand are provided (Sec 1.2.3). Complex binding models such as M,N, two equivalent and independent site model (Sec 2.1), M,N,O,R four state sequential model (Sec 2.2), M aggregation model (Sec 2.3) are provided. Comparison of thermogram simulated from ‘difference’ and ‘differential’ equations (Sec 1.6) is provided. Matlab code to simulate the thermogram and NDH data using dynamic and algebraic method, respectively, is provided in Sec (1.7). Dynamic model for simulating SPR profile is also provided (Sec 3.0). Details of the global sensitivity analysis for the fit parameters is provided (Sec 4.0). Matlab codes [37] to simulate ITC/SPR profiles as mentioned in Fig. 1, Fig. 2, Fig. 3, Fig. 4, Fig. 5 and SI Fig (S1,S2) are provided as data in brief. Only the customized matlab code to carry out global sensitivity analysis is provided within data in brief. Complete "Safe toolbox" codes can be obtained from https://www.safetoolbox.info/register-for-download/.

Declaration of competing interest

The authors declare no conflict of interest.