The study of how a ligand material, such as a hormone or antibody, binds to one or more kinds of molecular complexes, called sites, is of fundamental importance in biochemistry. Sites are often embedded in cell membranes, and the binding serves to control the behavior of the cell itself. Typically we are interested in the number of distinct kinds of sites and their frequency of occurrence, and also the equilibrium constants for each ligand-site binding reaction which indicates the absolute strength of each such binding reaction.

For example, quantitative analysis of hormone-receptor binding can be easily performed using appropriate software such as MLAB. MLAB is a computer program whose name is an acronym for "modeling laboratory"; it is an interactive system for mathematical modeling, originally developed at the National Institutes of Health. MLAB can fit multiple non-linear models to data points obtained from standard direct-binding or competitive displacement assays. Typical assays involve measuring the competition between radiolabelled and cold ligand in detergent-solubilized membrane preparations or on whole cells. Affinity constants and limit values of binding protein concentrations for single or multiple sites can be computed by fitting saturation curves in MLAB. Output can include Scatchard plots obtained by a suitable transformation.

There are two common categories of data. The first category results from cold displacement of bound labeled ligand where the labelled ligand concentration is held constant, the unlabelled ligand concentration is varied, and the ligand binding is calculated as a function of the proportion of labelled to total ligand concentration. The second category is a straight saturation analysis where the ligand binding is measured as a function of increasing labelled ligand concentration. This latter category of data also arises in fluorescent, ultraviolet, or electrical current studies. In general, MLAB can handle any measurements modeled by linear or nonlinear saturation plots.

In practice, the basic data must be adjusted to account for an experimentally-determined nonspecific binding level, specific activity of the labelled ligand, the values for maximal binding of the cellular preparations, and defective non-binding ligand. Replicate data points need not be averaged; they are entered as individual data points which are each included in the curve fitting process. In contrast, curve fitting based on dose averages may require careful attention to weights to obtain a reasonable fit. Experimental researchers rarely have a priori knowledge of the validity of individual data points, so each point is generally assigned an unbiased equal weight. However, weighting functions can be differentially applied in MLAB to adjust for differing error percentages, such as may occur in measurements of low levels of radioactivity. Because of the options for flexible weights, MLAB addresses some of the problems inherent in the treatment of experimental data where the number of samples is usually low and the error can be relatively high.

MLAB can be used to analyze down and up regulation in receptor number on the cell surface, where the binding capacity or receptor concentration is the limiting x-intercept on the Scatchard plot. The MLAB program is also useful in determining whether changes in binding by different ligand agonists/antagonists can be attributed to reduced binding affinity of the defined site/sites (seen as changes in the affinity constant), or to the total loss of binding ability of a class of sites (seen as changes in the binding capacity, and changes in the fit from an n+1 to n site model). Similarly, MLAB can aid in comparison studies of a single ligand against different binding proteins that vary by site-directed mutagenesis, or alternative splicing; such analyses can be useful in exploring mechanisms of ligand/receptor interaction, including amino acid and charge requirements.

In addition, post-translational processing and folding of the nascent protein in the endoplasmic reticulum can be studied by estimating the number of active receptor molecules that have reached the cell surface. These studies are of current interest to biochemists, and an accurate method of quantification is important since conclusions about the biological system, which in turn determine the future direction of research, are based on the calculated binding parameters.

In addition to saturation and Scatchard analyses, MLAB can be used in many other curve fitting applications such as kinetic, Michaelis-Menten, Lineweaver-Burke plots, and various statistical analyses. Kinetic analysis involves fitting differential equation models, which is an important core capability of MLAB.

The MLAB Mathematical Modeling System

MLAB has hundreds of useful functions, e.g., the discrete Fourier transform function dft and the parametric spline interpolation function splinep. One of the central components of the system is a curve fitting program which will adjust the parameters of a model function to minimize the weighted sum of the errors raised to a specified power. A repertoire of mathematical operators and functions, routines for solving differential equations, a collection of routines for onscreen drawing and for hardcopy plotting, and mechanisms for saving data between sessions provide a powerful and convenient environment for data manipulation, arithmetic calculations, and for building and testing models.

The user communicates with MLAB by typing commands which are executed at once or by providing a script to be executed. Should the user have questions, typing HELP will put the on-line system documentation at his disposal. The MLAB language is defined in the MLAB reference manual.

One of MLAB's main uses is to fit models to data. Curve-fitting is a useful analytical tool in many diverse disciplines. The basic notion is easily described. Given data, say various points in the plane (x1,y1), (x2,y2), ..., (xn,yn), and a function y = f(x) where f involves some parameters, say a and b, as for example f(x) = a*xb+1, we may wish to calculate values for the parameters a and b so that the function f well-predicts the observed data, that is, so that f(xi) = yi for 1 <= i <= n. In this case, we say we have fit the model f to the data by estimating the parameters a and b. The end result is merely the values obtained for the initially unknown parameters. The same principles apply in higher dimensions with arbitrarily many parameters. MLAB can simultaneously fit multiple non-linear model functions, some or all of which may be implicit functions, or may even be defined by a system of differential equations.

The curve-fitting and graphics display facilities of MLAB make it an ideal tool for the estimation of equilibrium constants from data, which typically consist of observed amounts of ligand bound for various specified amounts of ligand provided for binding.

The Mathematics of Multiple Site Binding

Suppose we have a ligand, F, which binds to each of N independently-acting sites, S1, S2, ..., SN, which are present in the various concentrations of S10mmol, S20mmol, ..., and SN0mmol respectively. Each reaction F + Si <--> Bi forms a bound complex Bi, and we shall define ki to be the associated equilibrium constant. Let F0 denote the concentration of ligand present initially, and let F be the concentration of free ligand at equilibrium. Similarly, let Si denote the concentration of free sites of type i at equilibrium, and let Bi denote the concentration of bound site i complex at equilibrium. Note we use the symbols F, Si, and Bi for the numerical quantities in millimole units of these materials, as well as for the names of the materials themselves. Then: ki= Bi/(FSi) and Bi + Si = Si0 for 1 <= i <= N, and F= F0 - B1 - B2 -...- BN.

Note that ki is specified in liters/mmole units, so that 106ki is the corresponding equilibrium constant in liters/mole units.

Often there is a fictitious (N+1)-st site, X, which binds F molecules, which is introduced to describe the non-specific binding of F molecules with weak affinity to many locations, other than the N specific sites of interest. Let k0 be the equilibrium constant for the fictitious non-specific binding reaction F + X <--> Y. Then we have:

