INTRODUCTION

Top Abstract Introduction Results Discussion Appendix References

Until recently, obtaining reliable equilibrium and rate constants for interacting biomolecules was often difficult and time consuming. Optical biosensors now offer a rapid way to determine equilibrium and rate constants without the need to label the interacting biomolecules (Garland, 1996; Silin and Plant, 1997). Because of these features these instruments, such as BIACORE (Biacore, Uppsala, Sweden), have gained wide use and have been employed in the study of the interactions of a variety of biomolecules, including proteins, nucleic acids, lipids, and carbohydrates (reviewed in Szabo et al., 1995; Raghavan and Bjorkman, 1995; van der Merwe and Barclay, 1996; Myszka, 1997; Schuck, 1997a; Fivash et al., 1998).

In BIACORE instruments, one of the reactants is immobilized on a sensor chip. We will call it the receptor, in analogy to a receptor on a cell surface, although in the literature it is often referred to as the immobilized ligand. The other reactant, called the analyte, flows past the chip. The immobilization of the receptor is usually accomplished by coupling it to a thin dextran layer that extends ~100 nm out from the sensor surface, 0.2% of the height of the flow chamber. Detection of binding is based on the optical phenomenon of surface plasmon resonance (SPR). The SPR response is used to detect changes in the index of refraction caused by mass changes at the sensor surface. These changes are brought about by the binding of analyte to receptor (Jönsson et al., 1991; Malmqvist, 1993; Garland, 1996). Thus continuous monitoring of the SPR signal allows the kinetics of binding to be followed in real time. After binding, buffer alone may be introduced to monitor the dissociation kinetics.

As the preceding description suggests, the SPR signal in a typical experiment is not simply a report of the progress of the chemical processes of association and dissociation at the surface. Rather, the signal is the result of a combination of these chemical processes and the transport processes of diffusion and flow. Thus obtaining estimates of the intrinsic forward and reverse rate constants of the analyte-receptor reaction requires some model that will, when needed, account for transport. The model most often used assumes that after a brief transient, during which analyte is transported to the sensor surface, no correction is needed, because the free analyte concentration remains uniform in space and constant in time, kept so by the continuous influx of new analyte. We will refer to this model as the rapid mixing model because, after the brief transient, the binding at the sensor surface is the same as for a well-mixed system with constant analyte concentration.

Systematic errors in the estimates of the rate constants can arise if the rapid mixing model is not a reasonable description of the binding kinetics (Chaiken et al., 1992; Malmborg et al., 1992; Ito and Kurosawa, 1993; Morton et al., 1994; Wohlhueter et al., 1994). A number of authors have discussed this problem and presented conditions under which the assumption is expected to break down (Glaser, 1993; Karlsson et al. 1994; Schuck, 1996; Schuck and Minton, 1996; Yarmush et al., 1996; Christensen, 1997). There are two possible sources for the breakdown: 1) the transport (diffusion and flow) of analyte in solution to the dextran layer and 2) the transport (primarily diffusion) of analyte within the dextran layer. If one or both of these influence the kinetics of binding, then the free analyte concentration cannot be taken to be constant and uniform.

In this paper we consider only experiments in which the properties of the dextran layer do not influence the binding kinetics, i.e., where the binding can be treated as if the receptors are coupled directly to the sensor surface. BIACORE binding studies have been performed both with receptors coupled to dextran layers and with receptors coupled directly to the sensor surfaces (Karlsson and Fält, 1997; Parsons and Stockley, 1997). These experiments demonstrate that under appropriate experimental conditions (i.e., low surface capacity), the dextran layer has no significant effect on the binding kinetics.

To improve the analysis of bulk transport effects on the SPR signal, a simple model has been proposed that treats binding as a two-step process: transport of analyte to the sensor surface, followed by reaction of analyte with receptors on the surface (Myszka et al., 1997). This two-compartment model is attractive because it can be formulated in terms of a simple system of ordinary differential equations and because it is only slightly more complicated than the rapid mixing model. The two-compartment model has been used to describe several experimental systems studied with BIACORE instruments (Myszka et al., 1996, 1997; Morton and Myszka, 1998). This is encouraging, but uncertainty remains about how well the two-compartment model describes BIACORE binding kinetics. What is provided in this paper is a demonstration that the two-compartment model is sufficiently accurate to analyze BIACORE data. In particular, we show that the numerical parameter values obtained by fitting the two-compartment model to noisy SPR data are accurate estimates of the underlying physical parameter values.

