Biomolecular Interaction Analysis Using an Optical Surface Plasmon Resonance Biosensor: The Marquardt Algorithm vs Newton Iteration Algorithm.

Kinetic analysis of biomolecular interactions are powerfully used to quantify the binding kinetic constants for the determination of a complex formed or dissociated within a given time span. Surface plasmon resonance biosensors provide an essential approach in the analysis of the biomolecular interactions including the interaction process of antigen-antibody and receptors-ligand. The binding affinity of the antibody to the antigen (or the receptor to the ligand) reflects the biological activities of the control antibodies (or receptors) and the corresponding immune signal responses in the pathologic process. Moreover, both the association rate and dissociation rate of the receptor to ligand are the substantial parameters for the study of signal transmission between cells. A number of experimental data may lead to complicated real-time curves that do not fit well to the kinetic model. This paper presented an analysis approach of biomolecular interactions established by utilizing the Marquardt algorithm. This algorithm was intensively considered to implement in the homemade bioanalyzer to perform the nonlinear curve-fitting of the association and disassociation process of the receptor to ligand. Compared with the results from the Newton iteration algorithm, it shows that the Marquardt algorithm does not only reduce the dependence of the initial value to avoid the divergence but also can greatly reduce the iterative regression times. The association and dissociation rate constants, ka, kd and the affinity parameters for the biomolecular interaction, KA, KD, were experimentally obtained 6.969×10(5) mL·g(-1)·s(-1), 0.00073 s(-1), 9.5466×10(8) mL·g(-1) and 1.0475×10(-9) g·mL(-1), respectively from the injection of the HBsAg solution with the concentration of 16 ng·mL(-1). The kinetic constants were evaluated distinctly by using the obtained data from the curve-fitting results.

Introduction

Kinetic analysis of biomolecular interactions that are affected by partial mass transfer is the most difficult task in the quantity of binding kinetic constants [1-2]. There are several approaches have been suggested to perform the analysis of biomolecular interactions based on the use of the kinetic model, differing in how data are selected [3]. The high quality kinetic data can be extracted from the curve where binding is closer to equilibrium. The kinetic analysis of biomolecular interaction can also be applied to detect the toxic molecules in food safety, environmental pollutants and life science [4-5]. In immunology, the binding strength between antibody and antigen reflects the biological activities of the control antibody and its immune response significance in the pathologic process [6-7]. Moreover, both the association rate and dissociation rate of antigen to antibody (receptor to ligand) are the important parameters for the study of signal transfer between cells [8-9]. The traditional methods involving the kinetic analysis of biomolecular interactions are mainly included ELISA (Enzyme-Linked Immunosorbent Assay), Equilibrium Dialysis, and Affinity Chromatography [10-12]. By comparison with these technologies, the most prominent characteristics of the surface plasmon resonance biosensor are the real-time monitor of the kinetic process without marking the biological molecules. In these experiments of biomolecular interactions using the traditional method, a limited binding partner with suitable spectroscopic properties such as fluorescent tags or other labels should be screened to give a useful fluorescence spectrum [13]. However, in recent advances, the development of optical surface plasmon resonance (SPR) biosensors for accurate kinetic and affinity analysis makes the biomolecular monitor powerful [14-16]. The ability of SPR biosensors to analyze biomolecular interactions in real time features the quantity of the affinity of ligand for its receptor and the kinetics parameters of the interaction [17]. Furthermore, SPR technology makes possible a detailed analysis of biomolecular interactions at the molecular level, as well as enabling the analysis of multimolecular complex assembly and function [18-19]. From the sensorgram obtained from the SPR biosensor, the relationship between RUs (response unit) and time(s) established in the process of association and dissociation of biomolecular interactions is a complicated nonlinear function. There is a number of instrumentation-based analysis software for use in obtaining kinetic constants. For example, sensorgrams may be figured out using one of several binding models provided with evaluation software from SensiQ and Autolab [20]. Another analysis approach is performed by using the powerful software OriginPro which can be applied to process the known data using the nonlinear curve fitting. Although the kinetic constants can be obtained easily from the known evaluation software, however, it can’t be embedded into the microcontroller in the design of the homemade bioanalyzer using the SPR biosensor due to intellectual property. The Newton iteration algorithm for the calculation of kinetic constants of biomolecular interactions has been described [21]. We have found that this method had a great dependence of the initial value. Moreover, the fitting results are likely to be divergent with different initial conditions. This paper proposes a powerful method to implement the kinetic data analysis for biomolecular interactions using the Marquardt algorithm based on Gauss-Newton algorithm. The pseudo first order kinetic model of biomolecular interaction was established firstly. Then, the data collected from the biomolecular interaction between the hepatitis B surface antigen (HBsAg) and the hepatitis B surface antibody (HBsAb) was obtained by the homemade SPR bioanalyzer [22]. Finally, we used this approach established by the Marquardt algorithm to perform the nonlinear curve-fitting for the calculation of the association and dissociation rate constants and the affinity constants. The results show that Marquardt algorithm does not only reduce the dependence of initial value to avoid the problem of data divergence but also greatly reduce the iterative regression times.