ki = Bi/(FSi) and Bi + Si = Si0 for 1 <= i <= N,
k0 = Y/(FX) and Y + X = X0, and F = F0 - B1 - B2 -...- BN - Y,

where Y is the concentration of non-specifically bound ligand, X0 is the concentration of the fictitious non-specific binding site, and X is the concentration of the free fictitious site at equilibrium. Then Bi = kiSi0F/(1+kiF) for 1 <= i <= N, and Y = k0X0F/(1+k0F).

Now, to capture the notion of non-specific binding as a weak sticking of F molecules almost everywhere, let k0 tend to zero and let X0 tend to infinity such that k0X0 = c, where c is a fixed constant. Then Y -> cF, and we have: Bi = kiSi0F/(1+kiF) for 1 <= i <= N, Y = cF, and F = F0 - B1-...-BN- Y.

Note that the mathematical form Bi = Si0F/(1/ki+F) is difficult to handle when ki approaches 0, or when 1/ki+F approaches 0, so that the form given above is more convenient.

Usually values of F and/or B1+B2+...+BN+Y are measured for different values of F0 and we wish to use this data, which is generally of the form (F0,F) or (F0,B1+...+BN+Y), to estimate k1, k2, ..., kN, c, and, if not already known, S10, S20, ..., SN0. When only the total bound concentration B1+B2+...+BN+Y is measured, rather than B1, B2, ..., BN and Y separately, only a few kinds of sites can be distinguished by curve-fitting. In order to determine how many kinds of sites appear to be present, we must try each of the 1-site, 2-site, etc. models and choose among them on the basis of how well the data is fit. If the equilibrium constants obtained for two species of sites are close, then there are no grounds for considering them to be distinct species, based on this analysis alone.

Often people use the model B(F)=k1 S10F/(1+k1 F)+...+kN SN0 F/(1+kN F)+cF and fit it to data points of the form (F,B) where F is the free ligand concentration at equilibrium and B is the total bound ligand concentration at equilibrium, with one such point for each experiment. The difficulty with this approach is that we must compute F as F0-B so that there is correlated error in both the dependent and independent values used to form the data points. Our approach uses data points of the form (F0,B), and since F0 can usually be accurately determined, our model avoids the aforementioned difficulty of error in the independent variable.

It is possible to model the simultaneous use of several distinct types of ligands having distinct binding-affinities interacting with multiple classes of sites. This can be done by a straightforward elaboration of the equations used above. Although we shall not provide an example here, this is an important situation to be kept in mind.

The general two-site model with non-specific binding for a single ligand type can be defined in MLAB with the following dialog. Here and hereafter, the text shown following the MLAB prompt asterisk is an MLAB command statement entered by the user.

* FUNCTION B(F)  = k1*S10*F/(1+k1*F) + k2*S20*F/(1+k2*F)
* FUNCTION F(F0) = ROOT(Z,0,F0,F0-B(Z)-Z*(1+c))
* FUNCTION Y(F0) = c*F(F0)
Here A(F0) = B1+...+BN for N = 2.

These commands exemplify the MLAB FUNCTION statement, which is used to define a function or differential equation. Note that arguments of functions must be explicitly specified. Variables, such as k1 and k2, which appear in the body of a function, but not in its argument list, are called parameters. Parameters must be assigned values before an associated function can be evaluated.

ROOT is an operator which is built-in in MLAB. ROOT(Z,A,B,E) is a value between A and B which, when taken as the value of the dummy variable, Z, makes the expression, E, which involves Z, equal to zero. Thus ROOT(Z,A,B,E) is a solution of E(Z) = 0. The model given above, involving a so-called implicit function, deserves careful study; it is easily extendible to more than two sites. The amount F of free-ligand at equilibrium as a function of F0 satisfies F0-F=B(F)+c* F.

An Example

Suppose we have measured F in mmol units as a function of F0 in millimole units, as follows:

       F0      F
     .58668   .036
    1.1734    .096
    2.3467    .385
    2.9334    .61
    4.1068   1.15
    4.6934   1.46
    5.8668   2.11
    7.0402   2.73
and suppose we have S10 = S20 = 1.7121mmol. (If S10 and S20 were unknown, we could include them as fitting parameters.)

Then, we can introduce the appropriate constraints (which should always be used for this particular model), guess k1, k2, and c and then estimate them, as follows.

* constraints q ={k1>=0,k2>=0,c>=0,S10>=0,S20>=0}

The CONSTRAINTS statement permits the user to specify successive linear inequalities or equations involving the parameters (or potential parameters). Now we may specify values for the parameters, guessing when necessary.

* k1 = 10; k2 = 1; c = 0; S10 = 1.7121; S20 = S10;

Here the ASSIGNMENT statement is exemplified. In this case all the above assignments are assigning values to scalar variables, but the assignment statement is used to assign values to matrices as well. This can be seen in the next assignment statement which defines a matrix, D.

* D = Kread(8,2)
.58668 .036
1.1734 .096
2.3467 .385
2.9334 .61
4.1068 1.15
4.6934 1.46
5.8668 2.11
7.0402 2.73

The KREAD operator takes optional array size arguments; in this case an 8 row by 2 column matrix is specified, and reads in numbers from the keyboard to construct an appropriately-dimensioned matrix as the result. An analogous function is used to read numbers from a file. This matrix is then, in this case, assigned to D. Note the entry of the numbers which follow. Now we may examine D by "typing it out" using the TYPE statement.

* type D

D:   8 by 2 matrix
1:	.58668	.36E1
2:	1.1734	.96E-1
3:	2.3467	.385
4:	2.9334	.61
5:	4.1068	1.15
6:	4.6934	1.46
7:	5.8668	2.11
8:	7.0402	2.73

We have established a model function, F, and entered data, D. Expect for error, we expect that f(D[i,1]) = D[i,2] would hold, if only the parameters k1, k2, and c were set to their "correct" values. The following FIT statement requests MLAB to estimate k1, k2, and c by assigning them values which minimize the sum-of-squares objective function S(k1,k2,c) =SUMi=1,...,8(F(D[i,1])-D[i,2])2.