The approach we use to demonstrate the validity of the two-compartment kinetic scheme for analyzing BIACORE binding data is as follows:

1. We present a computer model that can simulate BIACORE binding experiments involving monovalent analytes binding to monovalent receptors, i.e., that accounts for laminar flow and diffusion of the analyte, chemical reactions at the sensor surface, and geometry of the BIACORE. This computer model is similar to that developed by Glaser (1993), the difference being that we include diffusion in the flow direction.

2. We use the computer model to simulate BIACORE binding data, choosing parameter values that lead to transport-influenced binding kinetics as well as parameter values where such effects are negligible.

3. We fit the two-compartment binding model to the simulated data and show that we recover the original values of the rate constants (the values used to simulate the data).

4. We then add random noise at realistic levels to the simulations and demonstrate that good estimates of the input parameters are recovered when the two-compartment model is used to analyze the noisy simulated data.

	OBTAINING SIMULATED DATA SETS

The mathematical model

In BIACORE instruments, analyte is transported by diffusion and flow to the sensor surface, where it reacts with immobilized receptors. The flow channel of the BIACORE has a rectangular cross section (length l, height h, and width w). We let t be the time, and we introduce Cartesian coordinates with the origin at the inflow port, the x axis parallel to the direction of flow, and the y axis normal to the receptor-coated surface (illustrated in Fig. 1). (In BIACORE instruments, the sensor surface is actually on top of the flow cell. For convenience we draw it on the bottom.) Because the instrument's flow chamber is 10 times wider than it is high, variations in concentrations across the width of the chamber can be neglected.

View larger version (11K): [in this window] [in a new window]   FIGURE 1   Schematic of a BIACORE flow chamber, where the variation along the width of the chamber has been neglected and the sensor surface has been placed on the bottom rather than the top. For a standard instrument, h = 5 × 103 cm, l = 0.24 cm, and the width equals 5 × 102 cm. The figure is not to scale. The true length is 48 times longer than the height.

Along the flow chamber, laminar flow is fully developed essentially over its entire length (Brody et al., 1996). The velocity profile is therefore parabolic, equal to zero at the top (y = h) and bottom (y = 0) boundaries and rising to a maximum, v_c, in the center, i.e., the velocity v(y) at a height y above the sensor surface is (Batchelor, 1967)

(1)

(Because Eq. 1 is sometimes written in terms of the average velocity, , we note that  = 2v_c/3.)

The dependent variables of the model are the bulk concentration of free analyte, C(t, x, y), and the surface density of bound analyte, B(t, x) (with units of mass/volume and mass/area, respectively). The net analyte flux due to convection and diffusion is

(2)

where J is the flux vector, D is the analyte diffusion coefficient, and i and j are the unit vectors in the x and y directions. The component of J normal to any surface, integrated over that surface, equals the flux through the surface, e.g., the number of particles per second passing through the surface. From conservation of mass (the continuity equation), the following partial differential equation (PDE) is obtained:

(3)

The boundary conditions for Eq. 3 are as follows. Because the top boundary is impenetrable and unreactive, the mass flux of analyte must vanish at this boundary:

(4a)

In contrast, the analyte mass flux at the sensor surface should equal the time rate of change in the amount bound at the surface, i.e.,

(4b)

where C(t, x, 0) is the free analyte concentration at position x just above the sensor surface. R(t, x) = R_T  B(t, x) is the free surface receptor concentration at x. R_T is the total receptor concentration and is constant with respect to both position and time. The rate constants of the reaction are k_a and k_d. Note that these are the intrinsic microscopic forward and reverse rate constants for binding of the analyte to the receptor.

At the entry port to the flow cell (x = 0), the concentration is constant and equal to the injection concentration, C_T, i.e.,

(4c)

At the end of the sensor, where analyte exits (x = l), we adopt a simple continuation condition that assumes the exit of analyte is due entirely to flow:

(4d)

Hence they leave the computational domain and have a negligible effect on the dynamic processes happening inside the instrument.

It is useful to write Eq. 3 and the boundary conditions, Eqs. 4a-4d, in nondimensional form. To do this we introduce the following nondimensional time and lengths, concentrations, and rate constants:

(5a)

(5b)

(5c)

and nondimensional parameters:

(5d)

Equations 3 and 4a-4c become

(6)

(7a)

(7b)

(7c)

(7d)