Materials and Methods

Materials

The three-channel Spreeta modules (TSPR1K23) manufactured with the gold slide bonded to the sensor modules were from Nomadics, Inc. (Stillwater, USA). Hepatitis B surface antibody (HBsAb) was purchased from Zhengzhou Biocell Antibody Centre (Henan, China). The HBsAb was stored frozen, and its standard solutions were prepared daily with phosphate buffer solution (PBS). An ELISA diagnostic kit for Hepatitis B surface antigen (HBsAg) was purchased from Shanghai Rongsheng Biotech Co., Ltd. (Shanghai, China). The standard HBsAg solutions were diluted with PBS (pH 7.4) and stored at 4°C. Double distilled water was used throughout the whole experiment. A 0.01M PBS (pH 7.4) was prepared by dissolving 0.24 g KH2PO4, 8.0 g NaCl, 1.44 g K2HPO4 and 0.2 g KCl in 1000 mL double distilled water.

Kinetic model

If a single ligand binds to the receptor in a 1:1 stoichiometric ratio to form the receptor-ligand complexes, the association process is described by considering two substances ligand L and receptor R, which were combined to emerge a complex LR [23-25]. In a practical reaction, both the association and dissociation processes occur simultaneously. For reversible associations and dissociations in a chemical equilibrium, it can be described by the following expression: where, k (mol•L-1•s-1) is the association rate constant used to describe the binding kinetic constant between ligand L and receptor R. The dissociation rate constant k (s-1) is the ratio of the concentration of the dissociated complex to the undissociated complex.

It is equally valid to write the rate equations as follows: where the brackets denote concentrations of the free R, free L and the concentrations of the complex [RL] at equilibrium. From this equation, it can be seen that dissociation rate k and association rate k for a given system can be determined any time. The concentrations of [R], [L], and [RL] are measured under equilibrium conditions.

The net rate reached approximately to zero when the equilibrium condition was formed. That is and k C C = k C , therefore, it can be already expressed as where, K and K are the equilibrium association and dissociation constants.

In the ligand binding process, two reactions take place as follows: (a) the total number of associations per unit time interval in a particular region is proportional to the total number of receptors involved, because they all can create a complex with the same probability [26-27]. The relationship among the amount of the complexes formed per unit time Ca, the instantaneous concentration of the free analyte C , and the concentration of free receptors C -C is expressed as (b) on the other hand, for each compound, there is certain probability that it will be dissociated into ligand L and receptor R within a unit time interval. This probability is the same for all compounds at the given conditions. The dissociation leads to a decrease of the compound concentrations proportional to its instantaneous value described as: where, C is the amount of the complex LR associated per unit time.

The rate of consumption of ligand L depends on both the concentration of ligand L and the concentration of receptor R. The chemical equilibrium Eq (1) can be expressed by the pseudo first order reaction rate equation (Kinetic equation) [28]. The corresponding differential equation is derived as follows: The instantaneous concentrations of complex LR can be indicated by the response values (R) of the SPR biosensor. Furthermore, the concentrations of unfree receptor R obtained at equilibrium are represented by R , the concentration of free receptor R is R -R, accordingly, the Eq (8) can be rearranged to: If the initial value R is 0 at the initial time t (t = 0), the value R can be solved from the Eq (9) using the Integral Transformation Method, which is written as the following expression illustrated at the arbitrary time t, where, k = k C +k

Then, we use the value of LR instead of , the Eq (10) can be expressed as, When the ligands combine with the receptors completely in the area of association, the dissociation process of compounds occurs. Therefore, in the dissociation process with the concentration of ligand L of 0, the Eq (9) can be rewritten in the following form.

For solving Eq (12) by Integral Transformation Methods, we get the values of RU from the experiment performed by the SPR biosensor, that is where, t is the initial time of dissociation, t is arbitrary time between the initial time and the end time, and RU is the response value of the SPR biosensor at time t . The affinity constants can be determined from the data obtained from the SPR biosensor at the steady response state during the association phase.

Now, assume y = RU, a = LR and b = k , the Eq (11) can be simplified as: The kinetic model of dissociation process (13) can be simplified to where, m = k . The value of m which is the dissociation rate constant calculated from the Marquardt algorithm was evaluated firstly. Then the association rate constant can be obtained in accordance with the expression b = k = k C +k . Accordingly, the values of affinity constant K and K are calculated respectively.