* maxiter = 30; TOLSOS = .001
* fit(k1,k2,c), F to d, constraints q
final parameter values
value           error           dependency      parameter
13.86352736     2.331373145     0.665073064     K1
0.5321372226    0.06602915638   0.9432600296    K2
0.5874015019    0.02229743988   0.9171690302    C
5 iterations
best sum of squares = 9.78763e-04
root mean square error = 1.39912e-02
deviation fraction = 6.82654e-03
R squared = 9.99854e-01
no active constraints

The behavior of the fit statement depends upon the supplied constraints q, as well as upon the MLAB control variables: maxiter, the maximum number of iterations and tolsos, the requested convergence factor.

MLAB uses a carefully-tuned version of the Marquardt-Levenberg magnified-diagonal algorithm which is, in turn, a form of the Gauss-Newton procedure for minimizing a function which is in the form of a sum-of-squares. This process estimates the value of the parameter vector b=(k1,k2,c)' by successive approximations b(0), b(1), ..., b(n), where b(0) is the vector of initial guesses for k1, k2, and c, and b(j+1) = b(j) + beta(j), where

beta(j)=(X'V-1X+epsilon G)-1 X'V-1(y -(f(x1;b(j)), ..., f(x8;b(j)))'), with
Xst= partial f(xs; b(j))/partial bt and
Gst= if s=t then (X'V-1X)st else 0 and
xs=D[s,1] for 1 <= s <= 8, and
y=(D[1,2], ..., D[8,2])',

where V is the estimated covariance matrix of the observations. In our example, V=I, the identity matrix. In general V is determined from weight-values supplied by the experimenter.

An iteration consists of computing b(j+1) from b(j). Note that this requires the partial derivatives of the model function with respect to the parameters evaluated at b(j), since these values form the matrix X. In MLAB, these derivatives are automatically computed symbolically and evaluated to form X. The convenience thus obtained is considerable and the parameter estimation process is provided with more accurate derivative values. For example the derivative of F with respect to k1 can be explicitly displayed in MLAB as follows.

* type F'k1
FUNCTION F'K1(F0) = EVAL(Z,ROOT(Z,0,F0,(F0-B(Z))-Z*(1+c)),

Indeed derivatives are full-fledged members of the class of functions and can be used in graphics or curve-fitting in MLAB just as can any other user-defined function.

A sub-iteration consists of computing b(j+1) with a particular value of epsilon which specifies the amount of diagonal magnification. At each iteration, the value epsilon starts at 10-9 and is increased until the corresponding value of b(j+1) results in a smaller sum-of-squares value, whereupon this vector is taken to be the final b(j+1) iterate.

The parameter estimation process stops when the limit of the number of iterations is reached or, more usually, when the decrease in the sum-of-squares value between successive iterations is less than a specified fractional amount determined by the user-specified convergence factor in TOLSOS. For TOLSOS = .001, the sum-of-squares must change by less than .1 percent for the curve-fitting process to stop based on this criterion.

When the estimation process does stop, the parameters are reset to their computed estimates, and they and their estimated standard deviations are typed-out. Associated values called dependency values, which lie between 0 and 1, are also typed-out. It suffices here to remark that large dependency values above .99 usually indicate a non-unique solution; that is, other parameter estimates exist which would provide a nearly equally-small sum-of-squares.

MLAB also types-out the root-mean-square error which is the estimate of the standard deviation in each observation, given that they are identically distributed. This quantity should roughly equal the experimental error in the data. There are, of course, many caveats and restrictions which must hold to insure the validity of the supporting statistics provided.

The material typed out above shows that the vector of parameters (k1,k2,c) has been estimated to be (13.8633, .532151, .587396), with reasonably small dependency values and standard error values, and with an RMS error of about .014, which should be comparable with the experimental error in our data, D. The sum-of-squares was reduced from an initial value of 2.64314 to .000978763 at the final parameter values.

In order to visually see how our model with its parameters set to their best-estimated values corresponds to the data, we may draw a graph of the data points and the model function. Although we are drawing only the simplest and most direct kind of picture here, it should be noted that MLAB provides facilities for many types of point-symbols and types of lines, axes with arbitrarily-placed numeric labels in various formats, titles in the form of text strings in arbitrary sizes and various fonts with subscripts and superscripts, color, and a number of other special features. It is quite possible to prepare more or less elaborate publication-quality graphs with a modest amount of effort. Indeed this is one of MLAB's most-used facilities. The desired graph can be constructed as follows.

* draw D, linetype none,pointtype circle
* draw c2 = points(F,0:8:.2)

The first DRAW command above plots the data points, while the second draw command constructs a curve called C2, which is a graph consisting of solid straight-lines connecting the points which are the rows of the 2-column matrix which is the value of the expression POINTS(F,0:8:.2). This matrix has the values 0 through 8 in steps of .2 in its first column and corresponding values of the function F evaluated at 0 through 8 in steps of .2 in the second column. The POINTS operator is very useful for graphing functions. Both these curves are drawn in the default MLAB 2D-window called W (since no other window is specified) which has predefined labeled axes already present. The picture finally appears when its display is requested with a view statement and a plot can be obtained if desired using the PLOT statement.

Often a Scatchard graph of Bound/Free vs. Bound (i.e., B(F(F0))/F(F0) vs. B(F(F0))) is desired. Although such a formulation should not be used for curve-fitting due to the non-normal error introduced by computing Bound/Free, it is quite straightforward to draw the data and model Scatchard plots as follows. Note the current picture in W is saved and restored.

* FUNCTION R(X,Y)=IF Y=0 THEN (k1*s10+k2*s20) ELSE X/Y
* MF = F ON 0:8:.2
* M = B ON MF
* M = M&'(R ON M&'MF)
* MF = D COL 2
* M = (D COL 1)-(1+C)*MF
* M = M&'(R ON M&'MF)
* TOP TITLE "Bound/Free vs. Bound"

The &' operator denotes column concatenation, while the ON operator, as in F ON H, applies the function F to each row of the matrix H treated as an argument list for F and returns the column vector of result values.

It is an enlightening exercise to fit the Scatchard model to the corresponding transformed data and to then draw a graph of the original data points and the function f using the parameter values obtained. An example of this comparison is given in the following section.

If we know we have non-cooperative two-site specific binding, with our given S10 and S20 values, together with non-specific binding with no degradation (degradation could be due to the presence of various enzymes for example), then this curve-fit of the correct model gives us the best available estimates of k1, k2, and c. But if the model is possibly not correct, we have a more difficult problem of deciding what the correct model is. This often cannot be resolved without further experimentation.

The non-specific binding constant, c, can be independently estimated by adding a large amount of unlabeled ligand, L, which is distinguishable from F, to a mixture of F and the site-bearing material. Then virtually all the bound F which results will be non-specifically bound and taking this value as Y yields the approximation c = Y/(F0-Y). With a value of c known from one or more such experiments, c can be fixed in subsequent curve-fitting, or alternatively the data can be modified by using F-cF for F and fitting an appropriate model with c = 0.

If we have no such independent estimate of c, and we are unsure if non-specific binding is occurring, we may test the unrealistic possibility that there is no non-specific binding as follows.

* tolsos = .0002
* c = 0
* fit(k1,k2), F to d, constraints q
final parameter values
value           error           dependency      parameter
10.42883512     881.9654883     0.9989885554    K1
10.87461407     712.3997932     0.9989885554    K2
27 iterations
best sum of squares = 1.31885e+00
root mean square error = 4.68838e-01
deviation fraction = 2.57765e-01
R squared = 8.02989e-01
no active constraints
LSQRPT is a sum of zero or more of the values 1, 2, 4, and 8. When LSQRPT contains 1, the best sum of squares obtained at each iteration is typed-out. When LSQRPT=1, the fit command directs MLAB to type-out only the final results. All intermediate reports are suppressed. It is also possible to use LSQRPT=0 which suppresses all output, but the parameters are reset as usual. In any event, a matrix called COVP (the estimated covariance matrix of the parameters) and a scalar called SOSQ (the final sum-of-squares value) are created, replacing any previous data objects with these names. The user may examine these values, or employ them in further calculations as desired.

We may examine our fit by redefining the curve C2.

* draw c2 = points(F,0:8:.2); view

This is not an adequate fit, as we could have predicted by noting the large parameter standard deviations and the large RMS error. Although the dependency values are not extremely close to one, we may consider the possibility that our fit is bad because we fell into an unfortunate local minimum. A crude way to study this possibility is to draw a contour map of the sum-of-squares surface to see if several local minima are apparent. This is possible in this case, since we have fewer than three parameters, and MLAB can be used to view the sum-of-squares surface as follows. Drawing this contour map is time-consuming due to the large number of root calculations required.

* FUNCTION FH(F0,K1,K2)=ROOT(Z,0,F0,F0-(K1*S10*Z/(1+K1*Z) 
* FUNCTION S(k1,k2)=SUM(I,1,8,(FH(D[I,1],k1,k2)-D[I,2])2)
* M=contour(points(S,cross(0:20,0:20)),1.3:1.35:.005&1.4:1.8:.2&2:5)

The CROSS operator expression above constructs a 441 by 2 matrix whose rows are all the pairs (i,j) with 0 <= i,j <= 20. The CONTOUR operator expression constructs a specially-coded two-column matrix which is suitable for drawing with line-type VMARKER. Contour lines are formed for the surface height values in the list 1.3:1.35:.005 & 1.4:1.8:.2 & 2:5 where & denotes row concatenation. Evidently, our solution is unique.

We may try a one-site model with S10 equal to 3.4242 mmol by using our two-site model where both sites are identical as follows. If S10 were not known, we could add S10 to the parameter list to be adjusted.

* S10 = S10+S20; k2 = 0;
* fit(k1,c), F to d, CONSTRAINTS Q
final parameter values
value           error           dependency      parameter
2.325176617     0.2567377717    0.5948957102    K1
0.4642651058    0.03442515131   0.5948957102    C
15 iterations
best sum of squares = 1.55348e-02
root mean square error = 5.08835e-02
deviation fraction = 2.31134e-02
R squared = 9.97679e-01
no active constraints

*delete w
*draw D, line none, pointtype circle
*draw points(f,0:8:.2); VIEW

This fit is not as good as the fit obtained with the two-site model, and we can probably reject it. We should have an estimate of the variance of the data in order to test the goodness-of-fit. Even with curve-fitting evidence, it is best if there is independent evidence of multiple sites. Also, of course, a three or more site model, with or without non-specific binding may be the correct model.

Thus, in general, given the correct model, curve-fitting is a powerful device for estimating parameters; but curve-fitting is not always a very good way to distinguish the correct model from alternate impostor models. Usually we must resort to discriminating experiments.

MLAB can handle simultaneous curve-fitting, where several functions are to be fit to corresponding sets of data points and where, moreover, these several functions share some parameters in common. MLAB model functions can be functions of more than one formal argument, so that data points may lie in n-space for arbitrary n, but we do not provide an example here. However we shall provide an example of simultaneous curve-fitting.

Suppose we not only have the data given above for 2-site binding, but also are given the following observations of the concentrations Y of F-ligand which is non-specifically bound, for various amounts of ligand provided initially, denoted by F0.

   F0	   Y
  1.1734   .3
  2.5      .5
  3        .3
  3.98     .7
  6.5     1.6

Then we may enter this data as the rows of a matrix N and estimate k1, k2, and c based on all the available data.

* N=KREAD(5,2)
1.1734	.3
2.5	.5
3	.3
3.98	.7
6.5	1.6

In general, we should use weights which specify the relative accuracy of each data point. For each data matrix, a vector of weight values may be given as a clause in the FIT statement. The weight-value, gi, that corresponds to an observation (xi,yi) should be 1/var(yi). Of course, guesses or estimates must generally suffice. The value 1/yi2 for yi corresponds approximately to a constant percentage error in the observed value yi. Weight functions are permitted in MLAB, and a more accurate device is to use the reciprocal of the square of the model function itself as the weighting function; this results in an iterative reweighting process.

The use of weights permits data points of high reliability to exert greater influence on the parameter estimates than data points of lesser reliability. When simultaneous fitting is being done, weights have the additional function of compensating for differing units or magnitudes of the observations from various data sets. Without weights, the curve-fitting process would tend to favor one model component, fitting it at the expense of others, if the deviations there were much larger than the others, even though this may be an artifact of the use of different units.

MLAB includes an operator, EWT, which computes a vector of weight-values for a given set of data points, M, by estimating the standard deviation at each point by the difference between the data curve and a smoothed form of it, which effectively assumes the error is due to white noise. These standard-deviation estimates are, in turn, smoothed, and then used to form reciprocal variance weight values. In our example, we shall use EWT to obtain the necessary weight vectors, although this is not totally appropriate for chemical binding data, which is often taken in a range where the measurement accuracy improves as the amount of bound ligand increases. Also, EWT is unreliable for small numbers of data points. Nevertheless, we shall proceed as follows.

* DW = EWT(D); DN = EWT(N)

: a  8 by 3 matrix
1:  0.58668    0.036    340.092919
2:  1.1734     0.096    340.092919
3:  2.3467     0.385    271.140209
4:  2.9334     0.61     256.511957
5:  4.1068     1.15     193.017413
6:  4.6934     1.46     171.441311
7:  5.8668     2.11     93.14747  
8:  7.0402     2.73     93.14747  

 a  5 by 3 matrix
1:  1.1734     0.3      159.445918
2:  2.5        0.5      159.445918
3:  3          0.3      156.341715
4:  3.98       0.7      134.665471
5:  6.5        1.6      134.665471

* fit (k1,k2,c), F to d with weight dw, Y to n with weight dn, CONSTRAINTS q
Begin iteration 1 bestsosq=5.86711e+01
Begin iteration 2 bestsosq=2.76531e+01
Begin iteration 3 bestsosq=2.69881e+01
final parameter values
value           error           dependency      parameter
1.681292363     0.3984166729    0.6190759098    K1
1.905413519e-17 0.06987097833   0.6367292361    K2
0.6401358429    0.057666687     0.09280130233   C
3 iterations
best sum of squares = 2.69879e+01
root mean square error = 1.64280e+00
deviation fraction = 6.88260e-02
lagrange multiplier[1] = -0
lagrange multiplier[2] = -73.23991748
lagrange multiplier[3] = -0
* delete w
* draw d line none pt circle
* draw points(F,0:8:.2)
* draw n line NONE pointtype "+"
* draw points(Y,0:8:.2) linetype dashed

As mentioned above MLAB also allows weights to be the reciprocals or other functions of the actual squared deviations themselves. This device is called curve-fitting with iterative reweighting and is invoked in MLAB by specifying an appropriate weight function. Moreover MLAB employs an iterative reweighting technique to allow the general form of a sum of pth powers to be minimized, rather than just a sum-of-squares. This so-called Lp-norm fitting is sometimes useful for p = 1 or 1.5 when the data is not normally-distributed, and tends to contain some outliers. The value p is specified as the value of the MLAB control variable fitnorm.

Practical Direct and Competitive Ligand Binding Analysis

Given a radioactively-labeled ("hot") ligand material F, we wish to mix varying amounts of hot ligand each with the same fixed amount of site material, and obtain a table of the amounts of bound ligand corresponding to the total amounts of ligand added. This requires that we be able to separate the bound ligand from the free ligand and then measure the amount of bound ligand by the numbers of radioactive-decay events per minute for each distinct experiment.

There is an alternate device which sometimes has an advantage in simplicity over this direct method. Let us add a single fixed amount of hot ligand to the site material and then add varying amounts of unlabeled ("cold") ligand. The cold ligand form of F competes with the hot ligand in binding to the sites, and the effect is that as we add more cold ligand, the observed concentration of bound hot ligand decreases. We may sometimes be able to add an incremental amount of cold ligand to our mixture, separate and measure the amount of bound hot ligand, and then remix and add the next incremental amount of cold ligand; when this is possible, the variability due to differing concentrations of sites and hot ligand is reduced. Since replicate experiments suffer from using slightly-differing amounts of site material and of hot ligand, it is often better to expend effort getting more points for a single sample of site material than to run replicate experiments (except as a gross check on the experimental process.) In any event, the results of replicate experiments need not be averaged; the MLAB curve-fitting process will, in effect, give each point equal weight as is appropriate.

We have a further complication to deal with. Both the hot and cold ligand materials usually bind weakly and non-specifically throughout our mixture, in addition to specifically binding with the site material. The amount YH of such non-specifically bound hot ligand is generally taken to be proportional to the total amount of hot ligand present less the amount bound, i.e. YH = cFH for some constant c where FH denotes the concentration of free hot ligand at equilibrium. Let AH denote the concentration of bound hot ligand, both specifically and non-specifically, and let BH denote the concentration of specifically-bound hot ligand. The concentration of hot ligand in all forms is LH=AH+FH. Then AH=BH+YH. If we add a very large amount of cold ligand to our mixture, almost all the specific binding sites will be occupied by cold ligand. The amount of observed bound hot ligand AH may then taken to be the amount of non-specifically bound hot ligand YH, whence c is approximately equal to YH/FH = AH/(LH-AH). We can thus use our other observed AH values and compute corresponding BH values by subtracting c(LH-AH) from the AH values, where c has been estimated as indicated above. (Alternatively, it is possible, and generally preferable, to curve-fit for the unknown non-specific binding constant c as is done below.)

In any event, for each observation point, we know the constant Molar concentration LH of hot ligand present, the concentration LC of cold ligand present and the observed concentration AH=BH+YH of bound hot ligand (with due regard for the volume of the mixture). Computing the concentration AH requires a computation based on the following parameters.

  1. the molecular weights of hot and cold ligand (often these are taken to be the same.)
  2. the half-life of the radioactive-label material and the age of the hot ligand being used, measured in days.
  3. the counting efficiency of the radioactive decay detector device, measured as the ratio of observed counts per minute to the actual number of decay-events per minute.
  4. the counting time period (minutes) during which decay events are monitored.
  5. the background radioactivity (in counts per time-period) in the site material and the cold ligand taken together.
  6. the specific activity of the hot ligand measured as the expected number of decays per minute per mole of hot ligand of age 0.
  7. the total counts over the counting time period for the total amount of hot ligand material.
  8. the bindable fraction of the hot ligand defined as the fraction of the radioactive label material attached to non-inert ligand which is actually capable of binding.
  9. the volume (milliliters) of the mixture.
  10. the mass amounts (nanograms) of the cold ligand and the hot ligand in our mixture.
  11. the observed counts during the counting time-period for the bound ligand component.
Let WH be the molecular weight of hot ligand.
Let WC be the molecular weight of cold ligand.
Let HL be the half-life in days of the label material.
Let AL be the age in days of the label material.
Let CE be the counting efficiency in (cpm/dpm) units.
Let CT be the counting time period in minutes.
Let BR be the background radioactivity in counts per time period.
Let SA be the specific activity of the hot ligand in dpm/mole units.
Let TC be the total counts per time-period of the hot ligand material.
Let BD be the bindable fraction of hot ligand material.
Let V be the volume in milliliters.
Let MH be the mass in nanograms of the hot ligand material.
Let MC be the mass in nanograms of the cold ligand material.
Let AO be the observed counts per time-period for the bound ligand component.

Then we have LH = MH/[WH * V* 106] Molar units (M) of hot ligand. As a check, LH = (1/SA)[TC/(CT * CE)] * 2AL/HL * (103/V). This equation can alternately be used to determine any of the other parameters involved, given LH in Molar units. In particular, if the observed value of TC is not close to the predicted value from this equation, there is some difficulty in the experimental parameters. The corrected value of LH we use in subsequent calculations is then obtained by multiplying LH by the bindable fraction factor BD. The concentration in Molar units of the bound hot ligand is computed as AH = (1/SA)[(AO-BR)/(CT * CE)] * 2AL/HL* (103/V).

The concentration of cold ligand is LC = MC/[WC * V* 106] M. (If a fraction of cold ligand is inert, the corresponding correction should be made to LC.) A series of competition experiments consists of varying MC and observing the differing values of AO (and potentially BR) that result. The final outcome is a sequence of (LC,LH,AH) triples, where LC and AH vary, and LH is generally fixed. For a series of direct-binding experiments, MC is 0, and we vary MH directly. The result is a sequence of (LC,LH,AH) triples where LC = 0 and LH and AH vary.

Let us consider a single-site material S. Suppose our mixture contains a concentration of S0 M of this site material together with LC M of cold ligand and LH M of hot ligand. Our goal is to estimate the equilibrium constant KC of the reaction LC + S <--> BC. Let BC denote the concentration of bound cold ligand at equilibrium and let FC denote the concentration of free (unbound) cold ligand at equilibrium, just as BH and FH denote these quantities for hot ligand. Suppose the reactions LH + S <--> BH and LC + S BC have the not-necessarily-equal equilibrium constants KH and KC, so that KH = BH/(FH SF) and KC = BC/(FC SF), where SF denotes the concentration of free sites at equilibrium. Also, BH + FH+YH = LH, BC + FC+YC = LC, and SF + BH + BC = S0. All concentrations are in Molar units (M).

Our goal is to estimate KC, given the varying observed concentrations of bound hot ligand associated with the sequence of cold ligand concentrations used. This cannot be done without knowing KH, or alternately, the ratio p := KH/KC. Thus KC = KH/p, and either KH or p, but not both, is to be determined. Often we may assume p = 1.

A pre-computer-style Scatchard-plot analysis method is often employed to estimate the equilibrium constant KC. When LC = 0, suppose that KC = KH/p, where p is known. Also we have KH = BH/(FH SF) and BH + SF = S0. Then BH = S0 - 1/KH} * BH/FH, or BH/FH = -KH(BH - S0). This last relationship is a linear relation between BH/FH and BH. If we graph BH/FH vs. BH for varying amounts of LH, we obtain a straight-line graph with slope -KH and x-axis intercept S0 and y-axis intercept KH S0. Note that BH/FH is a dimensionless quantity, so we could compute BH and FH separately in any desired units (such as cpm), and obtain the same ratio when BH/FH is formed. This is the so-called Scatchard plot commonly used when a more advanced analysis is not available; the unknown parameters KH and S0 are estimated directly from looking at a straight-line fit of the experimental points (BH,BH/FH). Then KC = KH/p. As LH increases, both BH and FH increase, but FH increases faster; this is reflected in the fact that the Scatchard plot line has a negative slope. If the plotted points do not approximately lie on a negatively-sloped straight line, the one-site binding model that has been assumed is probably incorrect. (What do you think about the example presented below?)

For a sequence of competitive binding experiments, the crucial point is that when we have 100 alpha % hot ligand and 100(1-alpha ) % cold ligand, then both forms of ligand participate in every chemical state in this proportion (modified by the equilibrium proportionality constant p.)

When p=1, so that KC = KH, we then have BH/FH = -KH[BH+BC-S0]. And BH = alpha (BH+BC) and BC = (1- alpha)(BH+BC), where alpha := LH/(LH+LC). Thus BC = ((1- alpha)/ alpha) BH, so BH/FH = -KH[BH/alpha-S0], analogous to the non-competitive binding case!

In general, however, it is preferable to use the data (F0,B), where F0 = LC+LH and B = AH/alpha = (BH+YH)/ alpha = (BH+YH) F0/LH, and fit the model given above for multiple-site binding with a single ligand where the independent variable is the total concentration of (hot and cold) ligand supplied, and the dependent variable is the concentration of (hot and cold) ligand bound (specifically and non-specifically). When p is not 1, using this model is required.

Thus we may construct the data points (F0,B) where F0 = LC + LH and B = AH F0/LH for each experimental data point (LC,LH,AH), and then fit the function A to this data in order to estimate the parameters c, S10, K1, (and S20, K2, if a two-site model is desired). When p is not 1, we introduce the constraints K2 = K1/p and S10 = S20 and use the two-site model for the one-site situation with KC is not KH.

A general do-file based on the MLAB statements given above can be constructed which queries the user for the experimental parameters: HL, AL, CE, CT, BR, SA, TC, BD, V, MH and for the sequence of value-pairs (MC,BO), and then fits the desired model, reports the estimated values of c, S10, and K1(=KC), and graphs the corresponding saturation plot and Scatchard plot for the data. A simplified example of such a do-file (called is shown below, together with an example of its use.

" = competitive single site binding analysis"
echodo=0 /* echodo=0 inhibits screen and log-file output */

type "Specify the COLD ligand molecular weight"

type "Specify the HOT ligand molecular weight"

type "Specify the half-life in days of the HOT ligand label"

type "Specify the age in days of the HOT ligand label"

type "Specify the counting-efficiency (c+/-/d+/-) of the label"

type "Specify the counting time-period in minutes"

type "Specify the background radioactivity in cpp(counts per period)"

type "Specify the specific activity of the HOT ligand(in d+/-/mole)"

type "Specify the total counts per period of HOT ligand(in cpp)"

type "Specify the HOT ligand bindability fraction"

type "Specify the Volume (in milliliters)"

type "Specify the mass in ng of the HOT ligand
 (or enter -1 to have the mass computed for you.)"

/*Compute the mass in ng of the HOT ligand! */
type "specified mh value (in ng):", mh
type "computed mh (hot mass) value (in ng) based on tc:",mh

type "Specify the number of experimental observation points"
type "type-in the sequence of "+strval(nb)+" data points"
type  "as: ([COLD ligand mass in ng], [observed bound cpp]) pairs"
type  "on successive lines"
d=kread(nb,2); d=sort(d)

/* 2 classes-of-sites binding model.
   We use k2=0 and s20=0 for the 1-class-of-sites specialization*/
FCT A(F0)=G(F(F0))
FCT B0(F0)=B(F(F0));
FCT B(F)=K1*S10*F/(1+K1*F)+K2*S20*F/(1+K2*F);
FCT F(F0)=ROOT(Z,0,F0,F0-B(Z)-Z*(1+C));
FCT Y(F0)=C*F(F0);
CONSTRAINTS CQ={K1>=0,K2>=0,C>=0,S10>=0,S20>=0};

type "Molar conc. of HOT ligand based on supplied molecular wt.",lh
type "Molar conc. of HOT ligand based on total cpp supplied",lhck

lh=lh*bd /* compensate for inert hot ligand */

fct m0(x)=if x>0 then x else 0
fct cconc(mc)=mc/(wc*v*1e6)
fct bconc(a0)=((a0-br)/(ct*ce))*(1e3/sa)*(2^(al/hl))/v

dm col 1 = cconc on (d col 1)  /*total cold (M)*/
dm col 2= m0 on bconc on (d col 2)   /*bound + non-spec. bound (M)*/

/*type-out the raw and converted input data*/
type "D=input data [Cold Ligand (ng) | Hot counts]", d
type "DM=[Cold Ligand (M) | Hot Spec.+Non-Spec.bound(M)]",dm

/* compute da col 1= LH+LC, da col 2= AH/alpha */
da=lh+(dm col 1)            /*total (hot+cold) */
alpha=lh/'da                /* vector of scale-factors*/
da col 2=(dm col 2)/'alpha  /* scaled bound+non-spec. bound */

/* compute guesses for K1 and S10 from the Scatchard fit.
   estimate C from the largest cold ligand data point. */

i=maxrow(dm col 1);

fre=lh-(dm col 2)            /*free hot (FH) in Molar units */
yns=c*fre                    /*non-spec. bound hot */
sd =m0 on ((Dm COL 2)-yns);  /* spec. bound hot*/
sd2 = sd/'fre                /*bound/free */
sd col 2=sd2
sd col 1 = (sd col 1)/'alpha /* spec. bound (all =H+C) */

type "directly-estimated non-spec. binding coef.", c

ltot=(da col 1)
ftot=((da col 1)-(da col 2))
btot = sd col 1

type "MMT=[L Total | B Total | Y Total | F Total]", mmt
type "MMH=[L Hot | B Hot | Y Hot | F Hot]", mmh
type "MMC=[L Cold | B Cold | Y Cold | F Cold]", mmt

/*Remember BH/FH = -K1*BH/alpha +K1*S10 */

fct lin0(b)=-k1*b+k1s10 /* simple linear form */
fct lin(b) = -k1*(b-s10)
k1=1; k1s10=lh

fit(k1,k1s10), lin0 to sd
fit(k1,s10),lin to sd
type "Scatchard model linear-regession fit estimates",k1,s10
ok1=k1; os10=s10

draw sd lt none pt circle color orange
draw points(lin,0:s10!4) lt dashed
top title "Scatchard plot"
title "[dashed=regression-fit line]" at (.5,.75) fract size .14 inches
title "[hashed=euclidean-fit line]" at (.5,.7) fract size .14 inches

/* compute k1 and s10 for the best euclidean-fitting line */
z=eigen(m) row (1,3,5)
rno=if z[1,1]>z[1,2] then 2 else 3
v1=(z row rno)'
type "Scatchard model euclidean-fit estimates;", k1,s10
draw points(lin,0:s10!50) pt nbar ptsize .01
ws=w; blank ws;

k2=0; s20=0;
/*define our weighting function (fixed-percentage error assumed) */
fct wvf(x)=1/x^p>
wv1=.00000001*(wvf on (da col 2)); 

FIT(K1,s10,C),A TO da with wt (wv1), CONSTRAINTS CQ;

/* Draw the (full) relative saturation curve log-plot form ------------*/
FL = (da[1,1]*.9):((da[nrows(da),1])*1.1)!100;
am= POINTS(a,FL);
lda= (log on da col 1)&'((da col 2)/'(da col 1))
draw lda lt none pt circle ptsize .01
lam =(log on am col 1)&' ((am col 2)/'(am col 1))
draw lam

nk1=k1; ns10=s10; k1=ok1; s10=os10
am = points(a,fl)
k1=nk1; s10=ns10
lam =(log on am col 1)&' ((am col 2)/'(am col 1))
draw lam lt dashed
top title "Log-Percent Saturation Plot"
bottom title "log(TOTAL LIGAND)"
left title "(spec. bound)/(total)"
title "[dashed=regression-fit line]" at (.5,.8) fract size .15 inches
title "[solid=MLAB model-fit line]" at (.5,.75) fract size .15 inches