If diffusion in the x direction is ignored (i.e., the ² term in Eq. 6 is set equal to zero), the model reduces to what has now become a standard description for the BIACORE flow cell (Lok et al., 1983; Glaser, 1993; Christensen, 1997). Yarmush et al. (1996) have also used this approximation, but in the context of a nonnegligible dextran layer. For a standard flow cell,  = 0.021. This suggests that such an approximation is justified, although for y < , diffusion in the x direction will dominate flow. We have not made this approximation, because we are using the model as an absolute standard for generating simulation data, and therefore it is desirable to include as much detail as possible. Because we are solving Eq. 6 numerically, there is no necessity to do otherwise.

The numerical method

We have solved Eq. 6 by using a general, conservative, and well-tested hydrodynamics transport code based on the method of finite elements (Dembo, 1994a,b). Shown in Fig. 2 is a small computational grid (8 × 8) with 49 node points and 32 surface elements. In the x direction the length is divided uniformly into eight parts. In the y direction the height first is divided into four equal parts. The top and bottom sections are then divided in half, so there are six sections in all. Finally, the new top and bottom sections are divided in half, forming the eight sections shown in Fig. 2. The height of the bottom and top sections is h/16. Refinements of this grid are generated by subdividing each of the quadrilaterals of the course mesh into four equal parts. For example, the top and bottom elements of a 64 × 64 grid are h/128 in height and l/64 in length. In the BIACORE flow cell, except at very short times after the initiation of binding and dissociation, the concentration of analyte is uniform far from the sensor surface. Thus there is no reason for generating a fine grid near the top surface. This was done only because the code we adapted required a grid that was symmetrical about h/2. We have since developed a code for which this is not required (unpublished result).

View larger version (13K): [in this window] [in a new window]   FIGURE 2   An 8 × 8 computational grid consisting of 49 node points and 32 boundary points.

In trial simulations we compared the time course for the average concentration of bound analyte, calculated using a 32 × 32 and a 64 × 64 grid. The average concentration of the bound ligand was calculated by summing the bound ligand concentration over each bottom surface element and dividing the sum by the number of bottom surface elements, n + 1 for an n × n grid. There was no significant difference between the values obtained using the two grids. Nevertheless, in the Results section, where we test the two-compartment model's ability to recover accurate parameter values, all simulations were done on a 64 × 64 grid. For this grid the first set of grid points above the sensor surface is at y = h/128  400 nm, or about four times the thickness of the dextran layer for a standard CM5 sensor surface (Stenberg et al., 1991).

Simulations

To illustrate the information that can be obtained from solving Eq. 6 numerically, subject to the boundary conditions in Eqs. 7a-7d, we simulate a BIACORE binding study using parameter values previously determined for interleukin-2 (IL-2), flowing past, and interacting with, its immobilized low-affinity receptor, IL-2R (Myszka et al., 1996). In the simulation we took the association rate constant k_a = 8 × 10⁶ M¹ s¹ and the dissociation rate constant k_d = 0.2 s¹. Analyte was injected for 50 s, after which buffer was introduced and dissociation proceeded for an additional 50 s. Fig. 3 shows the average concentration of bound analyte and the average free analyte concentration just above the sensor surface plotted as a function of time. We see from Fig. 3 that the free analyte concentration near the sensor surface is not constant during the association phase or zero during the dissociation phase, as is often assumed. Fig. 4 presents "snapshots" of the free analyte concentration, as a function of position throughout the flow cell, at four successive time points (two during the association phase and two during the dissociation phase). Note that at 25 s and 75 s there is substantial variation in the free analyte concentration along the sensor surface. As a consequence, the bound analyte concentration will also be nonuniform, as illustrated in Fig. 5. In BIACORE instruments the SPR detector is centered in the middle of the flow cell and monitors ~75% of the total flow path. Thus the output is a measure of the average bound analyte over this region.

View larger version (23K): [in this window] [in a new window]   FIGURE 3   Predicted time course for the average concentration of bound analyte () and average free analyte concentration at the sensor surface (- - -). In the simulation, 25 nM analyte was injected from 0 to 50 s (CT = 25 nM), then the injected concentration was set to zero (CT = 0) and dissociation was followed for 50 s. vc = 5 cm/s, corresponding to a flow rate of 50 µl/min, and the analyte diffusion coefficient D = 1 × 106 cm2/s. The immobilized receptor concentration RT = 7.5 × 1011 receptors/cm2 = 1.25 nM cm. (Because 1 RU = 1010 g/cm2, when all sites are filled with IL-2, MW = 14,000, there is an increase of 175 RU over baseline.) A standard flow cell was assumed, with h = 5 × 103 cm and l = 0.24 cm. ka = 8 × 106 M1 s1 and kd = 0.2 s1, as reported in the literature for IL-2 binding to IL-2R (Myszka et al., 1996). A 32 × 32 grid was used.