Establishment of Curve-Fitting Algorithms

Gauss-Newton Algorithm

For the kinetic model of association y = a(1-e ), the corresponding y were obtained from the experiment of the SPR biosensor. Once parameters a, b are obtained, the kinetic model of association for a particular biomolecular interaction can be formed successfully. In order to solve the equations, the initial value of a, b should be given, named a , b , respectively. The actual values of a, b were obtained from the following expressions:a = a 0+Δ1 and b = b 0+Δ2, where,Δ1,Δ2 represent the increments of a 0, b 0, respectively. Then the values of Δ1,Δ2 will be obtained from the following procedures. The function of y = a(1-e ) is expanded using the Taylor series at the point (a 0, b 0). The results will be expressed in Eq (16) by ignoring the quadratic term.

The residual value Q between experimental and theoretical value is obtained by utilizing least square method. The expression is shown as following, where, y(a , b ) is determined by the known a 0, b 0. Both and are the function of independent variable x. Moreover, x is the experimental result. Hence, Eq (14) can be simplified to the linear relationship on Δ1,Δ2 as follows, Rearrange above Eq (18) to the following Eq (19).

Correspondingly, where, the following substitution will be done.

So, both expressions (17) and expression (19) can be arranged to: where The expression (26) is the equation involving both unknown parameters Δ1,Δ2. Once both Δ1,Δ2 were obtained, both a 1, b 1 can be obtained according to the following expressions a 1 = a 0+Δ1, b 1 = b 0+Δ2. Here, a 0,b0 are replaced by a 1,b1. The iterative process may be continued to do until the criterion for convergence is satisfied (e.g. max |Δ1|< = ε 1 and max |Δ2|< = ε 2). However, this method is more dependent on the initial values. Obviously, different initial values can cause the iterative divergence. For this reason, the Marquardt algorithm was introduced.

Marquardt algorithm

Marquardt algorithm is somewhat similar to the Gauss-Newton algorithm. Both the initial values a 0, b 0 are also previous given and the nonlinear model was implemented using the Taylor series expansion at the point (a 0, b 0). Both Δ1,Δ2 were figured out by conducting the tangential line. The iterative process may be continued to do until the situation for convergence is satisfied. In the Marquardt algorithm, residual value of Q is calculated by the following expression. Where, Converted the above expression (30) into the matrix form, where, The procedure of Marquardt algorithm to solve parameters a, b is shown as following. (a) Initial values a 0, b 0 are given at first, and the corresponding initial value of squared residual Q 0 was figured out. (b) Use iterative method to determine the parameter d. The initial value of d is set to be 0.01. Then, it was substituted into the expression (31) to calculate the value of Δ1,Δ2. The parameters a 1, b 1 and Q 1 were obtained in the same approach mentioned above. Compare the value Q 1 with Q 0, if Q 1>Q 0, adjust the value d to repeat the above process until Q 1

Experimental Validation for the Biomolecular Interactions between HBsAb and HBsAg

Retention of binding activity is the most important consideration when immobilizing a biomolecule on the Au film of SPR biosensor and can be measured by comparing the relative binding responses recorded as RU. (a) Immobilizing the HBsAb on the surface plasmon resonance sensor surface, then the PBS (pH 7.4) was used to clean the sensor surface in order to obtain a highly stable response baseline. (b) Flowing of sample solution of HBsAg diluted to 100-fold through the microfludic cell (estimated concentration 16ng·mL-1), and monitoring the binding of HBsAg to immobilized HBsAb on the Au film deposited on the SPR biosensor, it causes a modification of the refractive index at the surface and results a change of the resonant angle. After a certain reaction time, an equilibrium plateau was reached, where there are no changes of signals produced with time. We injected the PBS to the microfludic cell to block the association of HBsAg-HBsAb compounds. At this time, the association between HBsAg and HBsAb was not existed, so the concentration of HBsAg was 0. (c) The HCL (pH 3.0) solution was used to remove all the HBsAg molecules to regenerate the SPR biosensor. Then, the buffer solution PBS was injected to restore the baseline again and a new cycle was beginning.

Results and Analysis

The parameters and the initial values can be set manually according to the experimental data in the SPR biomolecular interaction analysis software based on the Marquardt algorithm, which was designed by our research group. Different results from various initial values can be compared each other visually so that we can select the most ideal initial values to reduce errors. From the homemade SPR analysis software, the experimental data can be imported conveniently. The x-axis of the graph represents time (s), while the y-axis of the graph represents the signal responses indicated with RU, which was computed based on the following formula RU = (1.334-RIx) ×30000, where 1.334 was the refractive index of deionized water. RIx is the refractive index of an unknown sample, which can be measured by using the SPR biosensor in real-time and 30,000 is a pre-determined factor for increasing the sensitivity of the calculated responses [22]. The initial time was set to be 250s, and the time of association phase and dissociation phase were 251s and 38s. Fig 1 is the response curve for the sequential injection of HBsAg solution, indicating the association phase and dissociation phase of HBsAg and HBsAb. The dotted lines mark the injection of the HBsAg solution with RU (response unit) values.