/*Now draw the Scatchard form of the data +fit-curves -------*/
w = ws; unblank w
draw points(lin,0:s10!4) color yellow
TOP TITLE "Scatchard Plot: Bound/Free vs. Bound" color green;
title "[solid=MLAB model-fit line]" at (.5,.8) fract size .14 inches
left title "(bound)/(free)"
bottom title "bound ligand"
lowy=minv(sd col 2); hiy=maxv(sd col 2)
window 0 to 2*s10, lowy to hiy adjust wnice

The MLAB log-file result of running the above do-file on some actual experimental data is shown below.

   MLAB Mathematical Modeling System, Revision: April 24, 1996
   Copyright: Civilized Software, Inc. (301)962-3711, email:

Sat May 18 15:29:38 1996
'* ' is the command prompt

* do liganal

Specify the COLD ligand molecular weight
wc=  1100
Specify the HOT ligand molecular weight
wh=  1130
Specify the half-life in days of the HOT ligand label
hl=  2
Specify the age in days of the HOT ligand label
al=  0
Specify the counting-efficiency (cpm/dpm) of the label
ce=  .5
Specify the counting time-period in minutes
ct=  5
Specify the background radioactivity in cpp(counts per period)
br=  0
Specify the specific activity of the HOT ligand(in dpm/mole)
sa=  2.5086e18
Specify the total counts per period of HOT ligand(in cpp)
tc=  65789
Specify the HOT ligand bindability fraction
bd=  1
Specify the Volume (in milliliters)
v=  .5
Specify the mass in ng of the HOT ligand
(or enter -1 to have the mass computed for you.)
mh=  -1
     specified mh value (in ng):
     computed mh (hot mass) value (in ng) based on tc:

Specify the number of experimental observation points
nb=  10
     type-in the sequence of 10 data points
     as: ([COLD ligand mass in ng], [observed bound cpp]) pairs
     on successive lines
  Type in numbers for KREAD.  End lines with ENTER
  0     8118
  .05   7694
  .1    7480
  .2    7080
  .5    6360
  1     6170
  3     5397
  10    5096
  30    5096
  1000  4714

     Molar conc. of HOT ligand based on supplied molecular wt.
     Molar conc. of HOT ligand based on total cpp supplied

     D=input data [Cold Ligand (ng) | Hot counts]

   D: a  10 by 2 matrix

  1: 0      8118
  2: 0.05   7694
  3: 0.1    7480
  4: 0.2    7080
  5: 0.5    6360
  6: 1      6170
  7: 3      5397
  8: 10     5096
  9: 30     5096
 10: 1000   4714

     DM=[Cold Ligand (M) | Hot Spec.+Non-Spec.bound(M)]

   DM: a  10 by 2 matrix

  1: 0                2.58885434E-12
  2: 9.09090909E-11   2.45363948E-12
  3: 1.81818182E-10   2.38539424E-12
  4: 3.63636364E-10   2.25783305E-12
  5: 9.09090909E-10   2.02822291E-12
  6: 1.81818182E-9    1.96763135E-12
  7: 5.45454545E-9    1.72111935E-12
  8: 1.81818182E-8    1.62512955E-12
  9: 5.45454545E-8    1.62512955E-12
 10: 1.81818182E-6    1.50330862E-12

     directly-estimated non-spec. binding coef.
  C = 7.71837904E-2

     MMT=[L Total | B Total | Y Total | F Total]

   MMT: a  10 by 4 matrix

  1: 2.09803077E-11   1.16933226E-12   1.41952208E-12   1.83914534E-11
  2: 1.11889399E-10   5.45935993E-12   7.62606519E-12   9.88039735E-11
  3: 2.0279849E-10    9.18445554E-12   1.38730875E-11   1.79740946E-10
  4: 3.84616671E-10   1.48997633E-11   2.64914424E-11   3.43225466E-10
  5: 9.30071217E-10   2.50658571E-11   6.48466349E-11   8.40158725E-10
  6: 1.83916213E-9    4.38447819E-11   1.28640441E-10   1.6666769E-9  
  7: 5.47552576E-9    6.12326499E-11   3.87952055E-10   5.02634106E-9 
  8: 1.82027985E-8    1.13851314E-10   1.29613305E-9    1.67928141E-8 
  9: 5.45664349E-8    3.41291496E-10   3.88541133E-9    5.0339732E-8  
 10: 1.8182028E-6     1.16933226E-12   1.30280259E-7    1.68792254E-6 

     MMH=[L Hot | B Hot | Y Hot | F Hot]

   MMH: a  10 by 4 matrix

  1: 2.09803077E-11   1.16933226E-12   1.41952208E-12   1.83914534E-11
  2: 5.96713722E-10   2.91151354E-11   4.06703209E-11   5.26928266E-10
  3: 1.9602776E-9     8.87781883E-11   1.34099139E-10   1.73740028E-9 
  4: 7.05089676E-9    2.73146488E-10   4.85648281E-10   6.29210199E-9 
  5: 4.12306854E-8    1.11118638E-9    2.87469514E-9    3.72448039E-8 
  6: 1.61223437E-7    3.84349283E-9    1.12767949E-8    1.46103149E-7 
  7: 1.42902491E-6    1.59807452E-8    1.012493E-7      1.31179486E-6 
  8: 1.57929939E-5    9.87789387E-8    1.12454254E-6    1.45696724E-5 
  9: 1.41918596E-4    8.87644758E-7    1.01053353E-5    1.30925616E-4 
 10: 0.15756973       1.01337083E-7    1.12903936E-2    .146279336    

     MMC=[L Cold | B Cold | Y Cold | F Cold]

   MMT: a  10 by 4 matrix

  1: 2.09803077E-11   1.16933226E-12   1.41952208E-12   1.83914534E-11
  2: 1.11889399E-10   5.45935993E-12   7.62606519E-12   9.88039735E-11
  3: 2.0279849E-10    9.18445554E-12   1.38730875E-11   1.79740946E-10
  4: 3.84616671E-10   1.48997633E-11   2.64914424E-11   3.43225466E-10
  5: 9.30071217E-10   2.50658571E-11   6.48466349E-11   8.40158725E-10
  6: 1.83916213E-9    4.38447819E-11   1.28640441E-10   1.6666769E-9  
  7: 5.47552576E-9    6.12326499E-11   3.87952055E-10   5.02634106E-9 
  8: 1.82027985E-8    1.13851314E-10   1.29613305E-9    1.67928141E-8 
  9: 5.45664349E-8    3.41291496E-10   3.88541133E-9    5.0339732E-8  
 10: 1.8182028E-6     1.16933226E-12   1.30280259E-7    1.68792254E-6 

     Scatchard model linear-regession fit estimates
  K1 = 138222119
  S10 = 3.0576542E-10

     Scatchard model euclidean-fit estimates;
  K1 = 287268836
  S10 = 1.82633621E-10

final parameter values
      value               error                dependency    parameter
   1187347821          228157865.8            0.9547016996     K1
   4.947046514e-11     8.636211996e-12        0.9633732247     S10
   0.0804516126        0.001324755233         0.5724726409     C
19 iterations
best weighted sum of squares = 4.378871e-11
weighted root mean square error = 2.501106e-06
weighted deviation fraction = 1.716122e-02
R squared = 9.982448e-01
no active constraints

Note that the log-percent saturation plot using the estimated values for K1 and S10 obtained by linear-regression with the bound/free vs. bound Scatchard-plot form of the data is decidedly inferior to the graph using the direct saturation model-based estimates. Of course, as the following picture shows, this situation is not so clear in the Scatchard plot itself. But the log-percent saturation plot is the more faithful form of the data with regard to minimal error in the x-values, and thus the associated estimated parameters are likely to be the most accurate.

You can see that MLAB is an extremely flexible and general tool for curve-fitting. Moreover, it has a broad range of other useful functions, only a few of which have been alluded to here.