View larger version (76K): [in this window] [in a new window]   FIGURE 4   Predicted free analyte concentration in the BIACORE flow cell for the simulation parameters given in Fig. 3. Different colors indicate different free analyte concentrations in 0.5 nM increments. The contours separating two colors are lines of fixed concentration. The value of the concentration at any point on the grid is obtained by interpolating between the calculated values of the two nodes the grid line connects. Note that in the association phase (t = 25 s and 50 s) the contours separating the two colors nearest the top surface are not smooth. This indicates that in this region, more node points are required to improve accuracy. This simulation was done on a 32 × 32 grid. Despite the apparent lack of spatial resolution far from the boundary, when we compared calculations using a 32 × 32 and a 64 × 64 grid, the quantity we are most interested in, the average bound analyte concentration, , changed only slightly. Between 10 s and 90 s, the difference was always less than 1%.

View larger version (16K): [in this window] [in a new window]   FIGURE 5   Spatial variation in the bound analyte concentration, B(t, x). Shown are time courses for the bound analyte concentration in RU at a distance of l/4 (), l/2 (- · -), and 3l/4 (·····) along the sensor surface. The simulation parameters used are the same as in Fig. 3. The results demonstrate that the bound analyte concentration varies with the distance from the entrance to the flow cell.

The parameter values used in Figs. 3-5 were chosen to demonstrate how transport affects the kinetics of binding and dissociation in a BIACORE instrument. In Fig. 6 we show that as the association rate is decreased, transport effects become negligible, i.e., the concentration of free analyte becomes constant in time, and the rapid mixing model becomes a good description of the binding.

View larger version (15K): [in this window] [in a new window]   FIGURE 6   The time course of the average free analyte concentration just above the sensor surface, for different values of the forward rate constant, ka. Shown are time courses for the average free analyte concentration for ka = 8 × 104 M1 s1 (), 8 × 105 M1 s1 (- · -), and 8 × 106 M1 s1 (·····). All other simulation parameters are the same as in Fig. 3.

In the next section we discuss a proposed method of analyzing BIACORE binding data that are influenced by mass transport. To test this method we simulate additional BIACORE binding experiments under a variety of conditions and see if the method can recover the kinetic parameters used in the simulations.

	ANALYZING MASS TRANSPORT-INFLUENCED BINDING DATA

A two-compartment model

Obtaining kinetic parameters from binding data requires fitting a model to the data. At present, it is impractical to do this by using the mathematical description of the previous section, a PDE (Eq. 6), and its boundary conditions, Eqs. 7a-7d. A simpler mathematical model is required, one that can accurately calculate the average concentration of bound analyte at the sensor surface. This model need do nothing more. For example, it does not have to describe the spatial variation of the analyte concentration throughout the flow cell or of the bound analyte over the sensor surface.

A two-compartment model, consisting of a set of coupled ordinary differential equations, has been used to analyze a variety of BIACORE binding experiments that have been shown to be influenced by mass transport (Myszka et al., 1996, 1997; Morton and Myszka, 1998). First we review this model, and then in the following section we test it, to see if accurate rate constants can be recovered from simulated BIACORE experiments.

The model handles the variation in analyte concentration in the BIACORE flow cell (see Fig. 4) by dividing the flow chamber into two compartments as shown in Fig. 7. Within each compartment the concentrations are uniform in space but may change in time. These concentrations are averages over the length of the BIACORE flow cell. It is assumed that the analyte concentration in the outer compartment is constant and equal to the injection concentration, C_T. The model ignores the brief transition periods that occur when binding or dissociation is initiated and the concentration in the outer compartment rises or falls to C_T. At the typical flow rates used in kinetic experiments (50-100 µl/min, corresponding to v_c = 5  10 cm/s), the time needed for an analyte molecule that is not near the walls of the flow cell to travel 0.24 cm, the length of the flow chamber, is less than 0.05 s. The concentration in the inner compartment, C, changes because analyte is transported between compartments and because analyte binds to, and dissociates from, immobilized receptors on the sensor surface. If V_i is the volume of the inner compartment, S is the surface area of the flow chamber, and k_M is the transport coefficient describing diffusive movement of analyte between the compartments, then

(8a)

(8b)

The initial conditions for the binding phase of the experiment are that at t = 0, C_T equals the injection concentration, C = 0 and B = 0. The initial conditions for the dissociation phase are that at the time of dissociation, t_diss, C_T = 0, and C and B equal the values they attained at t = t_diss. If the association phase is sufficiently long so that a steady state is reached, then at t = t_diss, C = C_T and B = KC_TR_T/(1 + KC_T), but this is not required.