Fig 1

Sensorgram showing the association and dissociation processes of biomolecular interaction between HBsAg and HbsAb.

The data marked with a triangle is obtained in the average of more than three sets of measurement results in RU. This sensorgram is showing that the HBsAg was binding on the specific HBsAb (association phase) starting from the injection point a and reaches an equilibrium after approximately 251s. From the dissociation starting point b, the dissociation phase was formed sequentially. The microfludic cell of this SPR bioanalyzer was kept at a constant temperature of 37°C.

The curve-fitting obtained from using both Newton Iteration algorithm and Marquardt algorithm was shown in Fig 2.

Fig 2

The fitting results obtained from both Newton Iteration algorithm and Marquardt algorithm.

The data marked with a triangle is obtained in the average of more than three sets of measurement results in RU, and the fitted curve was marked with a solid line. A. The curve-fitting using the Newton Iteration algorithm with the initial value of 0.0095, B. The curve-fitting using the Newton Iteration algorithm with the initial value of 0.011, C. The curve-fitting using the Marquardt algorithm with the initial value of 0.0095, D. The curve-fitting using the Marquardt algorithm with the initial value of 0.011.

In sub-figures of Fig 2, the dotted lines represent the original data collected from the data acquisition system of this homemade bioanalyzer designed by using the homemade SPR bioanalyzer [22], while the solid lines represent the nonlinear curve-fitting to match the affinity kinetic model. From Fig 2, the following conclusions can be obtained: (1) The fitting results obtained from the Newton Iteration algorithm is divergent to some non-proper initial values (Fig 2A and 2B). However, Marquardt algorithm can avoid this problem effectively; (2) In the same initial conditions, the results obtained from the Marquardt algorithm had much better than the results obtained from the Newton iteration algorithm (Fig 2C and 2D).

In this experiment, the introductions of damping coefficient in Marquardt algorithm can real-time amend the increments of parameters so that the possibility of divergence is greatly reduced. The experimental data was obtained from HBsAg biomolecules with concentration of 16ng·mL-1. The curve-fitting is shown in sub-Figured in Fig 2 and the a = 3358.23246, b = 0.01188 in the association phase, a = 3358.73684, m = -0.00073 in the dissociation phase were obtained, respectively. Hence, the R , k , k K A, and K are estimated according to the data narrated above (See Table 1).

Table 1

Kinetic constants of molecular interaction between HBsAg and HBsAb.

Fitting curves	Kinetic models	Kinetic constants
Association process	RU=LR_eq(1-e^-k_obt)	LReq = 3358.232, kob = 0.01188 s-1, R max = 3559.486, CL = 16 ng· mL-1
Dissociation process	RU=LR_eqe^-k_d(t_2-t_1)	kd = 0.00073 s-1, KD = 1.0475×10−9 g·mL-1, ka = 6.969×105 mL·g-1·s-1, KA = 9.5466×108 mL·g-1

Conclusions

Surface plasmon resonance biosensors are important tools in characterizing biomolecular interactions as well as understanding the biomolecular recognition membrane established on the surface of the Au film of the SPR biosensor. A number of applications using SPR biosensors have been found in food safety, environmental pollutants and life science, which were related to the experimental design and the calculation of kinetic constants involving the biomolecular interaction process between antigen and antibody or receptors and ligand. To understand clearly the interaction process, we confirms that the data does indeed obey the pseudo-first-order binding interaction model and validates the extracted kinetic and affinity constants. This article has established an approach based on Marquardt algorithm for the analysis of the biomolecular interaction using an optical surface plasmon resonance biosensor. The Marquardt algorithm was addressed experimentally to provide better understanding the results compared to the Newton iteration algorithm with reducing the possibility of divergence when the curve fitting was established. A global fitting of the dissociation rate constant (k ) was performed firstly. The next global fitting with ka and Rmax (fixed k as a constant) was sequential obtained. The association and dissociation rate constants ka and k were 6.969×105 mL·g-1·s-1 and 0.00073 s-1, giving an affinity constant (K ) of 1.0475×10−9 g·mL-1 from the HBsAg with concentration of 16ng·mL-1, respectively. With the careful curve-fitting, surface plasmon resonance biosensors may be applied to provide accurate kinetic constants.