home *** CD-ROM | disk | FTP | other *** search
/ C/C++ Interactive Guide / c-cplusplus-interactive-guide.iso / c_ref / csource4 / 209_01 / ldhfitr.doc < prev    next >
Encoding:
Text File  |  1990-03-05  |  22.4 KB  |  732 lines

  1. LDHFITR.DOC       VERS:- 01.00 DATE:- 09/26/86 TIME:- 10:02:47 PM 
  2.  
  3. description of computer reduction of data from kinetic 
  4.     measurements for the enzyme lactate dehydrogenase 
  5. information on how to run program LDHFIT
  6.  
  7.  
  8.  
  9.  
  10.  
  11.  
  12.         By J. A. Rupley, Tucson, Arizona
  13.  
  14.  
  15.          NOTES ON DATA REDUCTION BY COMPUTER
  16.          AND INFORMATION ON RUNNING THE PROGRAM LDHFIT
  17.     
  18.  
  19.         INTRODUCTION:
  20.  
  21.          In order to obtain conclusions from quantitative measurements,
  22.         there must be some form of data reduction. This can be as simple
  23.         as a comparison by eye of two curves drawn through the data. If,
  24.         however, the data set is large and complex, for example with more
  25.         than one independent variable, and if the questions posed are
  26.         detailed or involve a complicated nonlinear model, then visual or
  27.         graphical methods are less satisfactory than a computer-based
  28.         analysis. Procedures of the latter type are now widely used.
  29.  
  30.          This laboratory is a short introduction to data reduction by use
  31.         of a computer. The intent is to show that a sophisticated
  32.         computer program can be handled easily, that its use saves time
  33.         and effort, that it can treat a more complicated model than can
  34.         be treated graphically, and that it produces information such as
  35.         estimates of uncertainties in the parameters that is difficult or
  36.         impossible to obtain from graphical methods.
  37.  
  38.          The data to be analyzed are initial rate measurements made
  39.         on the lactate dehydrogenase catalyzed reaction of pyruvate with
  40.         NADH, in the presence and absence of lactate as inhibitor. The
  41.         results of the computer fit are the following:  (1) values of the
  42.         kinetic constants V, KmA, KmB, KmAB, KmQ/KmPQ, and KBInhib. The
  43.         first five constants are those that can be evaluated by the
  44.         standard graphical methods of primary and secondary reciprocal
  45.         plots.  The constant KBInhib is the dissociation constant for the
  46.         dead-end complex LDH-NADH-lactate, which is included in the
  47.         mechanism fit by the computer program but cannot be included in
  48.         the mechanism on which the graphical methods are based.  (2)
  49.         Estimates of the standard deviations of the kinetic constants.
  50.         These are needed for an understanding of the reliability and
  51.         significance of the values calculated for the kinetic constants.
  52.         (3) A list of the coordinates of points suitable for construction
  53.         of the lines of the reciprocal plots of the standard graphical
  54.         methods.
  55.  
  56.  
  57.  
  58.  
  59.  
  60.  
  61.  
  62.  
  63.  
  64.  
  65.  
  66.  
  67.  
  68.                                         1
  69.  
  70.  
  71.  
  72.  
  73.  
  74.  
  75.  
  76.  
  77.  
  78.         By J. A. Rupley, Tucson, Arizona
  79.  
  80.         THEORY:
  81.  
  82.         A. REMARKS ON FITTING OF A MODEL TO DATA
  83.  
  84.          In a typical data reduction, a particular model to be tested is
  85.         fit to a set of data points under some criterion for best fit.
  86.         The ith data point of a set of N data points consists of a single
  87.         value for the dependent variable Yobserved(i) measured for
  88.         corresponding single values for the one or more independent
  89.         variables Xobserved(i). The commonly-used least squares criterion
  90.         for quality of fit is the minimum value of the sum of the squares
  91.         of the deviations between the observed values of Y and the values
  92.         of Y calculated according to the model being tested.
  93.  
  94.          Working from the model to be fit to the data, one develops an
  95.         equation relating, for each of the N data points, the dependent
  96.         variable Y to the independent variables X and to a set of M
  97.         variable parameters p:
  98.  
  99.              Ymodel(i) = F(Xobserved(i); p(j), j=1,M)        eq. 1
  100.  
  101.         For example, if the model predicts a linear relationship between
  102.         Y and a single independent variable X :
  103.  
  104.              Ymodel(i) = p(1) + p(2) * Xobserved(i)          eq. 2
  105.  
  106.         The constants p(1) and p(2) of equation (2) are the Y axis
  107.         intercept and the slope, respectively, and of course are the same
  108.         for all data points (for all pairs of values Y(i) and X(i)).
  109.  
  110.          The fitting of a model to data consists of finding the values of
  111.         the M variable parameters p that give the best agreement between
  112.         the N pairs of values of Ymodel(i) and Yobserved(i). Best
  113.         agreement can be defined as the minimum value of the least
  114.         squares function y:
  115.  
  116.                   N
  117.              y = SUM (Yobserved(i) - Ymodel(i))^2 * W(i)     eq. 3
  118.                  i=1
  119.  
  120.  
  121.         The factor W(i) of equation (3) is the normalized reciprocal
  122.         variance (the statistical weighting) of the ith data point, and
  123.         it can be set at unity if the data points are all of equal
  124.         estimated uncertainty.
  125.  
  126.          Combining equations (1) and (3), one sees that the least squares
  127.         function y of equation (3) is a function of the full set of N
  128.         data points and a set of M variable parameters:
  129.  
  130.           y = f(Yobserved(i), Xobserved(i), i=1,N; p(j), j=1,M)   eq. 4
  131.  
  132.  
  133.  
  134.                                         2
  135.  
  136.  
  137.  
  138.  
  139.  
  140.  
  141.  
  142.  
  143.  
  144.          The fitting problem therefore consists of finding the minimum
  145.         value of the least squares function y, which for a given set of
  146.         data depends only on the M variable parameters p (the data points
  147.         Yobserved(i)---Xobserved(i) in equation (4) are constant in the
  148.         fitting). There are several methods commonly used to find the
  149.         minimum of y and thus evaluate the best fit values of the
  150.         parameters p. The more useful of these can handle nonlinear model
  151.         functions F (equation (1)) of arbitrary mathematical form. The
  152.         rate law for lactate dehydrogenase is an example of a nonlinear
  153.         model function.
  154.  
  155.          In the simplex method used here, one constructs an M dimensional
  156.         polyhedron with M + 1 vertices (the simplex). Each dimension of
  157.         the simplex corresponds to a variable parameter of equation (4).
  158.         Each vertex of the simplex is a point in the M dimensional space,
  159.         which is called "parameter space" or "factor space." The M
  160.         coordinates of each vertex are values of the M parameters. Thus
  161.         each vertex of the simplex has an associated value of the least
  162.         squares function y. The starting simplex is constructed to be so
  163.         large as to include within it the point corresponding to the
  164.         minimum value of y. This minimum point has as its coordinates of
  165.         the best fit values of the parameters.
  166.  
  167.          The minimization process shrinks the simplex about the minimum
  168.         point, even though the coordinates of the minimum are not known
  169.         beforehand, until the vertices of the simplex are so close
  170.         together and so nearly equal that an exit test is satisfied. The
  171.         exit test is set so that a desired level of accuracy is obtained.
  172.         The values of the M parameters averaged over all the vertices, ie
  173.         the parameter values for the centroid of the simplex, serve as
  174.         reliable estimates of the best fit parameter values (those for
  175.         the least squares function minimum), because the minimum point is
  176.         known to be inside the shrunken simplex and thus near the
  177.         centroid.
  178.  
  179.          We generally want to estimate the uncertainties in the parameter
  180.         values obtained for a model fit to a particular set of data
  181.         points. Standard deviations of the parameters are calculated by
  182.         the program used here.
  183.  
  184.          There are likely to be large uncertainties in the parameters if
  185.         there are few data points or if there are large deviations
  186.         between Ymodel and Yobserved. As a rule, one should have 5 to 10
  187.         times as many data points as parameters.
  188.  
  189.          The first try at estimating uncertainties of the parameters can
  190.         fail. The calculation involves matrix inversion, the use of
  191.         differences between nearly equal large numbers, and the
  192.         approximation of a complex surface by a simple quadratic
  193.         function. It may be necessary to change certain test values and
  194.         then to repeat the calculation of the standard deviations. In
  195.         particular, if a parameter is close to a bound, so that expansion
  196.         of the simplex in that dimension is not possible, then that
  197.         parameter should be fixed in the quadratic fit.
  198.  
  199.  
  200.                                         3
  201.  
  202.  
  203.  
  204.  
  205.  
  206.  
  207.  
  208.  
  209.  
  210.  
  211.          All fitting methods can fail. We will not discuss problems with
  212.         bounds, local minima, ill-behaved functions, poor quality data,
  213.         physically unreasonable best fits, etc. References given in
  214.         subsection C below discuss fitting methods in more detail.
  215.  
  216.          For additional discussion of fitting by use of the simplex
  217.         method, see "Notes on the fitting program", and the article by
  218.         Nelder and Mead (1965).
  219.  
  220.  
  221.         B. THE FUNCTION FIT TO THE DATA
  222.  
  223.          In the computer program you will use, the function fit to the
  224.         data, which you obtained in the kinetic experiment with lactate
  225.         dehydrogenase, is for an ordered ternary-complex pathway with
  226.         dead-end complexes (EAP and EQB). The following scheme describes
  227.         the pathway:
  228.  
  229.                             EAP
  230.                              |   KPInhib = ea * p/eap
  231.                              |
  232.                   ----------EA-----------
  233.            k1 * a | k-1             k-2 | k2 * b
  234.                   |                     |
  235.                   |                     |
  236.                   E               EAB <---> EPQ                 eq. 5
  237.                   |                     |
  238.                   |                     |
  239.                k4 | k-4 * q     k-3 * p | k3
  240.                   ----------EQ-----------
  241.                              |
  242.                              |   KBInhib = eq * b/eqb
  243.                             EQB
  244.  
  245.          This pathway is more complicated than the one, without dead- end
  246.         complexes, on which are based the standard graphical methods of 
  247.         analysis of two-substrate--two-product reactions:
  248.  
  249.                   -----------EA----------
  250.                   |                     |
  251.                   |                     |
  252.                   E               EAB <---> EPQ                 eq. 6
  253.                   |                     |
  254.                   |                     |
  255.                   -----------EQ----------
  256.  
  257.          For the direction of reaction (pyruvate reduction by NADH) and
  258.         the product inhibitor (lactate) studied in these measurements,
  259.         the rate law for the pathway of equation (5) is:
  260.  
  261.  
  262.  
  263.  
  264.  
  265.  
  266.                                         4
  267.  
  268.  
  269.  
  270.  
  271.  
  272.  
  273.  
  274.  
  275.  
  276.              vo  =  V / [ 1  +  KmQ/KmPQ * p  +  KmA/a       eq. 7
  277.  
  278.                        +  KmB * (1 + KmQ/KmPQ * p) * (1 + 1/KPInhib * p)
  279.  
  280.                        +  KmAB/(a * b) * (1 + KmQ/KmPQ * p)
  281.  
  282.                        +  k3/(k3 + k4) * 1/KBInhib * b ]
  283.  
  284.          The presence in the pathway of equation (5) of the dead end
  285.         complexes EAP and EQB leads to a significantly more complicated
  286.         rate law than is found for the simpler pathway of equation (6);
  287.         compare equation (7) with the following equation, for the "bare"
  288.         compulsory order pathway without dead-end complexes (the pathway
  289.         of equation (6)):
  290.  
  291.              vo  =  V / [ 1  +  KmQ/KmPQ * p  +  KmA/a       eq. 8
  292.  
  293.                        +  KmB * (1 + KmQ/KmPQ * p)
  294.  
  295.                        +  KmAB/(a * b) * (1 + KmQ/KmPQ * p) ]
  296.  
  297.          The computer program used to reduce the kinetic data fits
  298.         equation (7) to measurements of the initial rate v(0) as a
  299.         function of the concentrations of NADH, pyruvate, and lactate.
  300.  
  301.          The parameters that are evaluated in the fitting are listed in
  302.         Table I and are those of equation (7).
  303.  
  304.          Several points should be noted regarding equations (7) and (8):
  305.         (1) The last term of equation (7), containing the equilibrium
  306.         constant KBInhib, probably can be neglected under initial rate
  307.         conditions with q=0; in the fitting program this is recognized by
  308.         setting the last term at a small value, 1E-10, which removes it
  309.         from the equation.  (2) With this change, equations (7) and (8)
  310.         are identical except for the appearance in one term of equation
  311.         (7) of a factor containing the equilibrium constant KPInhib,
  312.         which is for dissociation of the dead-end complex EAP defined in
  313.         the pathway of equation (5).
  314.  
  315.  
  316.  
  317.  
  318.  
  319.  
  320.  
  321.  
  322.  
  323.  
  324.  
  325.  
  326.  
  327.  
  328.  
  329.  
  330.  
  331.  
  332.                                         5
  333.  
  334.  
  335.  
  336.  
  337.  
  338.  
  339.  
  340.  
  341.  
  342.         By J.A. Rupley, Tucson, Arizona
  343.  
  344.         REFERENCES:
  345.  
  346.         A simplex method for function minimization.
  347.         J.A. Nelder and R. Mead (1965. Computer J. 7, 308.
  348.  
  349.         Digital computer user's handbook.
  350.         M. Klerer and G.A. Korn (1967). Mcgraw-Hill, New York.
  351.  
  352.         Data analysis in biochemistry and biophysics.
  353.         M.E. Magar (1972). Academic, New York.
  354.  
  355.         The solution of the general least squares problem with special
  356.         reference to high-speed computers.
  357.         R.H. Moore and R.K. Ziegler (1960). Los Alamos Scientific
  358.         Laboratory Report LA-2367.
  359.  
  360.  
  361.  
  362.  
  363.  
  364.  
  365.  
  366.  
  367.  
  368.  
  369.  
  370.  
  371.  
  372.  
  373.  
  374.  
  375.  
  376.  
  377.  
  378.  
  379.  
  380.  
  381.  
  382.  
  383.  
  384.  
  385.  
  386.  
  387.  
  388.  
  389.  
  390.  
  391.  
  392.  
  393.  
  394.  
  395.  
  396.  
  397.  
  398.                                         6
  399.  
  400.  
  401.  
  402.  
  403.  
  404.  
  405.  
  406.  
  407.  
  408.         By J.A. Rupley, Tucson, Arizona
  409.  
  410.         TABLE I:
  411.  
  412.                                  EQUATION (7)
  413.         FITTING                  RATE LAW
  414.         PARAMETERS               PARAMETERS
  415.  
  416.         1                        V
  417.  
  418.         2                        KmA  =  KmNADH
  419.  
  420.         3                        KmB  =  KmPyruvate
  421.  
  422.         4                        KmAB =  KmNADH-Pyruvate
  423.  
  424.         5                        KmQ/KmPQ = KmNAD/KmLactate-NAD
  425.  
  426.         6                        1/KPInhib = 1/KInhibLactate
  427.  
  428.         7                        k3/(k3 + k4) * 1/KInhibPyruvate
  429.                                       parm(7) approx. equal to 0 at t=0
  430.  
  431.  
  432.         INDEPENDENT VARIABLES:
  433.  
  434.         a = [NADH]     b = [Pyruvate]    p = [Lactate]    q = [NAD]
  435.                                                           q = 0 at t=0
  436.  
  437.  
  438.         DEPENDENT VARIABLE:
  439.  
  440.         vo = initial rate of conversion of pyruvate to lactate
  441.  
  442.  
  443.  
  444.  
  445.  
  446.  
  447.  
  448.  
  449.  
  450.  
  451.  
  452.  
  453.  
  454.  
  455.  
  456.  
  457.  
  458.  
  459.  
  460.  
  461.  
  462.  
  463.  
  464.                                         7
  465.  
  466.  
  467.  
  468.  
  469.  
  470.  
  471.  
  472.  
  473.  
  474.         By J. A. Rupley, Tucson, Arizona
  475.  
  476.         PROCEDURES:
  477.  
  478.         A. SUMMARY OF STEPS IN THE DATA REDUCTION
  479.  
  480.          (1) Startup. Create a new file with your data, following the
  481.         template of the DATA SHEET. Setup the title and the starting
  482.         ranges of the parameters, from which the starting simplex is
  483.         constructed. Be sure to save the setup on disk.
  484.  
  485.          (2) Minimization. The minimum of the least squares function is
  486.         located by an iterative method, in which the simplex, a
  487.         polyhedron in parameter space, is shrunk about the point
  488.         representing the minimum. The starting simplex is large,
  489.         sufficiently so to certainly include the minimum point. After
  490.         shrinking, all the vertices of the simplex are so close together,
  491.         and thus so close to the minimum point, that the average of the
  492.         vertices, the centroid, is a good estimate of the minimum point.
  493.         The minimization can take more than two hours. These results
  494.         should be saved on a disk file, so if the fitting fails, it is
  495.         possible to quickly return to this point.
  496.  
  497.          (3) Quadratic fit and standard deviations. The least squares
  498.         function about the minimum is approximated by a quadratic
  499.         surface, the shape of which gives estimates of the standard
  500.         deviations of the parameters. This step takes about 10 minutes.
  501.         There is considerable output, the last of which is a list of the
  502.         standard deviations.
  503.  
  504.          The fitting program can be set to cycle between a set of
  505.         minimization iterations and the quadratic approximation. This
  506.         option can produce quicker convergence.
  507.  
  508.  
  509.         B. FAILURE OF THE FITTING
  510.  
  511.          Fitting can lead to physically unreasonable values of the
  512.         parameters, for example for a kinetics experiment a set of
  513.         parameters for which V is unreasonably large or small, or a set
  514.         with standard deviations such that one or more of the principal
  515.         parameters (V, KmA, KmB, KmAB) have no significance.
  516.  
  517.          The best fit may be unacceptable. To judge whether a fit is
  518.         acceptable, compare the root mean square error (RMS error) with
  519.         the expected uncertainty in your independent variable, the meas-
  520.         ured rates.
  521.  
  522.          The most likely causes of failure are:
  523.  
  524.          (1) A few data points (flyers) that are widely wrong. The
  525.         results should always be examined for flyers, which we will
  526.         somewhat arbitrarily define as measured values differing from
  527.         calculated values by more than 3 standard deviations: ABS(Ymodel
  528.  
  529.  
  530.                                         8
  531.  
  532.  
  533.  
  534.  
  535.  
  536.  
  537.  
  538.  
  539.  
  540.         - Yobserved) > 3 x RMS ERROR. If there are flyers, delete them
  541.         and repeat the calculation, or correct the values by
  542.         recalculation of vo from your experimental data.
  543.  
  544.          (2) Trapping of the fit in a local minimum in parameter space
  545.         not near the true minimum, or a broad and ill-defined minimum
  546.         region. To test for trapping or a broad minimum, start the
  547.         fitting with a different simplex. For example, set tighter
  548.         starting ranges on one or more of the principal parameters (V,
  549.         KmA, KmB, KmAB). Alternatively, set one principal parameter, e.g.
  550.         V, at a reasonable value.
  551.  
  552.          Generally, to test for trapping, one repeats the fitting with a
  553.         more expanded starting simplex, not with a tighter one as
  554.         suggested above. However, assuming the tighter ranges include
  555.         physically reasonable values of the parameters, the approach
  556.         above, using a restricted simplex, has the advantage of testing
  557.         whether a set of physically reasonable values can give an
  558.         acceptable fit of the data. To judge whether a fit is acceptable,
  559.         even though it may not be as good as a physically unreasonable
  560.         best fit, compare the root mean square error (RMS error) with the
  561.         expected uncertainty in your dependent variable, the measured
  562.         rates.
  563.  
  564.          (3) Systematic error in the measurements or poor data. Look for
  565.         systematic bias in the extraction of the slope. Look for errors
  566.         in calculating the value of vo.
  567.  
  568.          One can plot the measured values with the calculated lines, ie
  569.         construct the reciprocal plots. It should then be clear what the
  570.         problem is. Correction, if it is possible, generally requires
  571.         examination of the original measurements, ie the tracings of
  572.         transmittance versus time. A point to be made is that the
  573.         computer fitting, because it is unbiased, is more likely than the
  574.         experimenter to detect systematic error or poor data.
  575.  
  576.  
  577.  
  578.  
  579.  
  580.  
  581.  
  582.  
  583.  
  584.  
  585.  
  586.  
  587.  
  588.  
  589.  
  590.  
  591.  
  592.  
  593.  
  594.  
  595.  
  596.                                         9
  597.  
  598.  
  599.  
  600.  
  601.  
  602.  
  603.  
  604.  
  605.  
  606.         By J. A. Rupley, Tucson, Arizona
  607.  
  608.  
  609.         INSTRUCTIONS FOR FITTING KINETIC DATA FOR LACTATE DEHYDROGENASE
  610.         WITH THE PROGRAM LDHFIT.COM (C-LANGUAGE):
  611.  
  612.  
  613.         A. CREATE THE DATA FILE
  614.  
  615.          By use of a word processing program, modify the template file
  616.         LDHDATA, incorporating your vo values and so constructing your
  617.         own data file, eg TUE-C.DAT.
  618.  
  619.  
  620.         B. RUN THE FITTING PROGRAM
  621.  
  622.          Run the program LDHFIT with your data file (eg TUE-C.DAT) to
  623.         produce the results file (eg TUE-C.PRT)
  624.  
  625.          A>LDHFIT TUE-C.DAT TUE-C.PRT
  626.  
  627.          It takes about 50 minutes for 100 iterations. 200 to 400
  628.         iterations are commonly required for convergence. The fitting
  629.         will cycle between sessions of Nelder-Mead minimization (30
  630.         iterations) and quadratic fitting.
  631.  
  632.          You should watch the fitting. Let the program run freely until
  633.         the test value is < 10-4. Then:  (1) See if a parameter needs to
  634.         be eliminated. If the calculation of -pmin- fails at this stage
  635.         in the fitting, probably the fitting will not converge. If one of
  636.         the parameter values listed under -pmin- in the summary of the
  637.         quadratic fit results is negative, then this parameter should be
  638.         fixed at 1.e-10. This can be done by use of a program option (see
  639.         below).  (2) Examine the data listing for flyers: (Yobserved(i) -
  640.         Ycalculated(i)) > 3 * rms error. If these are present, you should
  641.         delete them from the data set, or set their weight to zero. This
  642.         is done by reworking the data file, with a word processor, then
  643.         rerunning LDHFIT. You will get faster convergence by setting the
  644.         starting simplex in the data file at an intermediate point in the
  645.         fitting.
  646.  
  647.          To exercise options, during the minimization type:  CTRL-X
  648.  
  649.          After some output, you will be asked to select an option:
  650.              display of the data and then quadratic fit (<CR>)
  651.              fix of a parameter (CTRL-F)
  652.              exit from the program to operating system (CTRL-C)
  653.              return to the fitting (CTRL-X)
  654.  
  655.  
  656.  
  657.  
  658.  
  659.  
  660.  
  661.  
  662.                                         10
  663.  
  664.  
  665.  
  666.  
  667.  
  668.  
  669.  
  670.  
  671.  
  672.         B. AFTER THE FITTING CONVERGES
  673.  
  674.          When the simplex has shrunk to where it satisfies the exit test,
  675.         the minimization pauses and expects input from the operator:
  676.  
  677.              respond with <RETURN> when asked if you want to  display data;
  678.  
  679.              respond with <RETURN> when asked at the next pause if you
  680.         want to carry out quadratic fit;
  681.  
  682.              at the next pause, record the number of the iteration, and
  683.         exit by typing CTRL-C;
  684.  
  685.  
  686.         D. CONDENSE THE RESULTS FILE
  687.  
  688.          The results file (eg TUE-C.PRT) will be too long to print.
  689.         Select the information from the last data display and quadratic
  690.         fitting by use of a word processor or a special search program.
  691.  
  692.  
  693.         E. FAILURE OF THE FIT
  694.  
  695.          Repeat the fit after checking your data and possibly resetting
  696.         the starting parameter values.
  697.  
  698.  
  699.  
  700.  
  701.  
  702.  
  703.  
  704.  
  705.  
  706.  
  707.  
  708.  
  709.  
  710.  
  711.  
  712.  
  713.  
  714.  
  715.  
  716.  
  717.  
  718.  
  719.  
  720.  
  721.  
  722.  
  723.  
  724.  
  725.  
  726.  
  727.  
  728.                                         11
  729.  
  730.  
  731.  
  732.