View larger version (10K): [in this window] [in a new window]   FIGURE 7   A two-compartment model of the BIACORE flow cell. In the binding phase the average analyte concentration near the sensor surface, C, equals zero initially and grows with time until it equals the injection concentration CT or dissociation is initiated and it decays to zero. The analyte concentration far from the sensor surface remains constant and equal to CT during binding, and constant and equal to zero during dissociation. B and R are the bound and free receptor concentrations, averaged over the sensor surface.

The rapid mixing model assumes that after a brief transient, the free analyte concentration is uniform throughout the flow chamber. This model is described by Eq. 8b with C = C_T.

As discussed in the Appendix, k_M is the diffusion-limited forward rate constant, averaged over the sensor surface. Transport effects will influence the kinetics of binding when the reaction rate is fast compared to the rate of transport, i.e., when k_aR_T  k_M. (To evaluate this inequality, the units must be consistent, so for example, if k_a is in M¹ s¹, then R_T should be in M cm and k_M in cm/s.) As discussed by Lok et al. (1983) and Sjölander and Urbaniczky (1991), it can be shown that, to a good approximation,

(9)

Because the model does not specify the height of the inner compartment, there are four unknown parameters that enter Eqs. 8a and 8b, k_a, k_d, k_M, and h_i = V_i/S. However, for the range of parameter values that describe BIACORE experiments, the solutions to Eqs. 8a and 8b are insensitive to the value of h_i. It can be shown (see Appendix and Fig. 3) that at the sensor surface, once binding or dissociation is initiated, C changes rapidly over a short period of time and slowly thereafter. Furthermore, while this rapid change in C is occurring, there is a negligible change in B. This behavior means that a quasi-steady-state approximation can be made where we set dC/dt = 0 (Segel and Slemrod, 1989). When this is done, we see from Eq. 8a that h_i drops out of the equations, which explains why the choice of h_i does not affect the solution. Once the quasi-equilibrium approximation is made, the model is identical to the often discussed effective rate constant model (Karlsson et al., 1994).

Returning to Eqs. 8a and 8b, we define  = B/h_i,  = R/h_i, _T = R_T/h_i, and _M = k_M/h_i. In terms of these quantities, Eqs. 8a and 8b take the familiar form

(10a)

(10b)

When these equations are used to analyze binding data, some care must be taken in the choice of units. If we measure R, R_T, and B in RU; C and C_T in M; k_a in M¹ s¹; k_d in s¹; and k_M in cm/s, this is equivalent to setting h_i = 1 RU/M. (Because the data are insensitive to the value of h_i, we are free to choose any value for it we wish. Taking h_i = 1 RU/M avoids having to divide the data values by h_i before fitting the data. It is a convenient choice having nothing to do with the physical height of the inner compartment.) To obtain k_M in cm/s from _M in 1/s, the parameter determined from the least-squares fit of the data, one must multiply by 10⁷ and divide by the molecular weight of the analyte, i.e.,

(11)

where we have used the conversion factor 1RU = 10¹⁰ g/cm².

Data fitting approach

We briefly summarize the data fitting approach we employ, because it has been discussed in detail elsewhere (Myszka, 1997; Myszka et al., 1997; Morton and Myszka, 1998). The rate constants k_a and k_d, _M (related to the transport coefficient k_M by Eq. 11), and the surface receptor density R_T are taken to be free parameters. The best fit values of these parameters are determined from nonlinear least-squares fitting of the data, where the calculated values of the average bound analyte concentration for each time point are determined by solving Eqs. 10a and 10b numerically (Morton et al., 1995). (The software program used in this analysis (CLAMP) may be downloaded for no charge at the following web site: http://www.hci.utah.edu/cores/biacore/docs/clamp.html.)

To obtain accurate values for the kinetic rate constants and the transport coefficient, the data that are being fit must be sensitive to these parameters. Recall that as the receptor density on the sensor surface is increased, the binding reaction at the sensor surface speeds up, and the binding kinetics become transport limited. As a result, at high surface receptor densities, the binding kinetics are dependent on k_M. In the other limit, as the receptor density is reduced, transport effects decrease and these data become insensitive to k_M and strongly dependent on the reaction rate constants. To try to ensure that the data depend on both the transport and the reaction rate constants, we simultaneously fit data sets for different surface receptor densities. The procedure is sometimes called global analysis and has been extensively used to analyze BIACORE binding data (Morton et al., 1995; Myszka et al., 1996, 1997; Morton and Myszka, 1998).

	RESULTS

Top Abstract Introduction Results Discussion Appendix References

In this section we report the results of fitting the two-compartment model to simulated data sets generated under a variety of conditions. Before considering the effects of noisy data on parameter estimation, we see how well the model does at recovering parameter values in the absence of noise.

We start by analyzing data where the reactions are not influenced by mass transport. Fig. 8 illustrates two sets of binding experiments simulated with identical rate constants (k_a = 5 × 10⁴ M¹ s¹ and k_d = 0.04 s¹) and flow rates (100 µl/min), but with different receptor densities (R_T1 = 50 RU and R_T2 = 15 RU). Note that when the reaction is not influenced by transport, for the same analyte concentration, the binding responses are identical for the two surfaces when normalized with respect to the maximum response, i.e., a plot of the fraction of receptor sites bound versus time is independent of R_T. Under these conditions, the rapid mixing model (Eq. 8b, with C = C_T) gives an excellent fit to the data (not shown) and returns correct parameter values, as shown in Table 1A, i.e., the best fit values of the parameters match the values used to simulate the data. The standard deviation of the residuals equals 0.054 RU, which is small but nonzero. Because the ODE models used to fit the data are always approximations to the PDE model used to simulate the data, even in the absence of noise, the standard deviation of the residuals is nonzero. As shown in Table 1B, applying the two-compartment model to these same data sets gives essentially the same results. The additional parameter, _M, has little effect on the quality of the fit (std residuals = 0.053 RU), the returned parameter values, standard deviations, or the correlation coefficients. However, the value returned for _M of 2.9 × 10¹² s¹, which from Eq. 11 corresponds to k_M = 9.7 cm/s for MW = 3 × 10⁴ Da, is grossly in error. From Eq. 9, the predicted value for k_M is 2.6 × 10³ cm/s. This error is not surprising. The association rate is so slow that the binding is unaffected by transport, and so the data contain no information for the determination of k_M.

View larger version (17K): [in this window] [in a new window]   FIGURE 8   Binding response simulated for a reaction occurring under conditions where the binding is reaction limited. The parameter values used in the simulations were ka = 5 × 104 M1 s1, kd = 0.04 s1, D = 1 × 106 cm2/s, and vc = 10 cm/s. The surface receptor concentrations were 0.167 nM cm (a) and 0.05 nM cm (b), which, for an analyte with MW = 3 × 104, correspond to a maximum response of 50 RU and 15 RU, respectively. The ligand concentrations used in the simulation were 8000, 2670, 890, 297, and 99 nM. The results of the fit are given in Table 1.

View this table: [in this window] [in a new window]   TABLE 1   Parameter estimates, standard deviations, and correlations for the fits of the rapid mixing model and the two-compartment model to data simulated in Fig. 8 with ka = 5 × 104 M1 s1 and kd = 4 × 102 s1

Next we carried out the same fitting procedures on data simulated with a 10-fold higher association rate constant (k_a = 5 × 10⁵ M¹ s¹). All other parameters were the same. Again a good fit was obtained with the rapid mixing model (std residuals = 0.062 RU), yielding parameter values that are close to the input values (Table 2A). As shown in Table 2B, using a two-compartment model improves the quality of the fit somewhat (std residuals = 0.055 RU). Furthermore, the best fit value _M = 8.9 × 10⁸ s¹, which corresponds to k_M = 3.0 × 10³ cm/s, is greatly improved, being within 13% of the expected value. The improvements in the residual standard deviation, rate constant estimates, and the value of k_M suggest that under these conditions, transport is beginning to influence the binding kinetics.

View this table: [in this window] [in a new window]   TABLE 2   Parameter estimates, standard deviations, and correlations for the fits of the rapid mixing model and the two-compartment model to data simulated with ka = 5 × 105 M1 s1 and kd = 4 × 102 s1

Fig. 9 illustrates what happens to the shapes of the binding responses when the association rate constant is increased to 1 × 10⁸ M¹ s¹ and the data become strongly influenced by transport. (Note that the off-rate constant was increased to 0.08 s¹, so that the shapes of the progress curves would remain similar.) Without doing any data fitting, one can see (compare Fig. 9 c with 9 d) that the binding is transport limited because, unlike in Fig. 8, the shapes of the binding curves for the two surfaces are not identical. Now equilibrium is approached more rapidly on the surface with the lower receptor density. When the binding is influenced by transport, the fraction of receptor sites bound is no longer independent of R_T. Receptors no longer act independently. In the association phase they compete for analyte, whereas in the dissociation phase rebinding occurs.

View larger version (38K): [in this window] [in a new window]   FIGURE 9   Binding response simulated for a reaction occurring under conditions where the binding is transport limited. The parameter values used in the simulations were ka = 1 × 108 M1 s1 and kd = 0.08 s1. All other parameters were the same as in Fig. 8. The surface receptor concentrations correspond to a maximum response of 50 RU (a, c) and 15 RU (b, d) for a 3 × 104 Da analyte. The analyte concentrations used in the simulation were 8, 2.67, 0.89, 0.297, and 0.099 nM. Shown are simultaneous fits (red lines) of the rapid mixing model (a, b) and the two-compartment model (c, d) to the simulations. The results of the fits are given in Table 3.

The fit of the rapid mixing model to the simulated data (Fig. 9, a and b) is poor (residual std = 1.52 RU), whereas the fit of the two-compartment model (Fig. 9, c and d) is excellent (residual std = 0.07 RU). Furthermore, in Table 3A we see that the two-compartment model returns parameter values that are close to those used in the simulations. The best fit value _M = 7.72 × 10⁸ RU s¹, or equivalently, k_M = 2.6 × 10³ cm/s, is within 1% of the predicted value. These results demonstrate that in the absence of noise, the two-compartment model returns accurate estimates for both the intrinsic reaction rate constants and the transport coefficient when the binding kinetics are influenced by transport.

View this table: [in this window] [in a new window]   TABLE 3   Parameter estimates, standard deviations, and correlations for the fits of the two-compartment model to data simulated with ka = 1 × 108 M1 s1 and kd = 8 × 102 s1 in the absence of noise (Fig. 9) and the presence of noise (Fig. 10)

The accurate determination of k_M suggests that BIACORE can be used to determine the diffusion coefficient of the analyte by calculating D from Eq. 9. The value of _M in Table 3A gives D = 9.85 × 10⁷ cm²/s, which is remarkably close to the value used in the simulation of 1.0 × 10⁶ cm²/s. However, because Eq. 9 is an approximate result, to determine D from k_M, one must ensure that the experiments are done in a parameter range in which Eq. 9 is sufficiently accurate (see Discussion).

Finally, to study the effects of noise on the best fit values of the parameters, we added 0.5 RU of pseudo-normally distributed noise to the data sets of Fig. 9. This level of noise is consistent with what is seen in typical experiments done on BIACORE 2000 (Morton and Myszka, 1998). Fig. 10 shows the results of fitting the two-compartment model to these data. The residual standard deviation was 0.504 RU, which is consistent with the level of noise added to the data. The best fit parameter values are given in Table 3B. Surprisingly, this level of random noise had little effect on the determined parameter values or the correlation coefficients, although as expected, the standard errors were increased. To demonstrate that the linear approximation statistics returned for the nonlinear analysis were valid, we performed Monte Carlo analysis to generate confidence probability distributions for the estimated parameters (Straume and Johnson, 1992). One thousand data sets were simulated with noise and solved as before with varied starting parameters. The results of this analysis confirmed that the confidence intervals were normally distributed (data not shown), and the standard deviations for the confidence intervals (Table 3C) were similar to those determined from the linear approximation statistics (Table 3B). The Monte Carlo analysis routines are readily available in the CLAMP program.

View larger version (23K): [in this window] [in a new window]   FIGURE 10   Simulations of the binding response of Fig. 9 in the presence of noise. Random noise with a standard deviation of 0.5 RU was generated and added to the simulations of Fig. 9. As in Fig. 9, the surface receptor concentrations correspond to a maximum response of 50 RU (a) and 15 RU (b) for a 3 × 104 Da analyte. Shown are simultaneous fits (red lines) of the two-compartment model to the simulations. The results of the fits are given in Table 3.

	DISCUSSION

Top Abstract Introduction Results Discussion Appendix References

The binding kinetics in a BIACORE flow cell are the result of diffusion and flow, which bring the analyte to, and take it away from, the sensor surface, and the chemical reaction that occurs between analyte and receptor at the sensor surface. We have modeled these processes for the case of a monovalent analyte interacting with a monovalent receptor, when the conditions under which the receptor is immobilized do not affect the binding kinetics. If, for example, the receptor is coupled to the sensor surface through a dextran layer, then our model is restricted to conditions where the layer has no influence on the binding kinetics. As discussed in the Appendix, numerous authors have considered this case and argued that a two-compartment model gives a reasonable description of the binding kinetics when the data are influenced by transport (Glaser, 1993; Karlsson et al., 1994; Schuck and Minton, 1996; Christensen, 1997), although its use in analyzing such data has been questioned (Schuck, 1997b). To see if, in fact, a two-compartment model can be used to extract accurate parameter values from transport-influenced BIACORE binding data, we generated test sets of simulated data. The obvious advantage of using simulated data sets is that the values of the underlying parameters are known exactly. We extended the model of Glaser (1993) to include diffusion in the direction of flow, obtaining a partial differential equation with appropriate boundary conditions that accounts for the processes that occur in a flow cell: convection, diffusion, and reversible chemical reaction at the sensor surface. The simulated data sets were generated by solving this model numerically, to obtain noiseless data, and then adding noise to the numerical solutions. These simulated data sets served as our standard for testing the two-compartment model. We showed that when the two-compartment model was fit to the simulated data, values for the rate constants were returned with surprisingly high accuracy.

The data we simulated corresponded to optimal experimental conditions currently available on BIACORE 2000, a flow rate of 100 µl/min and surface receptor densities as low as 15 RU. Our data analysis method involved simultaneously fitting the two-compartment model to two sets of simulated data, corresponding to two different receptor densities, 15 RU and 50 RU. The idea behind the approach, known as global analysis (reviewed in Myszka et al., 1997), is to obtain data that are sensitive to both the chemical rate constants (data from the low-capacity surface) and the transport coefficient (data from the high-capacity surface). In our fits we showed that even though the maximum response from one of the surfaces was only 15 RU, the addition of the expected level of experimental noise (0.5 RU) did not affect our ability to accurately recover the reaction rate constants.

It has been pointed out that if receptors are immobilized in a dextran layer, the binding kinetics can be influenced by the transport of the analyte within the layer (Schuck, 1996; Yarmush et al., 1996). It is expected that such effects will become pronounced at high receptor densities. At such densities the two-compartment model will no longer give a reasonable description of the binding kinetics. Elsewhere we discuss how to reformulate the model when the thickness of the layer cannot be neglected (Wofsy et al., manuscript in preparation). There we show that the effects of the layer will be negligible when the layer "looks" thin to the analyte, i.e., when the average distance the analyte travels in the layer before it binds to a receptor (its mean free path) is long compared to the height of the layer, d. The mean free path equals , where R_T/d is the three-dimensional receptor concentration and D_i is the diffusion coefficient in the layer. Thus one can ignore the layer and treat the binding as if it is directly to the sensor surface when 1  . (For the low receptor densities, we expect D_i  D. In our simulations the 50 RU surface corresponded to R_T = 0.167 nM cm and D = 1 × 10⁶ cm²/s. For a standard flow cell, d  1 × 10⁵ cm. For these values the inequality breaks down when k_a  1 × 10⁸ M¹ s¹.)

In addition to determining the rate constants, we showed that the two-compartment model can be used to determine D, the diffusion coefficient of the analyte. This was done by determining the transport coefficient, k_M, by fitting the simulated data, and then calculating D from Eq. 9. Because Eq. 9 is an approximate result, it can be used to determine D only under conditions in which it is highly accurate. Equation 9 overestimates k_M. It is derived by ignoring the parabolic flow velocity profile and taking the velocity to rise linearly with the height above the sensor surface. It also ignores the upper boundary of the flow cell. Analyte that diffuses to this boundary is not reflected. However, if the time it takes to travel the length of the flow cell, l/ = 3l/(2v_c), is short compared to the time it takes to diffuse to the boundary, h²/(4D), or equivalently, if 1  6lD/(h²v_c), then we expect Eq. 6 to be valid. Thus to determine D from Eq. 9, one should use a fast flow velocity for the binding studies.

We have shown that the two-compartment model can be used to determine accurate values of the rate constants and transport coefficients for data that are influenced by flow and diffusion. From the transport coefficient the analyte diffusion coefficient can be determined as well. Using the two-compartment model and analyzing the data globally extends the range of reaction rates that can be determined with instruments like BIACORE. At present, this approach seems to be the best available for analyzing experiments when the effects of mass transport cannot be eliminated from the binding kinetics. This model has been incorporated in BIAevalution 3.0 analysis software available from BIACORE, Inc. It has already been used to analyze such experimental data, yielding excellent fits to the data and improved accuracy of the returned rate constants (Myszka et al., 1996; Myszka et al., 1997; Morton and Myszka, 1998).